% This shows examples of generating random variables from a distribution


%% Example: Inversion method demonstrated with normal distribution

u=10;
sigma=10;


PDF=inline('1/sqrt(2*pi*sigma^2)*exp(-(X-u).^2/(2*sigma^2))','X','u','sigma')
CDFINV=inline('u+sigma*sqrt(2)*erfinv(2*Y-1)','Y','u','sigma')

N=10000; %number of values to select from the distribution
Y=rand(N,1);
X=CDFINV(Y,u,sigma)

dx=.1;
[counts bin]=hist(X,(u-4*sigma:dx:u+4*sigma))
pdfcurve=PDF(bin,u,sigma);

%normalize to match our normalized PDF
normcounts=counts/sum(counts)/dx;
%check normalization
sum(normcounts*dx)

figure
bar(bin,normcounts,'r')
hold on
plot(bin,pdfcurve,'k','linewidth', 2)
xlabel('x','Fontsize',16)
ylabel('Normalized Counts','Fontsize', 16)
legend('Generated Data','PDF(x)')


%% Example: Accept Reject method demonstrated with sin^2 distribution

% note: our PDF is that of the ground state of a 1D particle in a box 
% (psi* psi)

N=10000;
L=10;
PDF=inline('2/L*sin(pi*X/L).^2','X','L')

nn=1;
while nn<N+1
    
    %select a value on our interval (0:L) from a uniform distribution
    x_candidate=rand*L;
    %select a value in the range (0:max(PDF)) from a uniform distribution
    y_candidate=rand*2/L;
    %accept or reject
    if y_candidate<=PDF(x_candidate,L) %accept if true
        X(nn)=x_candidate;
        nn=nn+1;
    end
    
end

dx=.1;
[counts bin]=hist(X,(0:dx:L))
pdfcurve=PDF(bin,L);

%normalize to match our normalized PDF
normcounts=counts/sum(counts)/dx;
%check normalization
sum(normcounts*dx)

figure
bar(bin,normcounts,'r')
hold on
plot(bin,pdfcurve,'k','linewidth', 2)
xlabel('x','Fontsize',16)
ylabel('Normalized Counts','Fontsize', 16)
legend('Generated Data','PDF(x)')















