Probability Theory

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
int main() {
5
    ios::sync_with_stdio(false);
6
    cin.tie(nullptr);
7
8
    // Seed a high-quality Mersenne Twister engine
9
    random_device rd;
10
    mt19937 gen(rd());
11
12
    // Standard normal distribution N(0,1)
13
    normal_distribution<double> normal01(0.0, 1.0);
14
15
    const size_t N = 200000; // number of samples
16
17
    // Accumulators for mean and second moment (for variance)
18
    long double sum = 0.0L;
19
    long double sumsq = 0.0L;
20
21
    // Estimate probability P(|X| <= 1)
22
    size_t count_in = 0;
23
24
    for (size_t i = 0; i < N; ++i) {
25
        double x = normal01(gen);
26
        sum += x;
27
        sumsq += x * x;
28
        if (fabs(x) <= 1.0) ++count_in;
29
    }
30
31
    long double mean = sum / N;
32
    long double var = sumsq / N - mean * mean; // Var(X) = E[X^2] - (E[X])^2
33
    long double p_interval = static_cast<long double>(count_in) / N;
34
35
    // Theoretical values for comparison
36
    // For N(0,1): E[X]=0, Var(X)=1, and P(|X| <= 1) = erf(1/sqrt(2))
37
    long double theo_mean = 0.0L;
38
    long double theo_var = 1.0L;
39
    long double theo_p_interval = erf(1.0L / sqrtl(2.0L));
40
41
    cout.setf(std::ios::fixed); cout << setprecision(6);
42
    cout << "Samples: " << N << "\n";
43
    cout << "Estimated mean       : " << mean << " (theory: " << theo_mean << ")\n";
44
    cout << "Estimated variance    : " << var  << " (theory: " << theo_var  << ")\n";
45
    cout << "Estimated P(|X|<=1)  : " << p_interval << " (theory: " << theo_p_interval << ")\n";
46
47
    return 0;
48
}
49

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute posterior P(Disease | TestPositive) given prior, sensitivity, and specificity
5
// prior = P(D), sens = P(T+ | D), spec = P(T- | ~D)
6
double bayes_posterior(double prior, double sens, double spec) {
7
    double pD = prior;
8
    double pNotD = 1.0 - pD;
9
    double pTpos_given_D = sens;
10
    double pTpos_given_notD = 1.0 - spec; // false positive rate
11
    double pTpos = pTpos_given_D * pD + pTpos_given_notD * pNotD; // law of total probability
12
    return (pTpos_given_D * pD) / pTpos;
13
}
14
15
int main() {
16
    ios::sync_with_stdio(false);
17
    cin.tie(nullptr);
18
19
    double prior = 0.01;     // 1% prevalence
20
    double sens  = 0.95;     // sensitivity: 95%
21
    double spec  = 0.90;     // specificity: 90%
22
23
    double exact = bayes_posterior(prior, sens, spec);
24
25
    // Monte Carlo validation
26
    random_device rd; mt19937 gen(rd());
27
    bernoulli_distribution disease(prior);
28
    bernoulli_distribution tp(sens);          // true positive when disease present
29
    bernoulli_distribution tn(spec);          // true negative when disease absent
30
31
    const size_t N = 500000; // trials
32
    size_t test_pos = 0;
33
    size_t dise_and_pos = 0;
34
35
    for (size_t i = 0; i < N; ++i) {
36
        bool D = disease(gen);
37
        bool Tpos;
38
        if (D) {
39
            Tpos = tp(gen);
40
        } else {
41
            // if no disease, test is positive with prob (1 - specificity)
42
            bernoulli_distribution fp(1.0 - spec);
43
            Tpos = fp(gen);
44
        }
45
        if (Tpos) {
46
            ++test_pos;
47
            if (D) ++dise_and_pos;
48
        }
49
    }
50
51
    double mc = (test_pos > 0) ? (double)dise_and_pos / test_pos : 0.0;
52
53
    cout.setf(std::ios::fixed); cout << setprecision(6);
54
    cout << "Exact P(D|T+): " << exact << "\n";
55
    cout << "MC   P(D|T+): " << mc    << " (N=" << N << ")\n";
56
57
    return 0;
58
}
59

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
int main() {
5
    ios::sync_with_stdio(false);
6
    cin.tie(nullptr);
7
8
    random_device rd; mt19937 gen(rd());
9
    exponential_distribution<double> exp1(1.0); // mean=1, var=1
10
11
    const int n = 50;       // samples per trial (increase to see stronger CLT)
12
    const int trials = 10000; // number of trials
13
14
    // We standardize: Z = sqrt(n) * (Xbar - mu) / sigma, with mu=1, sigma=1 for Exp(1)
15
    vector<double> Z; Z.reserve(trials);
16
17
    for (int t = 0; t < trials; ++t) {
18
        double s = 0.0;
19
        for (int i = 0; i < n; ++i) s += exp1(gen);
20
        double xbar = s / n;
21
        double z = sqrt((double)n) * (xbar - 1.0) / 1.0;
22
        Z.push_back(z);
23
    }
24
25
    // Build histogram over [-4,4] with B bins
26
    const int B = 20;
27
    const double L = -4.0, R = 4.0;
28
    vector<int> hist(B, 0);
29
    for (double z : Z) {
30
        if (z < L || z >= R) continue; // ignore out of range for display
31
        int idx = (int)((z - L) / (R - L) * B);
32
        if (0 <= idx && idx < B) hist[idx]++;
33
    }
34
35
    // Print ASCII histogram
36
    int maxc = *max_element(hist.begin(), hist.end());
37
    cout << "Histogram of standardized means (approx N(0,1))\n";
38
    for (int b = 0; b < B; ++b) {
39
        double a = L + (R - L) * b / B;
40
        double c = L + (R - L) * (b + 1) / B;
41
        int bar = (maxc > 0) ? (50 * hist[b] / maxc) : 0;
42
        cout << fixed << setprecision(2) << "[" << setw(5) << a << ", " << setw(5) << c << ") ";
43
        cout << string(bar, '#') << "\n";
44
    }
45
46
    return 0;
47
}
48