Measure Theory

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// A simple finite probability space on {0,1,...,n-1}
5
// The sigma-algebra here is the full power set (all subsets),
6
// so every function on the base set is measurable.
7
8
struct FiniteProbabilitySpace {
9
    vector<double> p; // probability mass function over {0..n-1}
10
11
    explicit FiniteProbabilitySpace(vector<double> probs) : p(move(probs)) {
12
        // Normalize to sum to 1 and ensure positivity
13
        double s = accumulate(p.begin(), p.end(), 0.0);
14
        if (s <= 0) throw runtime_error("Probabilities must sum to positive value");
15
        for (double &x : p) x = max(x, 0.0) / s;
16
    }
17
18
    int n() const { return (int)p.size(); }
19
20
    // Measure of a subset given by indices
21
    double measureOf(const vector<int>& subset) const {
22
        double m = 0.0;
23
        for (int i : subset) {
24
            if (i < 0 || i >= n()) throw runtime_error("Index out of range");
25
            m += p[i];
26
        }
27
        return m;
28
    }
29
30
    // Complement of a subset within {0..n-1}
31
    vector<int> complement(const vector<int>& subset) const {
32
        vector<char> in(n(), 0);
33
        for (int i : subset) {
34
            if (i < 0 || i >= n()) throw runtime_error("Index out of range");
35
            in[i] = 1;
36
        }
37
        vector<int> comp;
38
        for (int i = 0; i < n(); ++i) if (!in[i]) comp.push_back(i);
39
        return comp;
40
    }
41
42
    // Union and intersection helpers
43
    static vector<int> setUnion(const vector<int>& A, const vector<int>& B) {
44
        vector<int> U = A; U.insert(U.end(), B.begin(), B.end());
45
        sort(U.begin(), U.end()); U.erase(unique(U.begin(), U.end()), U.end());
46
        return U;
47
    }
48
    static vector<int> setIntersection(const vector<int>& A, const vector<int>& B) {
49
        vector<int> I; I.reserve(min(A.size(), B.size()));
50
        vector<int> As = A, Bs = B; sort(As.begin(), As.end()); sort(Bs.begin(), Bs.end());
51
        set_intersection(As.begin(), As.end(), Bs.begin(), Bs.end(), back_inserter(I));
52
        return I;
53
    }
54
55
    // Expectation of a function f: {0..n-1} -> R
56
    double expectation(const vector<double>& f) const {
57
        if ((int)f.size() != n()) throw runtime_error("Function size mismatch");
58
        double e = 0.0;
59
        for (int i = 0; i < n(); ++i) e += f[i] * p[i];
60
        return e;
61
    }
62
};
63
64
int main() {
65
    ios::sync_with_stdio(false);
66
    cin.tie(nullptr);
67
68
    // Define a 6-point probability space
69
    vector<double> probs = {0.05, 0.10, 0.20, 0.25, 0.15, 0.25};
70
    FiniteProbabilitySpace space(probs);
71
72
    // Define subsets A and B
73
    vector<int> A = {0, 2, 4};
74
    vector<int> B = {1, 2, 3};
75
76
    // Compute measures and check additivity for disjoint parts
77
    vector<int> Ac = space.complement(A);
78
    vector<int> AuB = FiniteProbabilitySpace::setUnion(A, B);
79
    vector<int> AnB = FiniteProbabilitySpace::setIntersection(A, B);
80
81
    double muA = space.measureOf(A);
82
    double muB = space.measureOf(B);
83
    double muAc = space.measureOf(Ac);
84
    double muAuB = space.measureOf(AuB);
85
    double muAnB = space.measureOf(AnB);
86
87
    cout << fixed << setprecision(6);
88
    cout << "mu(A) = " << muA << "\n";
89
    cout << "mu(B) = " << muB << "\n";
90
    cout << "mu(A^c) = " << muAc << " (should be 1 - mu(A) = " << 1 - muA << ")\n";
91
    cout << "mu(A ∪ B) = " << muAuB << "\n";
92
    cout << "mu(A ∩ B) = " << muAnB << "\n";
93
    cout << "mu(A) + mu(B) - mu(A ∩ B) = " << (muA + muB - muAnB) << " (inclusion-exclusion)\n\n";
94
95
    // Expectation of a function f(i) = i^2
96
    vector<double> f(space.n());
97
    for (int i = 0; i < space.n(); ++i) f[i] = 1.0 * i * i;
98
    cout << "E[f] with f(i)=i^2 = " << space.expectation(f) << "\n";
99
100
    return 0;
101
}
102

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// We approximate ∫_0^1 f(x) dx via Lebesgue-style simple functions.
5
// Approach: sample M points x ~ Uniform[0,1], compute v=f(x),
6
// bin by value into K bins [k/K, (k+1)/K), and form the lower simple sum:
7
// S_K = sum_{k=0}^{K-1} (k/K) * (measure of {x: f(x) in bin k}).
8
// For nonnegative f with max ≤ 1, S_K ≤ ∫ f and increases with K (Monotone Convergence).
9
10
static inline double f(double x) {
11
    // Choose a nonnegative, bounded-on-[0,1] function with max ≤ 1
12
    // Example: f(x) = e^{-x} ∈ [e^{-1}, 1], then we scale to [0,1]
13
    return exp(-x); // in [e^{-1}, 1]
14
}
15
16
int main() {
17
    ios::sync_with_stdio(false);
18
    cin.tie(nullptr);
19
20
    int M = 100000;     // number of samples
21
    vector<int> Ks = {8, 16, 32, 64}; // number of value bins
22
23
    mt19937_64 rng(12345);
24
    uniform_real_distribution<double> U(0.0, 1.0);
25
26
    // Pre-sample once to reuse across K (common random numbers)
27
    vector<double> vals; vals.reserve(M);
28
    for (int i = 0; i < M; ++i) {
29
        double x = U(rng);
30
        vals.push_back(f(x));
31
    }
32
33
    // Scale values to [0,1] using known max=1 and min=e^{-1} if desired.
34
    // For lower simple sums, it's safe to just cap at [0,1].
35
    // Here max is 1, so we just use v ∈ (0,1].
36
37
    // Monte Carlo benchmark: standard estimator of ∫ f
38
    double mc = accumulate(vals.begin(), vals.end(), 0.0) / M; // unbiased estimator
39
40
    cout << fixed << setprecision(6);
41
    cout << "Monte Carlo estimate of ∫_0^1 e^{-x} dx = " << mc
42
         << " (true value = 1 - e^{-1} ≈ " << (1.0 - 1.0/exp(1.0)) << ")\n\n";
43
44
    for (int K : Ks) {
45
        vector<long long> count(K, 0);
46
        for (double v : vals) {
47
            // Place v in [0,1]; v∈(e^{-1},1], but we clamp to [0,1]
48
            double w = min(1.0, max(0.0, v));
49
            int k = (w == 1.0) ? (K - 1) : int(w * K); // bin index in [0,K-1]
50
            count[k]++;
51
        }
52
        // Lower simple integral: sum over bins (lower edge)*measure(preimage)
53
        double S = 0.0;
54
        for (int k = 0; k < K; ++k) {
55
            double lower = (double)k / K;            // lower value level
56
            double measure_est = (double)count[k] / M; // fraction of domain mapping into bin
57
            S += lower * measure_est;                // domain length is 1
58
        }
59
        cout << "K = " << setw(2) << K << ": lower simple sum ≈ " << S
60
             << " (≤ true integral)." << '\n';
61
    }
62
63
    cout << "\nNote: As K increases, the lower simple sums approach the true integral from below (Monotone Convergence).\n";
64
    return 0;
65
}
66

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Demonstrate: For ν ≪ μ on a finite set, with weights p (μ) and q (ν),
5
// the Radon–Nikodym derivative w = q/p satisfies E_ν[f] = E_μ[f * w].
6
7
int main() {
8
    ios::sync_with_stdio(false);
9
    cin.tie(nullptr);
10
11
    int n = 10;
12
    mt19937_64 rng(2024);
13
    gamma_distribution<double> G(2.0, 1.0); // positive weights
14
15
    vector<double> p(n), q(n), f(n);
16
    for (int i = 0; i < n; ++i) {
17
        p[i] = G(rng) + 1e-6; // ensure positivity
18
        q[i] = G(rng) + 1e-6;
19
        f[i] = sin(0.5 * i) + 0.2 * i; // any real-valued function
20
    }
21
22
    auto normalize = [](vector<double>& v){
23
        double s = accumulate(v.begin(), v.end(), 0.0);
24
        for (double &x : v) x /= s;
25
    };
26
    normalize(p); normalize(q);
27
28
    // Compute E_ν[f]
29
    double E_nu = 0.0;
30
    for (int i = 0; i < n; ++i) E_nu += f[i] * q[i];
31
32
    // Compute importance weights w = q / p and E_μ[f * w]
33
    vector<double> w(n);
34
    for (int i = 0; i < n; ++i) w[i] = q[i] / p[i];
35
36
    double E_mu_weighted = 0.0;
37
    for (int i = 0; i < n; ++i) E_mu_weighted += f[i] * w[i] * p[i];
38
39
    cout << fixed << setprecision(10);
40
    cout << "E_nu[f]           = " << E_nu << "\n";
41
    cout << "E_mu[f * (q/p)]   = " << E_mu_weighted << "\n";
42
    cout << "Absolute error    = " << fabs(E_nu - E_mu_weighted) << "\n";
43
44
    // If we sample from μ and use weights w, we get an unbiased estimator of E_ν[f].
45
    // Demonstrate with Monte Carlo sampling from μ.
46
    int N = 100000; // samples
47
    discrete_distribution<int> cat_mu(p.begin(), p.end());
48
    double est = 0.0;
49
    for (int t = 0; t < N; ++t) {
50
        int i = cat_mu(rng);
51
        est += f[i] * w[i]; // estimator term; expectation equals E_nu[f]
52
    }
53
    est /= N;
54
    cout << "Monte Carlo (IS)  = " << est << " (with N=" << N << ")\n";
55
56
    return 0;
57
}
58