%
% use acceptance-rejection method to generate random variables from a beta distribution with 
% alpha = 2 and beta = 1, and plot.
% 
% Here: f(x) = 2 x      0 < x < 1 
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clc; 
clear all; 
close all; 

c = 2; 
n = 50000; 
x    = zeros(1,n); 
xy   = zeros(1,n); 
rej  = zeros(1,n); 
rejy = zeros(1,n); 
irv  = 1; 
irej = 1; 
while( irv <= n )
  y = rand(1); 
  u = rand(1); 
  if( u <= (2*y)/(c*1) )   % <- 2*y is the "difficult" density to sample from ... 1 is the uniform density 
    x(irv)  = y; 
    xy(irv) = u*c;
    irv = irv+1; 
  else
    rej(irej)  = y;        % <- if more than n rejections happen Matlab will grow this vector as needed.
    rejy(irej) = u*c; 
    irej = irej+1; 
  end
end

[N,h] = hist( x, 100 ); 
% rescale N to make a PDF: 
N = N/(h(2)-h(1))/sum(N); 
figure; bar(h,N,1,'w'); hold on; 
p_truth = betapdf(h,2,1); 
plot( h, p_truth, '-go' ); 
xlabel( 'x' ); ylabel( 'PDF beta distribution' ); 
title( sprintf('number of samples = %10d',n) ); 
saveas( gcf, 'prob_4_9_plts', 'epsc' );