Rejection Sampling

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Target: p(x) = 2x on [0,1], normalized
5
// Proposal: q(x) = 1 on [0,1] (Uniform), easy to sample
6
// Envelope: M = sup_x p(x)/q(x) = 2
7
8
int main() {
9
    ios::sync_with_stdio(false);
10
    cin.tie(nullptr);
11
12
    std::mt19937 rng(12345);                 // Fixed seed for reproducibility
13
    std::uniform_real_distribution<double> U01(0.0, 1.0);
14
15
    const double M = 2.0;                    // Envelope constant
16
    const size_t N = 200000;                 // Number of accepted samples desired
17
18
    vector<double> samples; samples.reserve(N);
19
    size_t proposals = 0;
20
    while (samples.size() < N) {
21
        double x = U01(rng);                 // X ~ q (Uniform[0,1])
22
        double u = U01(rng);                 // U ~ Uniform[0,1]
23
        double px = 2.0 * x;                 // p(x)
24
        double qx = 1.0;                     // q(x)
25
        double accept_prob = px / (M * qx);  // = x, since M=2
26
        proposals++;
27
        if (u <= accept_prob) {
28
            samples.push_back(x);
29
        }
30
    }
31
32
    // Diagnostic: empirical mean should be E[X] = 2/3 for p(x)=2x
33
    double mean = accumulate(samples.begin(), samples.end(), 0.0) / samples.size();
34
    cout << fixed << setprecision(6);
35
    cout << "Accepted: " << samples.size() << "/" << proposals
36
         << " (rate = " << (double)samples.size()/proposals << ")\n";
37
    cout << "Sample mean ~ " << mean << " (theory 2/3 = 0.666667)\n";
38
39
    return 0;
40
}
41

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Target: Uniform on the unit disk {(x,y): x^2 + y^2 <= 1}
5
// Proposal: Uniform on square [-1,1] x [-1,1]
6
// Envelope: M = area(square)/area(circle) = 4 / pi; acceptance test reduces to x^2 + y^2 <= 1
7
8
int main(){
9
    ios::sync_with_stdio(false);
10
    cin.tie(nullptr);
11
12
    std::mt19937 rng(2024);
13
    std::uniform_real_distribution<double> U(-1.0, 1.0);
14
15
    const size_t N = 200000;  // desired accepted points
16
    size_t accepted = 0, proposals = 0;
17
    vector<pair<double,double>> pts; pts.reserve(N);
18
19
    while (accepted < N) {
20
        double x = U(rng);
21
        double y = U(rng);
22
        proposals++;
23
        if (x*x + y*y <= 1.0) { // accept if inside circle
24
            pts.emplace_back(x,y);
25
            accepted++;
26
        }
27
    }
28
29
    cout << fixed << setprecision(6);
30
    double rate = (double)accepted / proposals; // should be close to pi/4
31
    cout << "Accepted: " << accepted << "/" << proposals
32
         << " (rate ~ " << rate << ") vs theory pi/4 = " << (M_PI/4.0) << "\n";
33
34
    // Optional: estimate pi via 4 * acceptance_rate
35
    cout << "Estimated pi ~ " << 4.0 * rate << "\n";
36
37
    return 0;
38
}
39

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Goal: Sample from target with unnormalized density \tilde{p}(x) = exp(-x^4) on R.
5
// Proposal: q(x) = Normal(0, 1) with density phi(x) = (1/sqrt(2*pi)) * exp(-x^2/2).
6
// We need M >= sup_x \tilde{p}(x) / q(x). Analytic bound is possible but here we numerically
7
// overestimate on a wide grid and add a generous safety margin. WARNING: In production,
8
// prefer a proven bound; numerical search can under-estimate if the grid is too coarse.
9
10
struct NormalPDF {
11
    // log density of N(0,1)
12
    static double logpdf(double x) {
13
        static const double LOG_SQRT_2PI = 0.5*log(2.0*M_PI);
14
        return -LOG_SQRT_2PI - 0.5*x*x;
15
    }
16
};
17
18
// Unnormalized target log-density: log tilde{p}(x) = -x^4
19
static inline double log_tilde_p(double x) { return -(x*x)*(x*x); }
20
21
int main(){
22
    ios::sync_with_stdio(false);
23
    cin.tie(nullptr);
24
25
    std::mt19937 rng(7);
26
    std::normal_distribution<double> N01(0.0, 1.0);
27
    std::uniform_real_distribution<double> U01(0.0, 1.0);
28
29
    // 1) Find a conservative M by scanning r(x) = tilde{p}(x)/q(x) over [-R, R]
30
    // The true sup occurs near x=0 for this pair; tails are safe since exp(-x^4) decays faster than Gaussian.
31
    const double R = 6.0;          // search radius
32
    const int G = 200000;          // grid points
33
    double max_log_r = -numeric_limits<double>::infinity();
34
    for (int i = 0; i <= G; ++i) {
35
        double x = -R + (2.0*R)*i / G;
36
        double lr = log_tilde_p(x) - NormalPDF::logpdf(x); // log ratio
37
        if (lr > max_log_r) max_log_r = lr;
38
    }
39
    double M = exp(max_log_r) * 1.10; // add 10% safety margin
40
41
    cerr << fixed << setprecision(6);
42
    cerr << "Estimated M (with margin): " << M << "\n";
43
44
    // 2) Draw samples via rejection using a log-safe acceptance test
45
    const size_t N = 100000; // number of accepted samples desired
46
    vector<double> xs; xs.reserve(N);
47
    size_t proposals = 0;
48
49
    while (xs.size() < N) {
50
        double x = N01(rng);     // X ~ q
51
        double u = U01(rng);     // U ~ Uniform(0,1)
52
        // log acceptance = log tilde{p}(x) - log M - log q(x)
53
        double log_alpha = log_tilde_p(x) - log(M) - NormalPDF::logpdf(x);
54
        proposals++;
55
        if (log(u) <= log_alpha) {
56
            xs.push_back(x);
57
        }
58
    }
59
60
    // Diagnostics: symmetry implies mean ~ 0; variance is finite (~0.78 for this target)
61
    double mean = accumulate(xs.begin(), xs.end(), 0.0) / xs.size();
62
    double var = 0.0; for (double v : xs) var += (v-mean)*(v-mean); var /= xs.size();
63
64
    cout << fixed << setprecision(6);
65
    cout << "Accepted: " << xs.size() << "/" << proposals
66
         << " (rate = " << (double)xs.size()/proposals << ")\n";
67
    cout << "Mean ~ " << mean << ", Variance ~ " << var << "\n";
68
69
    return 0;
70
}
71