quantity.hpp


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80


#pragma once

#include <array>
#include <cmath>
#include <list>
#include <vector>

class Quantity {
private:
  double total;
  double total2;
  std::vector<double> C;

public:
  unsigned n;
  std::list<double> hist;

  Quantity(unsigned lag) : C(lag) {
    n = 0;
    total = 0;
    total2 = 0;
  }

  void add(double x) {
    hist.push_front(x);
    if (hist.size() > C.size()) {
      hist.pop_back();
      unsigned t = 0;
      for (double a : hist) {
        C[t] += a * x;
        t++;
      }
      total += x;
      total2 += pow(x, 2);
      n++;
    }
  }

  double avg() const { return total / n; }

  double avg2() const { return total2 / n; }

  std::vector<double> ρ() const {
    double C0 = C.front() / n;
    double avg2 = pow(total / n, 2);

    std::vector<double> ρtmp;

    for (double Ct : C) {
      ρtmp.push_back((Ct / n - avg2) / (C0 - avg2));
    }

    return ρtmp;
  }

  std::array<double, 2> τ() const {
    double τtmp = 0.5;
    unsigned M = 1;
    double c = 8.0;

    std::vector<double> ρ_tmp = this->ρ();

    while (c * τtmp > M && M < C.size()) {
      τtmp += ρ_tmp[M];
      M++;
    }

    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
  }

  double σ() const {
    return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
  }

  double serr() const { return sqrt(this->σ()); }

  unsigned num_added() const { return n; }
};