acceptreject.m

%% Acceptance-rejection method for distributionally random numbers
% The acceptance-rejection method is another way to generate random numbers
% according to a specific distribution. Conceptually it can be thought of
% in the following way.

% Let f_X(x) be the pdf of the distribution from which we want random
% numbers.
% We define f on a subset of the the real line.
% Let A(f) denote the area underneath our pdf on said interval.
% The reason we need to define f on a subset is because A(f) would equal 1
% (i.e. we are simply computing the CDF) if we allowed for the entire real
% line.
% Now we generate some uniform random variables U.
% So for every pair of coordinates ( x in X , U*f_X(x) ) this is randomly
% and uniformly distributed across A(f).

% Now let us introduce another pdf, g_Y(x), for the distribution Y 
% And let us introduce some constant c > 1
% Now doing the same as before WITH x from the original distribution X, we
% can plot the pairs ( x in X , U*c*f_Y(x) ) and these will be randomly and
% uniformly distributed across A(cg)

% Now comes the method:
% We simply need to run through the numbers generated by U*c*g_Y(x) and
% check:
%  if they are > f_X(x) reject x in X (i.e. set to 0 or exclude)
%  elseif they are <= f_X(x) then keep x in X

% The important thing to note is that whatever the nsample is, there WILL
% BE FEWER random deviates! 
% i.e. nsample = 1,000,000 : random deviates = 750,000 

% ------------------------------------
% A note on finding  the constant c
% ------------------------------------
% We could of course take c to be large and then we absolutely know that
% the c*g will be greater than f, but this will lead to a low acceptance
% ratio. Instead we can compute analytically a good value for c by
% taking (f/g)'=0 where ' is the derivative NOT transpose. Solving this will
% give us a value for c (although may not be bulletproof).

%% Example: Normally Dist Random Deviates from Laplace Distribution

% We will seek to find normally distributed random deviates by using the
% Laplace distribution as the majorant function.
%  Recall Laplace Distribution PDF is
%  (1/2b)*exp( -abs(x-mu)/b )
% Since we are looking for standard normally distributed deviates we can
% take:
% b = 1 ...AND... mu = 0

clear all
close all

nsample = 10^6 ; % Number of random deviates required

% Parameters for our majorant function
c = sqrt(2*exp(1)/pi) ; % Constant for multipying the majorant function

% Generate a uniformly distributed random number
U = rand([1,nsample]) ;

% Generate a Laplace distributed random number
X = -1*log(2*(1-U)).*(U>=0.5) + 1*log(2*U).*(U<0.5) ;

% Majorant Function:
% Take X (our Laplace distributed random number) and input it into the pdf
% g for our Laplace distribution
g = 0.5*exp(-1*X).*(X>=0) + 0.5*exp(X).*(X<0) ;

% Minor Function:
% Take X (our Laplace distributed random number) and input it into the pdf
% f for our Normal distribution
f = normpdf(X) ; 

% Scale our Majorant function as required
maj = c*rand([1,nsample]).*g ;

% Algorithm: test for those U*c*g <= f
for i = 1:nsample
    if maj(i) > f(i)
        N(i) = 0 ; 
    else
        N(i) = X(i) ;
    end
end

% Gather the non-zero elements
Z = nonzeros(N)' ;

%% Statistics and plots

a = -4 ; % Lower truncation limit
b = 4 ; % Upper truncation limit
deltax = 0.1 ; % Grid step
x = a:deltax:b ; % Grid

% To check the acceptance ratio
accept_ratio = length(Z)/nsample
format long

% We can output the value of c
analytical_c = c
numerical_c = max(f./g)

close all

figure(1)
histogram(Z,x,'Normalization','pdf') ;
hold on
fx = 1/sqrt(2*pi)*exp(-(x.^2)/2) ;
gx = 0.5*exp(-abs(x)) ;
plot(x,fx,x,c*gx,'g')
xlabel('x')
legend('Sampled f(x)','Theoretical f(x)','Majorant function cg(x)')
title('Standard normal distribution using the accept-reject algorithm')