Bayesian Inference

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute log Beta function using lgamma for numerical stability
5
double log_beta(double a, double b) {
6
    return lgamma(a) + lgamma(b) - lgamma(a + b);
7
}
8
9
// Compute log binomial coefficient C(n,k) using lgamma
10
double log_choose(int n, int k) {
11
    if (k < 0 || k > n) return -INFINITY;
12
    return lgamma(n + 1.0) - lgamma(k + 1.0) - lgamma(n - k + 1.0);
13
}
14
15
int main(){
16
    ios::sync_with_stdio(false);
17
    cin.tie(nullptr);
18
19
    // Prior: Beta(alpha, beta)
20
    double alpha = 2.0; // prior pseudo-heads
21
    double beta  = 2.0; // prior pseudo-tails
22
23
    // Observed data: n trials, k successes (heads)
24
    int n = 20;
25
    int k = 14;
26
27
    // Posterior hyperparameters
28
    double alpha_post = alpha + k;
29
    double beta_post  = beta  + (n - k);
30
31
    // Posterior summaries
32
    double mean_post = alpha_post / (alpha_post + beta_post);
33
    double var_post  = (alpha_post * beta_post) /
34
                       ((alpha_post + beta_post) * (alpha_post + beta_post) * (alpha_post + beta_post + 1.0));
35
    // MAP exists only if alpha_post > 1 and beta_post > 1
36
    double map_post;
37
    if (alpha_post > 1.0 && beta_post > 1.0) {
38
        map_post = (alpha_post - 1.0) / (alpha_post + beta_post - 2.0);
39
    } else {
40
        map_post = mean_post; // fallback
41
    }
42
43
    // 95% normal-approx credible interval for theta (approximation)
44
    double sd_post = sqrt(var_post);
45
    double ci_low = max(0.0, mean_post - 1.96 * sd_post);
46
    double ci_high = min(1.0, mean_post + 1.96 * sd_post);
47
48
    cout << fixed << setprecision(6);
49
    cout << "Posterior Beta(" << alpha_post << ", " << beta_post << ")\n";
50
    cout << "Posterior mean: " << mean_post << "\n";
51
    cout << "Posterior MAP:  " << map_post  << "\n";
52
    cout << "Approx 95% CI:  [" << ci_low << ", " << ci_high << "]\n\n";
53
54
    // Posterior predictive for m future trials has Beta-Binomial distribution
55
    int m = 5; // number of future trials
56
    cout << "Predictive P(K=k for m=" << m << ")" << "\n";
57
    for (int kfuture = 0; kfuture <= m; ++kfuture) {
58
        // log p = log C(m,k) + log B(alpha'+k, beta'+m-k) - log B(alpha', beta')
59
        double logp = log_choose(m, kfuture)
60
                    + log_beta(alpha_post + kfuture, beta_post + (m - kfuture))
61
                    - log_beta(alpha_post, beta_post);
62
        cout << "k=" << kfuture << ": " << exp(logp) << "\n";
63
    }
64
65
    return 0;
66
}
67

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct PosteriorNormal {
5
    double mu_n;    // posterior mean
6
    double tau2_n;  // posterior variance
7
};
8
9
PosteriorNormal update_normal_normal(double mu0, double tau2_0, double sigma2, const vector<double>& x){
10
    int n = (int)x.size();
11
    double sumx = accumulate(x.begin(), x.end(), 0.0);
12
    double prec0 = 1.0 / tau2_0;   // prior precision
13
    double precx = n / sigma2;     // data precision
14
    double precn = prec0 + precx;  // posterior precision
15
    double mu_n = (prec0 * mu0 + (sumx / sigma2)) / precn;
16
    double tau2_n = 1.0 / precn;
17
    return {mu_n, tau2_n};
18
}
19
20
int main(){
21
    ios::sync_with_stdio(false);
22
    cin.tie(nullptr);
23
24
    // Prior: theta ~ N(mu0, tau0^2)
25
    double mu0 = 0.0;
26
    double tau2_0 = 4.0; // prior variance (tau0^2)
27
28
    // Known observation variance sigma^2
29
    double sigma2 = 1.0;
30
31
    // Observed data
32
    vector<double> x = {0.2, -0.1, 0.4, 0.3, -0.2, 0.1};
33
34
    // Posterior update
35
    PosteriorNormal post = update_normal_normal(mu0, tau2_0, sigma2, x);
36
37
    // 95% credible interval (Normal posterior)
38
    double sd_post = sqrt(post.tau2_n);
39
    double ci_low = post.mu_n - 1.96 * sd_post;
40
    double ci_high = post.mu_n + 1.96 * sd_post;
41
42
    // Posterior predictive for a new x_new: Normal(mu_n, sigma2 + tau2_n)
43
    double pred_mean = post.mu_n;
44
    double pred_var  = sigma2 + post.tau2_n;
45
    double pred_sd   = sqrt(pred_var);
46
47
    cout << fixed << setprecision(6);
48
    cout << "Posterior mean: " << post.mu_n << "\n";
49
    cout << "Posterior var:  " << post.tau2_n << "\n";
50
    cout << "95% CI:         [" << ci_low << ", " << ci_high << "]\n";
51
    cout << "Predictive mean:" << pred_mean << "; sd=" << pred_sd << "\n";
52
53
    return 0;
54
}
55

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Log posterior for Beta(alpha', beta') (up to constant)
5
double log_posterior_beta(double theta, double a, double b) {
6
    if (theta <= 0.0 || theta >= 1.0) return -INFINITY;
7
    return (a - 1.0) * log(theta) + (b - 1.0) * log1p(-theta); // log1p(-theta) = log(1 - theta)
8
}
9
10
// Reflect a point into (0,1) interval for random-walk proposals
11
double reflect01(double x){
12
    // Handle multiple crossings if step is large
13
    while (x < 0.0 || x > 1.0) {
14
        if (x < 0.0) x = -x;       // reflect below 0
15
        if (x > 1.0) x = 2.0 - x;  // reflect above 1
16
    }
17
    // Avoid exactly 0 or 1 which are invalid for log
18
    const double eps = 1e-12;
19
    if (x <= 0.0) x = eps;
20
    if (x >= 1.0) x = 1.0 - eps;
21
    return x;
22
}
23
24
int main(){
25
    ios::sync_with_stdio(false);
26
    cin.tie(nullptr);
27
28
    // Data: k successes out of n trials with Beta(alpha, beta) prior
29
    int n = 50, k = 32;
30
    double alpha = 2.0, beta = 2.0;           // prior
31
    double a_post = alpha + k;                 // posterior alpha'
32
    double b_post = beta  + (n - k);           // posterior beta'
33
34
    // MCMC settings
35
    int burn_in = 5000;
36
    int iters   = 50000;
37
    int thin    = 10;
38
    double step_sd = 0.05; // proposal std dev for random-walk in theta-space
39
40
    mt19937_64 rng(random_device{}());
41
    normal_distribution<double> norm(0.0, step_sd);
42
    uniform_real_distribution<double> unif(0.0, 1.0);
43
44
    // Initialize at posterior mean
45
    double theta = a_post / (a_post + b_post);
46
    vector<double> samples;
47
    samples.reserve((iters - burn_in) / max(1, thin));
48
49
    int accepted = 0;
50
    for (int t = 0; t < iters; ++t) {
51
        double prop = reflect01(theta + norm(rng));
52
        double logp_cur = log_posterior_beta(theta, a_post, b_post);
53
        double logp_prop = log_posterior_beta(prop, a_post, b_post);
54
        double loga = logp_prop - logp_cur; // symmetric proposal -> Hastings ratio cancels
55
        if (loga >= 0.0 || log(unif(rng)) < loga) {
56
            theta = prop;
57
            if (t >= burn_in) ++accepted;
58
        }
59
        if (t >= burn_in && ((t - burn_in) % thin == 0)) {
60
            samples.push_back(theta);
61
        }
62
    }
63
64
    // Compute posterior mean estimate from samples
65
    double mean_est = accumulate(samples.begin(), samples.end(), 0.0) / samples.size();
66
67
    cout << fixed << setprecision(6);
68
    cout << "Samples kept: " << samples.size() << "\n";
69
    cout << "Acceptance rate (post-burn-in): "
70
         << (double)accepted / max(1, iters - burn_in) << "\n";
71
    cout << "Posterior mean (MC): " << mean_est << "\n";
72
    cout << "Posterior mean (analytic): " << a_post / (a_post + b_post) << "\n";
73
74
    return 0;
75
}
76