Probability Axioms & Rules

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// A simple finite probability space with outcomes of type int
5
struct ProbabilitySpace {
6
    unordered_map<int, double> pmf; // P({w}) for each outcome w
7
8
    // Check normalization and non-negativity within a tolerance
9
    bool validate(double eps = 1e-12) const {
10
        double sum = 0.0;
11
        for (auto &kv : pmf) {
12
            if (kv.second < -eps) return false; // negative mass
13
            sum += kv.second;
14
        }
15
        return fabs(sum - 1.0) <= 1e-10; // allow small floating error
16
    }
17
18
    // Compute probability of an event represented as a set of outcomes
19
    double prob(const unordered_set<int> &event) const {
20
        double p = 0.0;
21
        for (int w : event) {
22
            auto it = pmf.find(w);
23
            if (it != pmf.end()) p += it->second;
24
        }
25
        return p;
26
    }
27
};
28
29
// Set utilities
30
unordered_set<int> set_union_us(const unordered_set<int>& A, const unordered_set<int>& B) {
31
    unordered_set<int> U; U.reserve(A.size() + B.size());
32
    for (int x : A) U.insert(x);
33
    for (int x : B) U.insert(x);
34
    return U;
35
}
36
37
unordered_set<int> set_intersection_us(const unordered_set<int>& A, const unordered_set<int>& B) {
38
    const unordered_set<int> *S = &A, *L = &B;
39
    if (A.size() > B.size()) { S = &B; L = &A; }
40
    unordered_set<int> I; I.reserve(S->size());
41
    for (int x : *S) if (L->count(x)) I.insert(x);
42
    return I;
43
}
44
45
bool disjoint(const unordered_set<int>& A, const unordered_set<int>& B) {
46
    const unordered_set<int> *S = &A, *L = &B;
47
    if (A.size() > B.size()) { S = &B; L = &A; }
48
    for (int x : *S) if (L->count(x)) return false;
49
    return true;
50
}
51
52
int main() {
53
    ios::sync_with_stdio(false);
54
    cin.tie(nullptr);
55
56
    // Build a fair die probability space: outcomes 1..6 each with 1/6
57
    ProbabilitySpace ps;
58
    for (int w = 1; w <= 6; ++w) ps.pmf[w] = 1.0/6.0;
59
60
    cout << boolalpha;
61
    cout << "Valid PMF? " << ps.validate() << "\n"; // should be true
62
63
    // Define events
64
    unordered_set<int> A = {2,4,6};          // even
65
    unordered_set<int> B = {4,5,6};          // >= 4
66
    unordered_set<int> C = {1};              // singleton
67
    unordered_set<int> D = {2,3};            // disjoint with C
68
69
    // Compute basic probabilities
70
    double pA = ps.prob(A);
71
    double pB = ps.prob(B);
72
    auto AUB = set_union_us(A, B);
73
    auto AIB = set_intersection_us(A, B);
74
    double pAUB = ps.prob(AUB);
75
    double pAIB = ps.prob(AIB);
76
77
    cout << fixed << setprecision(6);
78
    cout << "P(A)=" << pA << ", P(B)=" << pB << "\n";
79
    cout << "P(A∩B)=" << pAIB << ", P(A∪B)=" << pAUB << "\n";
80
81
    // Verify inclusion–exclusion: P(A∪B) = P(A)+P(B)-P(A∩B)
82
    double lhs = pAUB;
83
    double rhs = pA + pB - pAIB;
84
    cout << "Inclusion–Exclusion holds? diff=" << fabs(lhs - rhs) << "\n";
85
86
    // Verify additivity for disjoint events C and D
87
    cout << "C and D disjoint? " << disjoint(C, D) << "\n";
88
    auto CUD = set_union_us(C, D);
89
    double pC = ps.prob(C), pD = ps.prob(D), pCUD = ps.prob(CUD);
90
    cout << "P(C∪D)=" << pCUD << ", P(C)+P(D)=" << (pC + pD) << ", diff=" << fabs(pCUD - (pC + pD)) << "\n";
91
92
    // Complement rule for A
93
    unordered_set<int> Omega = {1,2,3,4,5,6};
94
    unordered_set<int> Ac;
95
    for (int w : Omega) if (!A.count(w)) Ac.insert(w);
96
    double pAc = ps.prob(Ac);
97
    cout << "P(A^c)=" << pAc << ", 1-P(A)=" << (1.0 - pA) << ", diff=" << fabs(pAc - (1.0 - pA)) << "\n";
98
99
    return 0;
100
}
101

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct Estimator {
5
    uint64_t countA = 0, countB = 0, countAB = 0, trials = 0;
6
    void add(bool A, bool B) {
7
        trials++;
8
        if (A) countA++;
9
        if (B) countB++;
10
        if (A && B) countAB++;
11
    }
12
    tuple<double,double,double> probs() const {
13
        double N = (double)trials;
14
        return {countA/N, countB/N, countAB/N};
15
    }
16
};
17
18
int main() {
19
    ios::sync_with_stdio(false);
20
    cin.tie(nullptr);
21
22
    // RNG setup
23
    std::random_device rd;
24
    std::mt19937_64 gen(rd());
25
    std::uniform_int_distribution<int> die(1,6);
26
27
    const uint64_t N = 5'000'000ULL; // 5 million trials
28
    Estimator est;
29
30
    for (uint64_t i = 0; i < N; ++i) {
31
        int x = die(gen); // roll a fair die
32
        bool A = (x % 2 == 0); // even
33
        bool B = (x >= 4);     // at least 4
34
        est.add(A, B);
35
    }
36
37
    auto [pA, pB, pAB] = est.probs();
38
    double pAUB_via_counts = pA + pB - pAB; // inclusion–exclusion
39
40
    cout << fixed << setprecision(6);
41
    cout << "Estimated P(A)=" << pA << ", P(B)=" << pB << ", P(A∩B)=" << pAB << "\n";
42
    cout << "Estimated P(A∪B) via inclusion–exclusion = " << pAUB_via_counts << "\n";
43
44
    // Theoretical values for comparison
45
    cout << "Theoretical: P(A)=0.5, P(B)=1/2, P(A∩B)=1/3, P(A∪B)=2/3\n";
46
47
    return 0;
48
}
49

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct BayesResult {
5
    double prior;       // P(D)
6
    double sens;        // P(+ | D)
7
    double spec;        // P(- | D^c)
8
    double post_pos;    // P(D | +)
9
    double post_neg;    // P(D | -)
10
};
11
12
BayesResult compute_exact(double prior, double sens, double spec) {
13
    double pD = prior;
14
    double pDc = 1.0 - pD;
15
    double pPosGivenD = sens;
16
    double pNegGivenDc = spec;
17
    double pPosGivenDc = 1.0 - spec;
18
19
    // Law of total probability for P(+)
20
    double pPos = pPosGivenD * pD + pPosGivenDc * pDc;
21
    double pNeg = 1.0 - pPos;
22
23
    // Bayes’ theorem
24
    double pDgivenPos = (pPosGivenD * pD) / pPos;        // assume pPos>0
25
    double pDgivenNeg = ( (1.0 - pPosGivenD) * pD ) / pNeg; // P(D | -)
26
27
    return {pD, sens, spec, pDgivenPos, pDgivenNeg};
28
}
29
30
pair<double,double> simulate(uint64_t N, double prior, double sens, double spec, std::mt19937_64 &gen) {
31
    bernoulli_distribution hasDisease(prior);
32
    bernoulli_distribution pos_if_D(sens);
33
    bernoulli_distribution neg_if_notD(spec);
34
35
    uint64_t countPos = 0, countDandPos = 0;
36
    uint64_t countNeg = 0, countDandNeg = 0;
37
38
    for (uint64_t i = 0; i < N; ++i) {
39
        bool D = hasDisease(gen);
40
        bool testPos;
41
        if (D) {
42
            testPos = pos_if_D(gen);
43
        } else {
44
            // If not diseased, test is positive with (1 - specificity)
45
            bernoulli_distribution pos_if_notD(1.0 - spec);
46
            testPos = pos_if_notD(gen);
47
        }
48
        if (testPos) {
49
            countPos++;
50
            if (D) countDandPos++;
51
        } else {
52
            countNeg++;
53
            if (D) countDandNeg++;
54
        }
55
    }
56
57
    double pDgivenPos_hat = countPos ? (double)countDandPos / (double)countPos : numeric_limits<double>::quiet_NaN();
58
    double pDgivenNeg_hat = countNeg ? (double)countDandNeg / (double)countNeg : numeric_limits<double>::quiet_NaN();
59
60
    return {pDgivenPos_hat, pDgivenNeg_hat};
61
}
62
63
int main() {
64
    ios::sync_with_stdio(false);
65
    cin.tie(nullptr);
66
67
    // Example parameters: low prevalence disease
68
    double prior = 0.01;   // 1% prevalence
69
    double sens  = 0.95;   // sensitivity
70
    double spec  = 0.90;   // specificity
71
72
    auto exact = compute_exact(prior, sens, spec);
73
74
    cout << fixed << setprecision(6);
75
    cout << "Exact P(D|+)= " << exact.post_pos << ", P(D|-)= " << exact.post_neg << "\n";
76
77
    std::random_device rd; std::mt19937_64 gen(rd());
78
    uint64_t N = 10'000'000ULL; // 10 million patients
79
    auto sim = simulate(N, prior, sens, spec, gen);
80
    cout << "Simulated P(D|+)= " << sim.first << ", P(D|-)= " << sim.second << "\n";
81
82
    return 0;
83
}
84