Expectation, Variance & Moments

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <limits>
5
6
// A class that computes mean and variance in one pass using Welford's algorithm.
7
// It also exposes population and unbiased sample variance.
8
class OnlineStats {
9
public:
10
    OnlineStats() : n_(0), mean_(0.0), M2_(0.0) {}
11
12
    // Add a new observation x
13
    void add(double x) {
14
        ++n_;
15
        double delta = x - mean_;
16
        mean_ += delta / static_cast<double>(n_);
17
        double delta2 = x - mean_;
18
        M2_ += delta * delta2; // Accumulate sum of squared deviations
19
    }
20
21
    // Merge another OnlineStats into this one (Chan's method)
22
    void merge(const OnlineStats& other) {
23
        if (other.n_ == 0) return;
24
        if (n_ == 0) { *this = other; return; }
25
        double n_total = static_cast<double>(n_ + other.n_);
26
        double delta = other.mean_ - mean_;
27
        M2_ = M2_ + other.M2_ + delta * delta * (static_cast<double>(n_) * static_cast<double>(other.n_) / n_total);
28
        mean_ = (static_cast<double>(n_) * mean_ + static_cast<double>(other.n_) * other.mean_) / n_total;
29
        n_ += other.n_;
30
    }
31
32
    // Number of observations
33
    std::size_t count() const { return n_; }
34
35
    // Mean of the observations
36
    double mean() const { return mean_; }
37
38
    // Population variance: Var = M2 / n (defined if n > 0)
39
    double variance_population() const { return (n_ > 0) ? (M2_ / static_cast<double>(n_)) : std::numeric_limits<double>::quiet_NaN(); }
40
41
    // Unbiased sample variance: s^2 = M2 / (n - 1) (defined if n > 1)
42
    double variance_unbiased() const { return (n_ > 1) ? (M2_ / static_cast<double>(n_ - 1)) : std::numeric_limits<double>::quiet_NaN(); }
43
44
private:
45
    std::size_t n_;    // count
46
    double mean_;      // running mean
47
    double M2_;        // sum of squared deviations from the mean
48
};
49
50
int main() {
51
    // Example data: imagine these are i.i.d. measurements
52
    std::vector<double> data = {2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0};
53
54
    OnlineStats stats;
55
    for (double x : data) stats.add(x);
56
57
    std::cout << "Count: " << stats.count() << "\n";
58
    std::cout << "Mean:  " << stats.mean() << "\n";
59
    std::cout << "Var (population): " << stats.variance_population() << "\n";
60
    std::cout << "Var (unbiased sample): " << stats.variance_unbiased() << "\n";
61
62
    // Verify Var = E[X^2] - (E[X])^2 using the same data (two-pass for demonstration only)
63
    double mean = 0.0, mean2 = 0.0;
64
    for (double x : data) { mean += x; mean2 += x * x; }
65
    mean /= data.size();
66
    mean2 /= data.size();
67
    double var_identity = mean2 - mean * mean;
68
69
    std::cout << "Var via identity E[X^2] - (E[X])^2: " << var_identity << "\n";
70
    return 0;
71
}
72

