Markov Chain Monte Carlo (MCMC)

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute log N(x; mu, sigma^2) up to the exact constant (includes -0.5*log(2*pi))
5
double log_normal_pdf(double x, double mu, double sigma) {
6
    static const double LOG_SQRT_2PI = 0.5 * log(2.0 * M_PI);
7
    double z = (x - mu) / sigma;
8
    return -LOG_SQRT_2PI - log(sigma) - 0.5 * z * z;
9
}
10
11
// Stable log-sum-exp for two terms
12
inline double logsumexp2(double a, double b) {
13
    double m = max(a, b);
14
    return m + log(exp(a - m) + exp(b - m));
15
}
16
17
// Target: mixture 0.3*N(-3,1^2) + 0.7*N(2,0.5^2)
18
double log_target(double x) {
19
    double a = log(0.3) + log_normal_pdf(x, -3.0, 1.0);
20
    double b = log(0.7) + log_normal_pdf(x,  2.0, 0.5);
21
    return logsumexp2(a, b); // exact log density of the mixture
22
}
23
24
int main() {
25
    ios::sync_with_stdio(false);
26
    cin.tie(nullptr);
27
28
    // RNG setup (fixed seed for reproducibility)
29
    std::mt19937_64 rng(42);
30
    std::normal_distribution<double> std_normal(0.0, 1.0);
31
    std::uniform_real_distribution<double> unif01(0.0, 1.0);
32
33
    // MH parameters
34
    const int total_iters = 50000;   // total iterations including burn-in
35
    const int burn_in = 5000;        // discard initial steps
36
    const int thin = 10;             // keep every 'thin'-th sample after burn-in
37
    const double proposal_sigma = 1.0; // random-walk step size
38
39
    // Initialize chain
40
    double x = 0.0;                  // starting point
41
    double logp = log_target(x);
42
43
    int accepted = 0;
44
    vector<double> samples; samples.reserve((total_iters - burn_in) / max(1, thin));
45
46
    for (int t = 1; t <= total_iters; ++t) {
47
        // Propose y = x + Normal(0, proposal_sigma^2)
48
        double y = x + proposal_sigma * std_normal(rng);
49
        double logp_y = log_target(y);
50
51
        // Acceptance probability (symmetric proposal, so q cancels)
52
        double log_alpha = logp_y - logp; // accept with prob min(1, exp(log_alpha))
53
        bool accept = false;
54
        if (log_alpha >= 0.0) {
55
            accept = true; // always accept upward moves
56
        } else {
57
            double u = unif01(rng);
58
            accept = (log(u) < log_alpha);
59
        }
60
61
        if (accept) {
62
            x = y;
63
            logp = logp_y;
64
            ++accepted;
65
        }
66
67
        // Record sample after burn-in with thinning
68
        if (t > burn_in && ((t - burn_in) % thin == 0)) {
69
            samples.push_back(x);
70
        }
71
    }
72
73
    // Basic summaries: sample mean and variance
74
    double mean = 0.0;
75
    for (double v : samples) mean += v;
76
    mean /= samples.size();
77
78
    double var = 0.0;
79
    for (double v : samples) var += (v - mean) * (v - mean);
80
    var /= (samples.size() - 1);
81
82
    double acc_rate = static_cast<double>(accepted) / total_iters;
83
84
    cout.setf(std::ios::fixed); cout << setprecision(6);
85
    cout << "Kept samples: " << samples.size() << "\n";
86
    cout << "Acceptance rate: " << acc_rate << "\n";
87
    cout << "Sample mean: " << mean << "\n";
88
    cout << "Sample variance: " << var << "\n";
89
90
    // Print a few samples to inspect multimodality
91
    cout << "First 10 kept samples:" << "\n";
92
    for (size_t i = 0; i < min<size_t>(10, samples.size()); ++i) {
93
        cout << samples[i] << (i + 1 == min<size_t>(10, samples.size()) ? '\n' : ' ');
94
    }
95
96
    return 0;
97
}
98

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Target: (X, Y) ~ N(0, 0) with Cov = [[1, rho], [rho, 1]]
5
// Full conditionals: X|Y=y ~ N(rho*y, 1 - rho^2), Y|X=x ~ N(rho*x, 1 - rho^2)
6
7
int main() {
8
    ios::sync_with_stdio(false);
9
    cin.tie(nullptr);
10
11
    double rho = 0.8; // correlation
12
    if (rho <= -1.0 || rho >= 1.0) {
13
        cerr << "rho must be in (-1,1)" << '\n';
14
        return 1;
15
    }
16
17
    double var_cond = 1.0 - rho * rho;
18
    double sd_cond = sqrt(var_cond);
19
20
    // RNG
21
    std::mt19937_64 rng(123);
22
    std::normal_distribution<double> std_normal(0.0, 1.0);
23
24
    // Gibbs parameters
25
    const int total_iters = 50000;
26
    const int burn_in = 5000;
27
    const int thin = 10;
28
29
    // Initialize
30
    double x = 0.0, y = 0.0;
31
32
    vector<pair<double,double>> samples; samples.reserve((total_iters - burn_in) / max(1, thin));
33
34
    for (int t = 1; t <= total_iters; ++t) {
35
        // Sample X | Y=y
36
        double mean_x = rho * y;
37
        x = mean_x + sd_cond * std_normal(rng);
38
39
        // Sample Y | X=x
40
        double mean_y = rho * x;
41
        y = mean_y + sd_cond * std_normal(rng);
42
43
        if (t > burn_in && ((t - burn_in) % thin == 0)) {
44
            samples.emplace_back(x, y);
45
        }
46
    }
47
48
    // Compute sample means and covariance
49
    double mx = 0.0, my = 0.0;
50
    for (auto &p : samples) { mx += p.first; my += p.second; }
51
    mx /= samples.size();
52
    my /= samples.size();
53
54
    double sxx = 0.0, syy = 0.0, sxy = 0.0;
55
    for (auto &p : samples) {
56
        double dx = p.first - mx;
57
        double dy = p.second - my;
58
        sxx += dx * dx;
59
        syy += dy * dy;
60
        sxy += dx * dy;
61
    }
62
    sxx /= (samples.size() - 1);
63
    syy /= (samples.size() - 1);
64
    sxy /= (samples.size() - 1);
65
66
    cout.setf(std::ios::fixed); cout << setprecision(6);
67
    cout << "Kept samples: " << samples.size() << "\n";
68
    cout << "Sample mean (x, y): (" << mx << ", " << my << ")\n";
69
    cout << "Sample var(x): " << sxx << ", var(y): " << syy << "\n";
70
    cout << "Sample cov(x,y): " << sxy << "  (target cov: " << rho << ")\n";
71
72
    // Print a few samples
73
    cout << "First 5 kept samples (x, y):\n";
74
    for (size_t i = 0; i < min<size_t>(5, samples.size()); ++i) {
75
        cout << samples[i].first << " " << samples[i].second << "\n";
76
    }
77
78
    return 0;
79
}
80