Rejection Sampling

Rejection sampling is a technique used to generate random values from a target probability distribution. The key idea is to sample values from a proposal distribution that is easy to sample from, and reject values that don’t satisfy the desired distribution.

The algorithm works as follows:

Sample a value x from the proposal distribution
Sample a value u from the uniform distribution U(0,1)
If u ≤ f(x)/M, accept x; otherwise reject x (f is the target PDF, M is a constant s.t. f(x)/M ≤ 1)
Repeat until desired number of samples are accepted

This ensures the accepted values follow the target distribution f(x). The acceptance rate depends on the proposal distribution; a close match gives higher acceptance.

Example in Java:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
import java.util.Random;

class RejectionSampling {

  public static double sample(double min, double max) {
    Random rand = new Random();
    double x, u;
    do {
      x = rand.nextDouble() * (max - min) + min;
      u = rand.nextDouble();
    } while (u > Math.exp(-x*x/2)/0.39894); 
    return x;
  }

  public static void main(String[] args) {
    double x = sample(-4, 4);
    System.out.println(x);
  }

}

This samples from the standard normal distribution by using the uniform distribution as the proposal.

Examples in C++ and Python follow a similar structure. The key is using a proposal distribution, sampling values and accepting based on the target distribution.

Rejection Sampling is a technique in probability and statistics, in the realm of Monte Carlo methods, for generating samples from a complex distribution. If you want to generate samples from a target distribution but find it difficult to do so directly, you can sample from an easier distribution and then apply a rejection criterion to keep or discard each sample. Essentially, you “reject” samples that don’t fit the criteria of the distribution you are interested in.

Example Code

Java

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
import java.util.Random;

public class RejectionSampling {
    public static void main(String[] args) {
        Random random = new Random();
        int accepted = 0;
        int trials = 10000;

        for (int i = 0; i < trials; ++i) {
            double x = random.nextDouble();  // Sample x uniformly between 0 and 1
            double y = random.nextDouble();  // Sample y uniformly between 0 and 1

            if (y <= Math.pow(x, 2)) {  // Rejection criterion
                accepted++;
            }
        }
        double estimatedArea = (double) accepted / trials;
        System.out.println("Estimated area under x^2 curve: " + estimatedArea);
    }
}

C++

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <iostream>
#include <random>

int main() {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<> dis(0.0, 1.0);

    int accepted = 0;
    int trials = 10000;

    for (int i = 0; i < trials; ++i) {
        double x = dis(gen);  // Sample x uniformly between 0 and 1
        double y = dis(gen);  // Sample y uniformly between 0 and 1

        if (y <= x * x) {  // Rejection criterion
            accepted++;
        }
    }
    double estimatedArea = static_cast<double>(accepted) / trials;
    std::cout << "Estimated area under x^2 curve: " << estimatedArea << std::endl;
    return 0;
}

Python

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
import random

accepted = 0
trials = 10000

for _ in range(trials):
    x = random.random()  # Sample x uniformly between 0 and 1
    y = random.random()  # Sample y uniformly between 0 and 1

    if y <= x ** 2:  # Rejection criterion
        accepted += 1

estimated_area = accepted / trials
print(f"Estimated area under x^2 curve: {estimated_area}")

The codes for all three languages follow the same logic: generate random points (x, y) and accept them based on whether they satisfy the criterion y <= x^2. The ratio of accepted points gives an estimation for the area under the curve y = x^2, in the range [0,1].

Very high frequency problems on rejection sampling:

Implement Rand10() Using Rand7()
Generate Random Point in a Circle