1
#include <iostream>
2
#include <random>
3
#include <cmath>
4
#include <limits>
5
6
// Estimate expectations by averaging many random samples.
7
int main() {
8
    std::mt19937_64 rng(123456789ULL);               // Fixed seed for reproducibility
9
    std::normal_distribution<double> dist(0.0, 1.0); // Standard normal N(0,1)
10
11
    const std::size_t n = 1000000; // Number of samples (1e6)
12
13
    long double sum = 0.0L;        // For E[X]
14
    long double sumsq = 0.0L;      // For E[X^2]
15
    long double sum4 = 0.0L;       // For E[X^4] (should be 3 for N(0,1))
16
17
    for (std::size_t i = 0; i < n; ++i) {
18
        double x = dist(rng);
19
        sum += x;
20
        sumsq += x * x;
21
        double x2 = x * x;
22
        sum4 += x2 * x2;
23
    }
24
25
    long double EX = sum / n;
26
    long double EX2 = sumsq / n;
27
    long double EX4 = sum4 / n;
28
    long double VarX = EX2 - EX * EX; // Identity Var = E[X^2] - (E[X])^2
29
30
    // Standard error estimates for the mean and for E[X^2]
31
    long double s2_for_mean = VarX;              // Var(X)
32
    long double se_mean = std::sqrt(s2_for_mean / n);
33
34
    // For E[X^2], variance is Var(X^2)/n; estimate Var(X^2) using E[X^4] - (E[X^2])^2
35
    long double VarX2 = EX4 - EX2 * EX2;
36
    long double se_EX2 = std::sqrt(VarX2 / n);
37
38
    std::cout.setf(std::ios::fixed); std::cout.precision(6);
39
    std::cout << "Monte Carlo with n = " << n << "\n";
40
    std::cout << "E[X]        ≈ " << static_cast<double>(EX)  << " (theory: 0)\n";
41
    std::cout << "E[X^2]      ≈ " << static_cast<double>(EX2) << " (theory: 1)\n";
42
    std::cout << "Var(X)      ≈ " << static_cast<double>(VarX) << " (theory: 1)\n";
43
    std::cout << "E[X^4]      ≈ " << static_cast<double>(EX4) << " (theory: 3)\n";
44
    std::cout << "SE(E[X])    ≈ " << static_cast<double>(se_mean) << "\n";
45
    std::cout << "SE(E[X^2])  ≈ " << static_cast<double>(se_EX2) << "\n";
46
    return 0;
47
}
48

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
5
struct Stats {
6
    std::size_t n = 0;
7
    long double mean = 0.0L;
8
    long double M2 = 0.0L; // sum of squared deviations
9
};
10
11
// Update with a new value (Welford step)
12
void update(Stats &s, long double x) {
13
    s.n += 1;
14
    long double delta = x - s.mean;
15
    s.mean += delta / s.n;
16
    long double delta2 = x - s.mean;
17
    s.M2 += delta * delta2;
18
}
19
20
// Merge two partial stats (Chan's formulas)
21
Stats merge(const Stats &a, const Stats &b) {
22
    if (a.n == 0) return b;
23
    if (b.n == 0) return a;
24
    Stats out;
25
    out.n = a.n + b.n;
26
    long double delta = b.mean - a.mean;
27
    out.mean = (a.n * a.mean + b.n * b.mean) / out.n;
28
    out.M2 = a.M2 + b.M2 + delta * delta * (static_cast<long double>(a.n) * static_cast<long double>(b.n) / out.n);
29
    return out;
30
}
31
32
int main() {
33
    std::vector<long double> data;
34
    for (int i = 1; i <= 100; ++i) data.push_back(i); // data = 1..100
35
36
    // Split into two shards and compute partial stats
37
    Stats A, B;
38
    for (std::size_t i = 0; i < data.size(); ++i) {
39
        if (i % 2 == 0) update(A, data[i]); else update(B, data[i]);
40
    }
41
42
    // Merge and compare with single-pass
43
    Stats merged = merge(A, B);
44
45
    Stats single;
46
    for (auto x : data) update(single, x);
47
48
    auto pop_var = [](const Stats &s){ return s.M2 / (s.n); };
49
    auto samp_var = [](const Stats &s){ return s.n > 1 ? s.M2 / (s.n - 1) : std::nanl(""); };
50
51
    std::cout.setf(std::ios::fixed); std::cout.precision(6);
52
    std::cout << "Merged mean:         " << static_cast<double>(merged.mean) << "\n";
53
    std::cout << "Single-pass mean:    " << static_cast<double>(single.mean) << "\n";
54
    std::cout << "Merged var (pop):    " << static_cast<double>(pop_var(merged)) << "\n";
55
    std::cout << "Single var (pop):    " << static_cast<double>(pop_var(single)) << "\n";
56
    std::cout << "Merged var (sample): " << static_cast<double>(samp_var(merged)) << "\n";
57
    std::cout << "Single var (sample): " << static_cast<double>(samp_var(single)) << "\n";
58
    return 0;
59
}
60