Lebesgue Integration

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Approximate \int_a^b f(x) dx by a Lebesgue-style simple function:
5
// 1) Sample N points X ~ Uniform[a,b]
6
// 2) Bin f(X) into B range-bins [y_j, y_{j+1})
7
// 3) Use bin mid as simple-function value and estimate measure of preimage by counts/N * (b-a)
8
// Returns the estimate and prints diagnostics.
9
struct LebesgueHistogramResult {
10
    double estimate;
11
    double y_min, y_max;
12
    vector<size_t> counts;
13
};
14
15
LebesgueHistogramResult lebesgue_hist_integral(
16
    const function<double(double)>& f,
17
    double a, double b,
18
    size_t samples, size_t bins,
19
    uint64_t seed = 42ULL
20
){
21
    mt19937_64 rng(seed);
22
    uniform_real_distribution<double> U(a, b);
23
24
    vector<double> vals;
25
    vals.reserve(samples);
26
    for (size_t i = 0; i < samples; ++i) {
27
        double x = U(rng);
28
        vals.push_back(f(x));
29
    }
30
31
    // Determine range from samples (robustified by small padding)
32
    auto [mn_it, mx_it] = minmax_element(vals.begin(), vals.end());
33
    double y_min = *mn_it, y_max = *mx_it;
34
    if (y_min == y_max) {
35
        // Constant function case
36
        double estimate = (b - a) * y_min;
37
        return {estimate, y_min, y_max, vector<size_t>(bins, 0)};
38
    }
39
    double pad = 1e-12 * max(1.0, max(fabs(y_min), fabs(y_max)));
40
    y_min -= pad; y_max += pad;
41
42
    vector<size_t> counts(bins, 0);
43
    auto bin_of = [&](double y){
44
        double t = (y - y_min) / (y_max - y_min);
45
        size_t k = (size_t) floor(t * bins);
46
        if (k >= bins) k = bins - 1; // clamp
47
        return k;
48
    };
49
50
    for (double y : vals) counts[bin_of(y)]++;
51
52
    // Simple-function value: bin midpoint; measure of preimage: count/N * (b-a)
53
    double estimate = 0.0;
54
    for (size_t k = 0; k < bins; ++k) {
55
        double y_left = y_min + (y_max - y_min) * (double)k / (double)bins;
56
        double y_right = y_min + (y_max - y_min) * (double)(k+1) / (double)bins;
57
        double y_mid = 0.5 * (y_left + y_right);
58
        double measure = (double)counts[k] / (double)samples * (b - a);
59
        estimate += y_mid * measure;
60
    }
61
62
    return {estimate, y_min, y_max, counts};
63
}
64
65
// Example function with narrow spikes of large height (bounded):
66
// On each interval [i*0.01, i*0.01 + 0.0001), value = 1000; otherwise 0.
67
// Total measure of spikes = 100 * 0.0001 = 0.01, integral = 1000 * 0.01 = 10.
68
double spiky(double x) {
69
    int i = (int) floor(x * 100.0); // i in {0,...,99} when x in [0,1)
70
    if (i < 0) i = 0; if (i > 99) i = 99;
71
    double left = i * 0.01;
72
    double right_spike = left + 0.0001; // very narrow
73
    return (x >= left && x < right_spike) ? 1000.0 : 0.0;
74
}
75
76
int main(){
77
    double a = 0.0, b = 1.0;
78
    size_t N = 200000;   // number of uniform samples
79
    size_t B = 50;       // number of y-bins
80
81
    auto res = lebesgue_hist_integral(spiky, a, b, N, B, 12345ULL);
82
83
    cout.setf(ios::fixed); cout << setprecision(6);
84
    cout << "Lebesgue histogram estimate of "+string("\\int_0^1 spiky(x) dx")+": " << res.estimate << "\n";
85
    cout << "(Ground truth is 10)\n";
86
    cout << "y-range estimated from samples: [" << res.y_min << ", " << res.y_max << "]\n";
87
    // Optional: show a few nonempty bins
88
    size_t shown = 0;
89
    for (size_t k = 0; k < res.counts.size() && shown < 5; ++k) {
90
        if (res.counts[k] > 0) {
91
            double y_left = res.y_min + (res.y_max - res.y_min) * (double)k / (double)B;
92
            double y_right = res.y_min + (res.y_max - res.y_min) * (double)(k+1) / (double)B;
93
            cout << "bin [" << y_left << ", " << y_right << ") count = " << res.counts[k] << "\n";
94
            shown++;
95
        }
96
    }
97
    return 0;
98
}
99

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
double spiky(double x) {
5
    int i = (int) floor(x * 100.0);
6
    if (i < 0) i = 0; if (i > 99) i = 99;
7
    double left = i * 0.01;
8
    double right_spike = left + 0.0001;
9
    return (x >= left && x < right_spike) ? 1000.0 : 0.0;
10
}
11
12
// Classical midpoint Riemann sum
13
double riemann_midpoint(const function<double(double)>& f, double a, double b, size_t M){
14
    double h = (b - a) / (double)M;
15
    double sum = 0.0;
16
    for (size_t i = 0; i < M; ++i){
17
        double x_mid = a + (i + 0.5) * h;
18
        sum += f(x_mid);
19
    }
20
    return sum * h;
21
}
22
23
// Monte Carlo estimate using expectation: \int_a^b f = (b-a) * E[f(U)]
24
struct MCResult { double estimate; double stderr; double ci95_halfwidth; };
25
26
MCResult monte_carlo_integral(const function<double(double)>& f, double a, double b, size_t N, uint64_t seed=2024ULL){
27
    mt19937_64 rng(seed);
28
    uniform_real_distribution<double> U(a, b);
29
    long double mean = 0.0L, M2 = 0.0L; // Welford for variance of f(U)
30
    for (size_t i = 0; i < N; ++i){
31
        double x = U(rng);
32
        double y = f(x);
33
        long double delta = y - mean;
34
        mean += delta / (long double)(i + 1);
35
        M2 += delta * (y - mean);
36
    }
37
    long double var = (N > 1) ? (M2 / (long double)(N - 1)) : 0.0L;
38
    long double est = (b - a) * mean;
39
    long double stderr = (b - a) * sqrt((double)var / (double)N);
40
    long double ci95 = 1.96L * stderr;
41
    return {(double)est, (double)stderr, (double)ci95};
42
}
43
44
int main(){
45
    double a = 0.0, b = 1.0;
46
    cout.setf(ios::fixed); cout << setprecision(6);
47
48
    // Ground truth integral = 10
49
    for (size_t M : {50UL, 100UL, 200UL, 1000UL}){
50
        double rmid = riemann_midpoint(spiky, a, b, M);
51
        cout << "Riemann midpoint (M=" << M << ") = " << rmid << "\n";
52
    }
53
54
    for (size_t N : {1000UL, 10000UL, 100000UL}){
55
        auto mc = monte_carlo_integral(spiky, a, b, N);
56
        cout << "Monte Carlo (N=" << N << ") = " << mc.estimate
57
             << ", stderr = " << mc.stderr
58
             << ", 95% CI half-width = " << mc.ci95_halfwidth << "\n";
59
    }
60
    return 0;
61
}
62

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// The Dirichlet function on [0,1]: f(x)=1 if x is rational, 0 if irrational.
5
// Mathematically, every interval contains rationals and irrationals, so
6
//  - lower Darboux sum over ANY partition is 0 (infimum on each subinterval is 0),
7
//  - upper Darboux sum is 1 (supremum on each subinterval is 1).
8
// Thus Riemann integral does not exist, but Lebesgue integral is 0 (since rationals have measure 0).
9
// This code demonstrates the theoretical Darboux sums for any partition size M.
10
11
pair<double,double> dirichlet_darboux_sums(size_t M){
12
    (void)M; // independent of partition
13
    double lower_sum = 0.0;
14
    double upper_sum = 1.0;
15
    return {lower_sum, upper_sum};
16
}
17
18
int main(){
19
    cout.setf(ios::fixed); cout << setprecision(6);
20
    for (size_t M : {4UL, 10UL, 100UL, 1000UL}){
21
        auto [L,U] = dirichlet_darboux_sums(M);
22
        cout << "Partition M=" << M << ": lower sum = " << L << ", upper sum = " << U << "\n";
23
    }
24
    cout << "Conclusion: Riemann integral fails (lower!=upper), but Lebesgue integral equals 0 because the set where f=1 has measure 0." << "\n";
25
    return 0;
26
}
27