Acceptance Rejection Method for Simulation

In simulation, use the acceptance rejection method to simulate a random variable X with density f(x), for which there is no direct method to use the inverse transformation to generate it. However, suppose that there is a known algorithm to generate random variable Y with density g(x) such that f(x) / g(x) is bounded, that is f(x) / g(x) ≤ c. To find c, take the derivative of f(x) / g(x) to find the maximum value. Here's the algorithm:

1. Generate Y according to the known algorithm
2. Generate U, the uniform distribution from 0 to 1
3. If U ≤ f(Y) / c*g(Y), then set X = Y
4. Otherwise, start all over

Note that to produce each copy of X, the loop runs until condition 3 is satisfied. In the end, c become the average number of runs to produce a copy of X. In the following, a random variable with f(x) = 2x*e^(-x^2) is generated using MATLAB. The chosen g(x) = e^(-x), so f(x) / g(x) = 2x*e(x-x^2), and the maximum value c turns out to be exactly 2. 

There's no way to verify that this method does indeed generate a random variable with the given f(x). Instead, to test the accuracy of the algorithm here, two values (a and b) are inputted at the end. The program first calculates the integral of f(x) from [a, b], or the probability that X falls between the range of [a, b]. Finally, the program uses the simulated copies of X to count the proportion of the copies of X that fell in that interval of [a, b]. The two values should be close to each other.

func = @(x) 2*x.*exp(-x.^2);
n = 100000;
results = zeros(1,n);
numTimes = zeros(1,n);

for i=1:n
    trialDone = false;
    count = 0;
    while ~trialDone
        count = count + 1;
        y = -log(rand);
        u = rand;
        if u <= y*exp(y-y^2)
            results(i) = y;
            trialDone = true;
        end
    end
    numTimes(i) = count;
end

a = input('Lower bound for x: ');
b = input('Upper bound for x: ');
disp 'Evaluated integral using f(x): ';
quad(func, a, b)
disp 'Simulated integral using acceptance / rejection method: ';
length(results(results>a & results
disp 'Average trial runs per simulated copy: '
mean(numTimes)

Choose values between 0 and 2.5 for a and b, for over 99.8% of the data are in that range. Upon running the script using different values, the calculated integral closely matches the frequency proportion obtained from the simulation technique using the acceptance rejection method. While this does not verify the accuracy of the method, it at least illustrates that the resultant output resembles what the expected output should be. Finally, the resultant average trial runs per simulated copy hovers closely to 2 as expected.