Gibbs Sampling

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct RunningStats {
5
    long long n = 0;
6
    double mean_x = 0.0, mean_y = 0.0;
7
    double M2_x = 0.0, M2_y = 0.0, C_xy = 0.0; // for variance/covariance (Welford)
8
    void add(double x, double y) {
9
        n++;
10
        double dx = x - mean_x;
11
        double dy = y - mean_y;
12
        mean_x += dx / n;
13
        mean_y += dy / n;
14
        M2_x += dx * (x - mean_x);
15
        M2_y += dy * (y - mean_y);
16
        C_xy += (x - mean_x) * (y - mean_y);
17
    }
18
    double var_x() const { return n > 1 ? M2_x / (n - 1) : 0.0; }
19
    double var_y() const { return n > 1 ? M2_y / (n - 1) : 0.0; }
20
    double cov_xy() const { return n > 1 ? C_xy / (n - 1) : 0.0; }
21
    double corr() const {
22
        double vx = var_x(), vy = var_y();
23
        if (vx <= 0 || vy <= 0) return 0.0;
24
        return cov_xy() / sqrt(vx * vy);
25
    }
26
};
27
28
int main() {
29
    // Target: (X, Y) ~ N(0, 0; 1, 1, corr = rho)
30
    double rho = 0.8;                 // correlation
31
    int burn_in = 5000;               // discard initial iterations
32
    int thin = 5;                     // keep every 'thin'-th sample
33
    int kept = 20000;                 // number of kept samples after burn-in and thinning
34
35
    // Conditionals: X|Y=y ~ N(rho*y, 1 - rho^2); Y|X=x ~ N(rho*x, 1 - rho^2)
36
    double cond_var = 1.0 - rho * rho;
37
    double cond_sd = sqrt(max(1e-15, cond_var));
38
39
    // RNG setup
40
    std::random_device rd;
41
    std::mt19937_64 gen(rd());
42
    std::normal_distribution<double> std_normal(0.0, 1.0);
43
    std::uniform_real_distribution<double> uni01(0.0, 1.0);
44
45
    // Initialize chain (arbitrary start)
46
    double x = 0.0, y = 0.0;
47
48
    // Burn-in phase
49
    for (int t = 0; t < burn_in; ++t) {
50
        // Sample X given Y
51
        double mean_x = rho * y;
52
        x = mean_x + cond_sd * std_normal(gen);
53
        // Sample Y given X
54
        double mean_y = rho * x;
55
        y = mean_y + cond_sd * std_normal(gen);
56
    }
57
58
    // Sampling with thinning; collect running statistics
59
    RunningStats stats;
60
    int collected = 0;
61
    for (int s = 0; collected < kept; ++s) {
62
        // One Gibbs step
63
        double mean_x = rho * y;
64
        x = mean_x + cond_sd * std_normal(gen);
65
        double mean_y = rho * x;
66
        y = mean_y + cond_sd * std_normal(gen);
67
68
        if ((s + 1) % thin == 0) {
69
            stats.add(x, y);
70
            collected++;
71
        }
72
    }
73
74
    cout << fixed << setprecision(4);
75
    cout << "Estimated E[X] = " << stats.mean_x << " (true: 0)\n";
76
    cout << "Estimated E[Y] = " << stats.mean_y << " (true: 0)\n";
77
    cout << "Estimated Var[X] = " << stats.var_x() << " (true: 1)\n";
78
    cout << "Estimated Var[Y] = " << stats.var_y() << " (true: 1)\n";
79
    cout << "Estimated Corr[X,Y] = " << stats.corr() << " (true: " << rho << ")\n";
80
81
    return 0;
82
}
83

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// 2D Ising model with nearest-neighbor interactions on an N x N grid.
5
// Spins s(i,j) in {-1, +1}. Gibbs update uses logistic probability based on neighbors.
6
7
struct Ising {
8
    int N;               // grid size
9
    double beta;         // inverse temperature (coupling strength)
10
    double J;            // interaction strength (usually 1.0)
11
    double h;            // external field (bias)
12
    vector<int> spins;   // flattened N*N array of -1 or +1
13
14
    Ising(int N_, double beta_, double J_=1.0, double h_=0.0)
15
        : N(N_), beta(beta_), J(J_), h(h_), spins(N_*N_, 1) {}
16
17
    inline int idx(int i, int j) const { return i * N + j; }
18
19
    // Periodic boundary conditions helper
20
    inline int mod(int a, int m) const { int r = a % m; return r < 0 ? r + m : r; }
21
22
    // Sum of 4-neighbor spins around (i,j) with periodic boundaries
23
    int neighbor_sum(int i, int j) const {
24
        int up    = spins[idx(mod(i-1,N), j)];
25
        int down  = spins[idx(mod(i+1,N), j)];
26
        int left  = spins[idx(i, mod(j-1,N))];
27
        int right = spins[idx(i, mod(j+1,N))];
28
        return up + down + left + right;
29
    }
30
31
    // One Gibbs sweep over the grid (systematic scan)
32
    void sweep(std::mt19937_64 &gen) {
33
        std::uniform_real_distribution<double> uni01(0.0, 1.0);
34
        for (int i = 0; i < N; ++i) {
35
            for (int j = 0; j < N; ++j) {
36
                int s_ij = spins[idx(i,j)];
37
                int nb = neighbor_sum(i, j);
38
                // Local field for spin +1 vs -1. Probability of +1 is logistic(z), z = 2*beta*(h + J*nb)
39
                double z = 2.0 * beta * (h + J * static_cast<double>(nb));
40
                // Stable logistic: p = 1/(1+exp(-z))
41
                double p;
42
                if (z >= 0) {
43
                    double ez = exp(-z);
44
                    p = 1.0 / (1.0 + ez);
45
                } else {
46
                    double ez = exp(z);
47
                    p = ez / (1.0 + ez);
48
                }
49
                double u = uni01(gen);
50
                spins[idx(i,j)] = (u < p) ? +1 : -1;
51
            }
52
        }
53
    }
54
55
    // Magnetization: average spin
56
    double magnetization() const {
57
        long long sum = 0;
58
        for (int v : spins) sum += v;
59
        return static_cast<double>(sum) / static_cast<double>(N * N);
60
    }
61
};
62
63
int main() {
64
    int N = 40;            // Grid size (40x40)
65
    double beta = 0.4;     // Inverse temperature (< ~0.44 is disordered phase for J=1, h=0)
66
    double J = 1.0;
67
    double h = 0.0;
68
69
    std::random_device rd;
70
    std::mt19937_64 gen(rd());
71
72
    Ising model(N, beta, J, h);
73
74
    // Initialize spins randomly to speed up mixing from a neutral state
75
    std::uniform_real_distribution<double> uni01(0.0, 1.0);
76
    for (int i = 0; i < N; ++i) {
77
        for (int j = 0; j < N; ++j) {
78
            model.spins[model.idx(i,j)] = (uni01(gen) < 0.5) ? -1 : +1;
79
        }
80
    }
81
82
    int burn_in_sweeps = 200;
83
    for (int b = 0; b < burn_in_sweeps; ++b) model.sweep(gen);
84
85
    // Collect magnetization samples
86
    int sweeps = 1000, thin = 5;
87
    double mean_mag = 0.0;
88
    int kept = 0;
89
    for (int t = 0; t < sweeps; ++t) {
90
        model.sweep(gen);
91
        if ((t + 1) % thin == 0) {
92
            mean_mag += model.magnetization();
93
            kept++;
94
        }
95
    }
96
    mean_mag /= max(1, kept);
97
98
    cout << fixed << setprecision(4);
99
    cout << "Average magnetization (beta=" << beta << ") = " << mean_mag << "\n";
100
    cout << "(Try increasing beta towards 0.5 to see magnetization increase due to ordering.)\n";
101
    return 0;
102
}
103