Central Limit Theorem

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct Coverage {
5
    double p1, p2, p3; // fractions within |Z|<=1,2,3
6
    double meanZ, varZ;
7
};
8
9
// Compute coverage of standardized sums for a chosen distribution
10
Coverage simulateCLT(size_t n, size_t M, const string& distName, mt19937_64& rng) {
11
    // Choose distribution and its true mean/variance
12
    function<double()> gen;
13
    double mu = 0.0, sigma = 1.0;
14
15
    if (distName == "uniform01") {
16
        uniform_real_distribution<double> U(0.0, 1.0);
17
        gen = [&](){ return U(rng); };
18
        mu = 0.5; sigma = sqrt(1.0/12.0);
19
    } else if (distName == "exponential1") {
20
        exponential_distribution<double> E(1.0);
21
        gen = [&](){ return E(rng); };
22
        mu = 1.0; sigma = 1.0;
23
    } else if (distName == "bernoulli03") {
24
        bernoulli_distribution B(0.3);
25
        gen = [&](){ return B(rng) ? 1.0 : 0.0; };
26
        mu = 0.3; sigma = sqrt(0.3 * 0.7);
27
    } else {
28
        throw runtime_error("Unknown distribution");
29
    }
30
31
    size_t c1 = 0, c2 = 0, c3 = 0;
32
    long double sumZ = 0.0L, sumZ2 = 0.0L;
33
34
    for (size_t m = 0; m < M; ++m) {
35
        long double S = 0.0L;
36
        for (size_t i = 0; i < n; ++i) S += gen();
37
        long double Z = (S - n*mu) / (sigma * sqrt((long double)n));
38
        sumZ += Z; sumZ2 += Z*Z;
39
        if (fabsl(Z) <= 1.0L) ++c1;
40
        if (fabsl(Z) <= 2.0L) ++c2;
41
        if (fabsl(Z) <= 3.0L) ++c3;
42
    }
43
44
    Coverage res;
45
    res.p1 = (double)c1 / (double)M;
46
    res.p2 = (double)c2 / (double)M;
47
    res.p3 = (double)c3 / (double)M;
48
    long double meanZ = sumZ / (long double)M;
49
    long double varZ = sumZ2 / (long double)M - meanZ*meanZ;
50
    res.meanZ = (double)meanZ; res.varZ = (double)varZ;
51
    return res;
52
}
53
54
int main() {
55
    ios::sync_with_stdio(false);
56
    cin.tie(nullptr);
57
58
    mt19937_64 rng(1234567);
59
    size_t n = 50;      // sample size per sum (try 10, 50, 200)
60
    size_t M = 20000;   // number of replications (increase for smoother results)
61
62
    vector<string> dists = {"uniform01", "exponential1", "bernoulli03"};
63
64
    cout.setf(ios::fixed); cout << setprecision(4);
65
    cout << "Normal theory: within 1,2,3 sigmas ≈ 0.6827, 0.9545, 0.9973\n\n";
66
67
    for (const auto& name : dists) {
68
        Coverage cv = simulateCLT(n, M, name, rng);
69
        cout << "Distribution: " << name << " (n=" << n << ", M=" << M << ")\n";
70
        cout << "  Empirical coverage |Z|<=1: " << cv.p1 << ", <=2: " << cv.p2 << ", <=3: " << cv.p3 << "\n";
71
        cout << "  Empirical mean(Z): " << cv.meanZ << ", var(Z): " << cv.varZ << "\n\n";
72
    }
73
74
    return 0;
75
}
76

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct OnlineStats {
5
    long long n = 0;
6
    long double mean = 0.0L;
7
    long double M2 = 0.0L; // sum of squared deviations
8
    void add(double x) {
9
        n++;
10
        long double delta = x - mean;
11
        mean += delta / n;
12
        long double delta2 = x - mean;
13
        M2 += delta * delta2; // stable update
14
    }
15
    long double variance() const { return n > 1 ? M2 / (n - 1) : NAN; }
16
};
17
18
int main(){
19
    ios::sync_with_stdio(false);
20
    cin.tie(nullptr);
21
22
    // Example: generate synthetic skewed data (Exponential(1))
23
    mt19937_64 rng(42);
24
    exponential_distribution<double> E(1.0);
25
26
    OnlineStats stats;
27
    const size_t N = 5000; // try varying N to see CI width shrink like 1/sqrt(N)
28
    for (size_t i = 0; i < N; ++i) stats.add(E(rng));
29
30
    long double xbar = stats.mean;
31
    long double s2 = stats.variance();
32
    long double s = sqrt(s2);
33
34
    // 95% z-interval using CLT (asymptotic). z_{0.975} ≈ 1.96
35
    const long double z975 = 1.959963984540054; 
36
    long double half_width = z975 * s / sqrt((long double)stats.n);
37
38
    cout.setf(ios::fixed); cout << setprecision(6);
39
    cout << "n = " << stats.n << "\n";
40
    cout << "Sample mean xbar = " << (double)xbar << "\n";
41
    cout << "Sample std s = " << (double)s << "\n";
42
    cout << "95% CI for mean (CLT): [" << (double)(xbar - half_width) << ", " << (double)(xbar + half_width) << "]\n";
43
44
    // For Exponential(1), true mean = 1; check if interval contains 1
45
    cout << "Contains true mean 1.0? " << ((xbar - half_width <= 1.0 && 1.0 <= xbar + half_width) ? "yes" : "no") << "\n";
46
47
    return 0;
48
}
49

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Simulate 1D linear regression gradient noise at a fixed theta
5
// Data: (x_i, y_i) with y_i = a*x_i + noise. Gradient per-sample: g_i = (theta*x_i - y_i)*x_i
6
// We compare the batch gradient to the full gradient and examine standardized noise.
7
8
int main(){
9
    ios::sync_with_stdio(false);
10
    cin.tie(nullptr);
11
12
    mt19937_64 rng(2024);
13
    const size_t N = 20000; // dataset size
14
    normal_distribution<double> X(0.0, 1.0);
15
    normal_distribution<double> eps(0.0, 0.5); // observation noise
16
17
    // True underlying relation y = a*x + eps
18
    const double a = 1.5;
19
20
    vector<double> x(N), y(N);
21
    for(size_t i=0;i<N;++i){
22
        x[i] = X(rng);
23
        y[i] = a * x[i] + eps(rng);
24
    }
25
26
    // Pick a fixed parameter theta (not necessarily optimal)
27
    const double theta = 0.8;
28
29
    // Compute per-example gradients and the full gradient
30
    vector<double> g(N);
31
    long double g_full = 0.0L;
32
    for(size_t i=0;i<N;++i){
33
        double gi = (theta * x[i] - y[i]) * x[i];
34
        g[i] = gi;
35
        g_full += gi;
36
    }
37
    g_full /= (long double)N; // mean gradient over full data
38
39
    // Repeatedly sample mini-batches; record standardized gradient noise
40
    const size_t B = 64;           // batch size
41
    const size_t M = 20000;        // number of mini-batches to sample
42
43
    uniform_int_distribution<size_t> U(0, N-1);
44
45
    // Estimate per-example gradient variance (needed to standardize)
46
    long double mu_g = g_full;
47
    long double m2 = 0.0L;
48
    for(size_t i=0;i<N;++i){
49
        long double d = g[i] - mu_g;
50
        m2 += d*d;
51
    }
52
    long double var_g = m2 / (N - 1);
53
54
    size_t c1=0,c2=0,c3=0; // coverage counts for |Z|<=1,2,3
55
    long double sumZ=0.0L, sumZ2=0.0L;
56
57
    for(size_t m=0;m<M;++m){
58
        long double g_batch = 0.0L;
59
        for(size_t b=0;b<B;++b){
60
            size_t idx = U(rng);
61
            g_batch += g[idx];
62
        }
63
        g_batch /= (long double)B;
64
        // Gradient noise (batch - full). Standard error ≈ sqrt(var_g / B)
65
        long double se = sqrt(var_g / (long double)B);
66
        long double Z = (g_batch - g_full) / se;
67
        sumZ += Z; sumZ2 += Z*Z;
68
        if (fabsl(Z) <= 1.0L) ++c1;
69
        if (fabsl(Z) <= 2.0L) ++c2;
70
        if (fabsl(Z) <= 3.0L) ++c3;
71
    }
72
73
    cout.setf(ios::fixed); cout << setprecision(4);
74
    cout << "Mini-batch size B = " << B << ", replications M = " << M << "\n";
75
    cout << "Coverage |Z|<=1: " << (double)c1/M << ", <=2: " << (double)c2/M << ", <=3: " << (double)c3/M << "\n";
76
    long double meanZ = sumZ / (long double)M;
77
    long double varZ = sumZ2 / (long double)M - meanZ*meanZ;
78
    cout << "Empirical mean(Z): " << (double)meanZ << ", var(Z): " << (double)varZ << " (should be near 0 and 1)\n";
79
80
    return 0;
81
}
82