%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% mcint1d_imp.m
%
% demonstration: Monte Carlo integration in 1D with importance sampling
%
% Wei Cai
% for LLNL CCMS Summer Institute Lectures Aug.1-3, 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%integrate g(x)=sqrt(1-x^2) over [0,1] to estiamte pi/4
%using random variables distributed over [0,1] with density distribution
%  f(x) = (1 - beta x^2) / (1-beta/3)
%plot the error as a function of N (number of random variables)

MAXN = 1e6;
%N = [10,20,50,100,200,500,1000];
N = [1:MAXN]';
nplot = [1:100,100:10:1000,1000:100:MAXN]';

%beta = 0.74; sig1=sqrt(0.0029);
%beta = 0.5;  sig1=sqrt(0.011);
%beta = 1e-6;  sig1=sqrt(0.05);
beta = 0.95;  sig1=sqrt(0.0221);

%generate random variable x according to 
% f(x) = (1 - beta x^2) / (1-beta/3)
% for a given beta
%
%Notice that f(x) = a g(x) + b h(x)
% where g(x) = (1 - beta x) / (1 - beta/2), a = (1-beta/2)/(1-beta/3)
%       h(x) = 6 x (1 - x)                , b = (beta/6)  /(1-beta/3)
%
%So that we can first generate random variables y and z
% according to g(y) and h(z), respectively
% x is then selected from y with probability a and from z with probabilty b

%generate y with g(y) = (1 - beta x) / (1 - beta/2)
xi=rand(MAXN,1);
y=(1-sqrt(1-xi*beta*(2-beta)))/beta;

%generate z with h(z) = 6z(1-z)
% z is the middle value of three random variable uniformly distributed in [0,1]
xi=rand(MAXN,3);
xi=sort(xi,2);
z=xi(:,2);

% selection:   x = y with probability a
%              x = z with probability 1-a
r = [y,z];               
a = (1-beta/2)/(1-beta/3);
xi=rand(MAXN,1);
ind=(xi<(1-beta/2)/(1-beta/3));
x=y.*ind+z.*(1-ind);

%To see the histogram of x, y and z, type
% hist(x,50), hist(y,50), hist(z,50)

g=sqrt(1-x.^2)*(1-beta/3)./(1-beta*x.^2);
G=cumsum(g)./[1:MAXN]';
G0=pi/4;
err=abs(G-G0);
errest=sig1./sqrt(N);

figure(1);
loglog(nplot,err(nplot),'b.',nplot,errest(nplot),'r-');
fs=18;
set(gca,'FontSize',fs);
xlabel('N');
ylabel('error');

%plot histogram of random variable x
figure(2);
hist(x,50);
fs=18;
set(gca,'FontSize',fs);
xlabel('x');
ylabel('counts');