Minimum Description Length (MDL)

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Solve linear system A x = b via Gaussian elimination with partial pivoting
5
vector<double> gaussianSolve(vector<vector<double>> A, vector<double> b) {
6
    int n = (int)A.size();
7
    for (int i = 0; i < n; ++i) A[i].push_back(b[i]); // augment
8
    for (int col = 0; col < n; ++col) {
9
        // Pivot
10
        int piv = col;
11
        for (int r = col + 1; r < n; ++r)
12
            if (fabs(A[r][col]) > fabs(A[piv][col])) piv = r;
13
        if (fabs(A[piv][col]) < 1e-12) {
14
            // Singular; return zeros as a fallback
15
            return vector<double>(n, 0.0);
16
        }
17
        swap(A[piv], A[col]);
18
        // Normalize row
19
        double diag = A[col][col];
20
        for (int j = col; j <= n; ++j) A[col][j] /= diag;
21
        // Eliminate other rows
22
        for (int r = 0; r < n; ++r) if (r != col) {
23
            double factor = A[r][col];
24
            if (fabs(factor) < 1e-18) continue;
25
            for (int j = col; j <= n; ++j) A[r][j] -= factor * A[col][j];
26
        }
27
    }
28
    vector<double> x(n);
29
    for (int i = 0; i < n; ++i) x[i] = A[i][n];
30
    return x;
31
}
32
33
// Fit polynomial of degree d to (x, y) via normal equations, return coefficients c[0..d]
34
vector<double> fitPolynomial(const vector<double>& x, const vector<double>& y, int d) {
35
    int n = (int)x.size();
36
    int k = d + 1;
37
    // Precompute sums S_m = sum x_i^m for m = 0..2d
38
    vector<double> S(2*d + 1, 0.0);
39
    for (int i = 0; i < n; ++i) {
40
        double xp = 1.0;
41
        for (int m = 0; m <= 2*d; ++m) {
42
            S[m] += xp;
43
            xp *= x[i];
44
        }
45
    }
46
    // Build normal matrix G and right-hand side h
47
    vector<vector<double>> G(k, vector<double>(k, 0.0));
48
    vector<double> h(k, 0.0);
49
    for (int r = 0; r < k; ++r) {
50
        for (int c = 0; c < k; ++c) {
51
            G[r][c] = S[r + c]; // sum x^(r+c)
52
        }
53
        for (int i = 0; i < n; ++i) {
54
            double xp = 1.0;
55
            for (int t = 0; t < r; ++t) xp *= x[i];
56
            h[r] += xp * y[i];
57
        }
58
    }
59
    return gaussianSolve(G, h);
60
}
61
62
// Compute MDL (in nats) for polynomial of degree d using Gaussian likelihood and BIC penalty
63
// MDL = -log L(\hat{theta}) + (k/2) log n, with k = d + 1 (coefficients) + 1 (variance)
64
// Here we fold variance into the likelihood via sigma^2_hat = RSS/n; the penalty uses k_params = d + 1.
65
double mdlPolynomial(const vector<double>& x, const vector<double>& y, int d) {
66
    int n = (int)x.size();
67
    vector<double> c = fitPolynomial(x, y, d);
68
    // Compute predictions and RSS
69
    double RSS = 0.0;
70
    for (int i = 0; i < n; ++i) {
71
        double xp = 1.0, yi_hat = 0.0;
72
        for (int j = 0; j <= d; ++j) {
73
            yi_hat += c[j] * xp;
74
            xp *= x[i];
75
        }
76
        double r = y[i] - yi_hat;
77
        RSS += r * r;
78
    }
79
    double sigma2_hat = RSS / max(1, n);
80
    if (sigma2_hat <= 1e-30) sigma2_hat = 1e-30; // numerical guard
81
    // Negative log-likelihood at MLE for Gaussian
82
    double nll = 0.5 * n * (1.0 + log(2.0 * M_PI * sigma2_hat));
83
    int k_params = d + 1; // effective free parameters for mean model; variance is accounted in nll
84
    double penalty = 0.5 * k_params * log((double)max(1, n));
85
    return nll + penalty;
86
}
87
88
int main() {
89
    ios::sync_with_stdio(false);
90
    cin.tie(nullptr);
91
92
    // Generate synthetic data: y = 2 - 0.5 x + 0.3 x^2 + noise
93
    int n = 200;
94
    vector<double> x(n), y(n);
95
    mt19937 rng(42);
96
    normal_distribution<double> noise(0.0, 0.5);
97
    for (int i = 0; i < n; ++i) {
98
        x[i] = -3.0 + 6.0 * (double)i / (n - 1);
99
        double yi = 2.0 - 0.5 * x[i] + 0.3 * x[i] * x[i];
100
        y[i] = yi + noise(rng);
101
    }
102
103
    int best_d = 0;
104
    double best_mdl = numeric_limits<double>::infinity();
105
    int d_max = 10;
106
    for (int d = 0; d <= d_max; ++d) {
107
        double mdl = mdlPolynomial(x, y, d);
108
        cout << "degree=" << d << " MDL(nats)=" << mdl << "\n";
109
        if (mdl < best_mdl) { best_mdl = mdl; best_d = d; }
110
    }
111
    cout << "Best degree by MDL: " << best_d << "\n";
112
    cout << "(To convert nats to bits, divide by log(2).)\n";
113
    return 0;
114
}
115

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute MDL in nats for k-bin histogram over data in [minV, maxV]
5
// MDL = -sum_j c_j log(c_j / n) + ((k-1)/2) log n
6
// Assumes bins are equal-width and probabilities estimated by MLE c_j / n.
7
double mdlHistogram(const vector<double>& data, int k, double minV, double maxV) {
8
    int n = (int)data.size();
9
    if (k <= 0) return numeric_limits<double>::infinity();
10
    vector<int> cnt(k, 0);
11
    double width = (maxV - minV) / k;
12
    if (width <= 0) return numeric_limits<double>::infinity();
13
    for (double v : data) {
14
        int b = (int)floor((v - minV) / width);
15
        if (b < 0) b = 0;
16
        if (b >= k) b = k - 1; // include right endpoint
17
        cnt[b]++;
18
    }
19
    // Negative log-likelihood at MLE for multinomial: -sum c_j log(c_j / n)
20
    double nll = 0.0;
21
    for (int c : cnt) {
22
        if (c > 0) nll += - (double)c * (log((double)c) - log((double)n)); // -c log(c/n)
23
    }
24
    // Penalty for (k-1) free parameters
25
    double penalty = 0.5 * max(0, k - 1) * log((double)max(1, n));
26
    return nll + penalty;
27
}
28
29
int main() {
30
    ios::sync_with_stdio(false);
31
    cin.tie(nullptr);
32
33
    // Generate data: mixture of two Gaussians on [0, 1]
34
    int n = 5000;
35
    mt19937 rng(123);
36
    normal_distribution<double> g1(0.3, 0.05), g2(0.7, 0.05);
37
    uniform_real_distribution<double> U(0.0, 1.0);
38
    vector<double> data; data.reserve(n);
39
    for (int i = 0; i < n; ++i) {
40
        double v = (U(rng) < 0.5) ? g1(rng) : g2(rng);
41
        v = min(1.0, max(0.0, v));
42
        data.push_back(v);
43
    }
44
45
    double minV = 0.0, maxV = 1.0;
46
    int k_best = 1; double best_mdl = numeric_limits<double>::infinity();
47
    int k_max = 100;
48
    for (int k = 1; k <= k_max; ++k) {
49
        double mdl = mdlHistogram(data, k, minV, maxV);
50
        if (mdl < best_mdl) { best_mdl = mdl; k_best = k; }
51
        if (k % 10 == 0) {
52
            cout << "k=" << k << " MDL(nats)=" << mdl << "\n";
53
        }
54
    }
55
    cout << "Best bin count by MDL: " << k_best << "\n";
56
    cout << "(Equal-width bins; convert to bits by dividing by log(2).)\n";
57
    return 0;
58
}
59