Sufficient Statistics

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute the Bernoulli log-likelihood in two ways:
5
// (1) From raw data (O(n) per evaluation)
6
// (2) From sufficient statistic T = sum x_i (O(1) per evaluation)
7
// Demonstrate they are identical for any p in (0,1).
8
9
double loglik_raw(const vector<int>& x, double p) {
10
    // Guard for p in (0,1)
11
    if (p <= 0.0 || p >= 1.0) return -INFINITY;
12
    double ll = 0.0;
13
    for (int xi : x) {
14
        // xi must be 0 or 1
15
        if (xi == 1) ll += log(p);
16
        else if (xi == 0) ll += log(1.0 - p);
17
        else return -INFINITY; // invalid Bernoulli observation
18
    }
19
    return ll;
20
}
21
22
double loglik_from_T(int T, int n, double p) {
23
    if (p <= 0.0 || p >= 1.0) return -INFINITY;
24
    // L(p|x) = p^T (1-p)^{n-T} => log L = T log p + (n-T) log(1-p)
25
    return T * log(p) + (n - T) * log(1.0 - p);
26
}
27
28
int main() {
29
    ios::sync_with_stdio(false);
30
    cin.tie(nullptr);
31
32
    // Example data: 10 trials with 6 successes (order arbitrary)
33
    vector<int> x = {1,0,1,1,0,1,0,1,0,1};
34
    int n = (int)x.size();
35
    int T = accumulate(x.begin(), x.end(), 0);
36
37
    cout << "n=" << n << ", T=sum x_i=" << T << "\n";
38
39
    // Compare log-likelihoods over a grid of p values
40
    vector<double> grid;
41
    for (int i = 1; i <= 9; ++i) grid.push_back(i / 10.0); // 0.1..0.9
42
43
    cout.setf(std::ios::fixed); cout << setprecision(6);
44
    for (double p : grid) {
45
        double ll1 = loglik_raw(x, p);
46
        double ll2 = loglik_from_T(T, n, p);
47
        cout << "p=" << p << ": ll_raw=" << ll1 << ", ll_T=" << ll2
48
             << ", diff=" << fabs(ll1 - ll2) << "\n";
49
    }
50
51
    // MLE from T: p_hat = T/n
52
    double p_hat = (double)T / n;
53
    cout << "MLE p_hat = T/n = " << p_hat << "\n";
54
    return 0;
55
}
56

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// We assume X_i ~ N(mu, sigma2) with known sigma2.
5
// Sufficient statistic for mu is the sum S = sum X_i (or equivalently the mean).
6
// This aggregator supports streaming adds and merges for distributed computation.
7
8
struct NormalKnownVarAgg {
9
    long long n = 0;     // number of observations
10
    long double S = 0.0; // sum of observations
11
12
    void add(long double x) {
13
        ++n; S += x; // For better FP accuracy, one could use Kahan summation.
14
    }
15
16
    void merge(const NormalKnownVarAgg& other) {
17
        n += other.n; S += other.S;
18
    }
19
20
    long double mean() const { return (n == 0) ? NAN : (S / n); }
21
22
    // Log-likelihood using only S and n:
23
    long double loglik(long double mu, long double sigma2, long double sumsq_residual_part) const {
24
        // log L(mu|x) = -n/2 log(2pi sigma2) - (1/(2 sigma2)) sum (x_i - mu)^2
25
        // Expand: sum (x_i - mu)^2 = sum x_i^2 - 2 mu S + n mu^2
26
        // If we precompute R = sum x_i^2 (theta-free), then we only need S and n for mu dependence.
27
        long double R = sumsq_residual_part; // equals sum x_i^2
28
        long double term = R - 2.0L * mu * S + (long double)n * mu * mu;
29
        return -0.5L * n * log(2.0L * M_PIl * sigma2) - (1.0L / (2.0L * sigma2)) * term;
30
    }
31
};
32
33
int main(){
34
    ios::sync_with_stdio(false);
35
    cin.tie(nullptr);
36
37
    // Simulate or provide data
38
    vector<long double> x = {2.1L, 1.9L, 2.3L, 2.0L, 1.8L};
39
    long double sigma2 = 0.04L; // known variance (e.g., sigma=0.2)
40
41
    NormalKnownVarAgg agg;
42
    long double R = 0.0L; // sum of squares
43
    for (auto v : x) { agg.add(v); R += v * v; }
44
45
    cout.setf(std::ios::fixed); cout << setprecision(6);
46
    cout << "n=" << agg.n << ", S=" << (double)agg.S << ", mean=" << (double)agg.mean() << "\n";
47
48
    // Compute MLE for mu: mu_hat = S/n
49
    long double mu_hat = agg.mean();
50
    cout << "MLE mu_hat = S/n = " << (double)mu_hat << "\n";
51
52
    // Compare log-likelihood at several mu values using only (n,S,R)
53
    for (long double mu : {1.8L, 1.9L, 2.0L, 2.1L, 2.2L}) {
54
        long double ll = agg.loglik(mu, sigma2, R);
55
        cout << "mu=" << (double)mu << ", logL=" << (double)ll << "\n";
56
    }
57
58
    return 0;
59
}
60

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// For iid Poisson(lambda), T = sum X_i is sufficient for lambda.
5
// We implement a streaming aggregator that can be merged across shards.
6
7
struct PoissonAgg {
8
    long long n = 0;       // number of observations
9
    long double S = 0.0L;  // sum of counts
10
11
    void add(long long x) {
12
        if (x < 0) throw runtime_error("Poisson count must be nonnegative");
13
        ++n; S += (long double)x;
14
    }
15
16
    void merge(const PoissonAgg& other) {
17
        n += other.n; S += other.S;
18
    }
19
20
    long double mle_lambda() const { return (n == 0) ? NAN : (S / n); }
21
22
    long double loglik(long double lambda, const vector<int>& x) const {
23
        if (lambda <= 0.0L) return -INFINITY;
24
        // Raw-data log-likelihood: sum( -lambda + x_i log lambda - log x_i! )
25
        long double ll = 0.0L;
26
        for (int xi : x) {
27
            if (xi < 0) return -INFINITY;
28
            ll += -lambda + xi * log(lambda) - lgammal((long double)xi + 1.0L); // log(x!) via lgamma
29
        }
30
        return ll;
31
    }
32
33
    long double loglik_from_T(long double lambda) const {
34
        if (lambda <= 0.0L) return -INFINITY;
35
        // Using factorization: log L = -n lambda + S log lambda - sum log x_i!
36
        // The last term is theta-free; we omit it when comparing across lambda.
37
        return -n * lambda + S * log(lambda);
38
    }
39
};
40
41
int main(){
42
    ios::sync_with_stdio(false);
43
    cin.tie(nullptr);
44
45
    vector<int> x = {2, 0, 1, 3, 2, 4, 1, 0};
46
47
    PoissonAgg a, b;
48
    // Simulate two shards
49
    for (size_t i = 0; i < x.size(); ++i) {
50
        if (i < x.size()/2) a.add(x[i]); else b.add(x[i]);
51
    }
52
53
    // Merge summaries
54
    PoissonAgg all = a; all.merge(b);
55
56
    cout.setf(std::ios::fixed); cout << setprecision(6);
57
    cout << "n=" << all.n << ", S=" << (double)all.S << ", lambda_hat=" << (double)all.mle_lambda() << "\n";
58
59
    // Compare raw vs. sufficient-stat log-likelihood up to a constant
60
    for (long double lambda : {0.5L, 1.0L, 1.5L, 2.0L}) {
61
        long double ll_raw = all.loglik(lambda, x);
62
        long double ll_T = all.loglik_from_T(lambda);
63
        cout << "lambda=" << (double)lambda
64
             << ": ll_raw (up to const)=" << (double)ll_raw
65
             << ", ll_T (no const)=" << (double)ll_T << "\n";
66
    }
67
68
    return 0;
69
}
70