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
81
82
83


#pragma once

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

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

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

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

  void add(double x) {
    if (n > wait) {
      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 - wait); }

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

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

    std::vector<double> ρtmp;

    for (double Ct : C) {
      ρtmp.push_back((Ct / (n - wait) - 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 - wait)};
  }

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

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

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