Random Variables & Distributions

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct DiscreteRV {
5
    vector<double> prefix; // prefix[i] = P(X <= i)
6
    explicit DiscreteRV(vector<double> probs) {
7
        if (probs.empty()) throw runtime_error("Empty probability vector");
8
        // Ensure non-negative and finite
9
        for (double p : probs) {
10
            if (!(p >= 0.0) || !isfinite(p)) throw runtime_error("Invalid probability");
11
        }
12
        // Normalize to sum to 1 (robust against rounding)
13
        long double sum = 0.0L;
14
        for (double p : probs) sum += p;
15
        if (sum <= 0.0L) throw runtime_error("Sum of probabilities must be positive");
16
        vector<double> pn;
17
        pn.reserve(probs.size());
18
        for (double p : probs) pn.push_back(static_cast<double>(p / sum));
19
        // Build prefix CDF; clamp the last entry to 1 exactly
20
        prefix.resize(pn.size());
21
        partial_sum(pn.begin(), pn.end(), prefix.begin());
22
        prefix.back() = 1.0; // guard against tiny drift like 0.9999999997
23
    }
24
    // Return P(X = k)
25
    double pmf(size_t k) const {
26
        if (k >= prefix.size()) return 0.0;
27
        double prev = (k == 0 ? 0.0 : prefix[k-1]);
28
        return max(0.0, prefix[k] - prev);
29
    }
30
    // Return P(X <= k)
31
    double cdf(size_t k) const {
32
        if (k >= prefix.size()) return 1.0;
33
        return prefix[k];
34
    }
35
    // Sample using u ~ Uniform[0,1)
36
    size_t sample(mt19937 &gen) const {
37
        uniform_real_distribution<double> U(0.0, 1.0);
38
        double u = U(gen);
39
        // upper_bound finds first prefix[i] > u, which is the correct index
40
        auto it = upper_bound(prefix.begin(), prefix.end(), u);
41
        size_t idx = static_cast<size_t>(it - prefix.begin());
42
        if (idx >= prefix.size()) idx = prefix.size() - 1; // safety
43
        return idx;
44
    }
45
};
46
47
int main() {
48
    ios::sync_with_stdio(false);
49
    cin.tie(nullptr);
50
51
    // Loaded die probabilities for faces 1..6
52
    vector<double> probs = {0.05, 0.10, 0.20, 0.25, 0.25, 0.15};
53
    DiscreteRV die(probs);
54
55
    // RNG: seed once
56
    random_device rd;
57
    mt19937 gen(rd());
58
59
    // Empirical test
60
    const int N = 200000; // number of samples
61
    vector<int> counts(6, 0);
62
    for (int i = 0; i < N; ++i) {
63
        size_t face_idx = die.sample(gen); // 0..5
64
        ++counts[face_idx];
65
    }
66
67
    cout << fixed << setprecision(4);
68
    cout << "Face\tPMF(th)\tPMF(emp)\tCDF(th)\n";
69
    for (size_t k = 0; k < 6; ++k) {
70
        double pmf_th = die.pmf(k);
71
        double pmf_emp = static_cast<double>(counts[k]) / N;
72
        double cdf_th = die.cdf(k);
73
        cout << (k+1) << '\t' << pmf_th << '\t' << pmf_emp << '\t' << cdf_th << "\n";
74
    }
75
76
    return 0;
77
}
78

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct ExponentialRV {
5
    double lambda; // rate > 0
6
    explicit ExponentialRV(double lam) : lambda(lam) {
7
        if (!(lambda > 0.0) || !isfinite(lambda)) throw runtime_error("lambda must be positive and finite");
8
    }
9
    // PDF and CDF
10
    double pdf(double x) const { return (x < 0 ? 0.0 : lambda * exp(-lambda * x)); }
11
    double cdf(double x) const { return (x < 0 ? 0.0 : 1.0 - exp(-lambda * x)); }
12
    // Quantile (inverse CDF). Use 1-u to avoid log(0) when u=0.
13
    double inv_cdf(double u) const {
14
        if (!(u >= 0.0 && u <= 1.0)) throw runtime_error("u must be in [0,1]");
15
        if (u == 1.0) return numeric_limits<double>::infinity();
16
        return -log(1.0 - u) / lambda;
17
    }
18
    // Sample using inverse transform with U ~ Uniform(0,1)
19
    double sample(mt19937 &gen) const {
20
        uniform_real_distribution<double> U(0.0, 1.0);
21
        double u = U(gen);
22
        // ensure u in (0,1) to avoid returning -log(1) = -0
23
        if (u == 0.0) u = nextafter(0.0, 1.0);
24
        return inv_cdf(u);
25
    }
26
};
27
28
int main() {
29
    ios::sync_with_stdio(false);
30
    cin.tie(nullptr);
31
32
    double lambda = 2.0; // mean = 1/lambda = 0.5
33
    ExponentialRV expo(lambda);
34
35
    random_device rd;
36
    mt19937 gen(rd());
37
38
    const int N = 200000; // samples
39
    double sum = 0.0, sumsq = 0.0;
40
    int count_gt_1 = 0;
41
    for (int i = 0; i < N; ++i) {
42
        double x = expo.sample(gen);
43
        sum += x;
44
        sumsq += x * x;
45
        if (x > 1.0) ++count_gt_1;
46
    }
47
48
    double mean_emp = sum / N;
49
    double var_emp = sumsq / N - mean_emp * mean_emp;
50
51
    cout << fixed << setprecision(6);
52
    cout << "Theoretical mean: " << 1.0 / lambda << "\n";
53
    cout << "Empirical mean:   " << mean_emp << "\n";
54
    cout << "Theoretical P(X>1): " << exp(-lambda * 1.0) << "\n";
55
    cout << "Empirical P(X>1):   " << static_cast<double>(count_gt_1) / N << "\n";
56
    cout << "Empirical variance (for reference): " << var_emp << " (theoretical = " << 1.0/(lambda*lambda) << ")\n";
57
58
    return 0;
59
}
60

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct NormalBM {
5
    double mu, sigma;
6
    bool has_spare = false; // cache the second variate
7
    double spare = 0.0;
8
    explicit NormalBM(double m, double s) : mu(m), sigma(s) {
9
        if (!(sigma > 0.0) || !isfinite(mu) || !isfinite(sigma)) throw runtime_error("Invalid parameters");
10
    }
11
    double sample(mt19937 &gen) {
12
        if (has_spare) { has_spare = false; return mu + sigma * spare; }
13
        uniform_real_distribution<double> U(0.0, 1.0);
14
        double u1 = U(gen);
15
        double u2 = U(gen);
16
        // Avoid log(0)
17
        if (u1 == 0.0) u1 = nextafter(0.0, 1.0);
18
        double r = sqrt(-2.0 * log(u1));
19
        double theta = 2.0 * M_PI * u2;
20
        double z0 = r * cos(theta);
21
        double z1 = r * sin(theta);
22
        spare = z1; has_spare = true;
23
        return mu + sigma * z0;
24
    }
25
    static double pdf(double x, double mu, double sigma) {
26
        const double inv = 1.0 / (sigma * sqrt(2.0 * M_PI));
27
        double z = (x - mu) / sigma;
28
        return inv * exp(-0.5 * z * z);
29
    }
30
};
31
32
int main() {
33
    ios::sync_with_stdio(false);
34
    cin.tie(nullptr);
35
36
    double mu = 0.0, sigma = 1.0; // standard normal
37
    NormalBM norm(mu, sigma);
38
    random_device rd; mt19937 gen(rd());
39
40
    const int N = 200000; // samples
41
    const int B = 21;     // histogram bins
42
    double lo = mu - 4*sigma, hi = mu + 4*sigma;
43
    double width = (hi - lo) / B;
44
45
    vector<int> hist(B, 0);
46
    double sum = 0.0, sumsq = 0.0;
47
48
    for (int i = 0; i < N; ++i) {
49
        double x = norm.sample(gen);
50
        sum += x; sumsq += x * x;
51
        if (x >= lo && x < hi) {
52
            int b = static_cast<int>((x - lo) / width);
53
            if (b == B) b = B - 1; // edge case
54
            ++hist[b];
55
        }
56
    }
57
58
    double mean_emp = sum / N;
59
    double var_emp = sumsq / N - mean_emp * mean_emp;
60
61
    cout << fixed << setprecision(6);
62
    cout << "Empirical mean: " << mean_emp << ", variance: " << var_emp << " (theoretical: 0, 1)\n\n";
63
64
    cout << "bin_center\temp_pdf\tth_pdf\n";
65
    for (int b = 0; b < B; ++b) {
66
        double center = lo + (b + 0.5) * width;
67
        double emp_pdf = static_cast<double>(hist[b]) / (N * width); // normalize counts to density
68
        double th_pdf = NormalBM::pdf(center, mu, sigma);
69
        cout << center << '\t' << emp_pdf << '\t' << th_pdf << '\n';
70
    }
71
72
    return 0;
73
}
74