%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example for the AR Sampler, following KPT, Ex. 11.22 and 11.24
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

tic;  % start stop watch

[fid]=fopen('c:\klaus\aaec6564\mlab\logs\mod14_AR.txt','w');
if fid==-1;
   error('File could not be opened');
  else;
   disp('File opened successfully');
end;

R=5000; %desired number of accepted draws

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ex1: Approximating the triangular over [-1,1] via a U[-1,1]

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'AR for triang(-1,1) via U(-1,1) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

M1=2; %"Mwiggle - see lecture notes"
a=-1;
b=1;

keep=0; %running counter for kept draws
keepmat1=zeros(1,R);
i=0; %"i" will count total number of draws, accepted or not
while keep<R 
    thcan = a+(b-a)*rand; %draw candidate from U[a,b]
    U=rand; %draw threshold
    fwig=(1-abs(thcan)); %target kernel
    swig=1/2; %source kernel
    if U < fwig/(M1*swig);
         keep = keep+1;
         keepmat1(keep) = thcan;
    end;
   i=i+1;
        
end;

ara=1/2; %analytical acceptance rate
are=R/i;%empirical acceptance rate
fprintf(fid,'Analytical and empirical acceptance rate \n'); 
fprintf(fid,'ara\t\are\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[ara are]); 
fprintf(fid,'\n');

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'AR for triang(-1,1) via n(0,1) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

M2=1;
ara = 1/(sqrt(2*pi)); %analytical acceptance rate

keep=0; %running counter for kept draws
keepmat2=zeros(1,R);
i=0; %"i" will count total number of draws, accepted or not
    while keep < R;
        thcan=randn; %draw from the source density
        U=rand; %threshold uniform
        fwig = max([(1-abs(thcan)) 0]);%true kernel- 
        % here it's important to add a rejection rule for all draws outside
        % of the target support of [-1,1], since the source kernel is
        % unbounded;         
        swig = exp(-(1/2)*thcan^2);%source kernel
        if U <= fwig/(M2*swig)
            keep = keep+1;
            keepmat2(keep) = thcan;
        end;
        i=i+1;    
    end;
are=R/i;
fprintf(fid,'Analytical and empirical acceptance rate \n'); 
fprintf(fid,'ara\t\tare\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[ara are]); 
fprintf(fid,'\n');

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'AR for triang(-1,1) via n(0,1/6) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

sig2=1/6;
thmax = (1 + sqrt(1/3))/2; %see lecture notes
fwig_ord = 1 - abs(thmax);
swig_ord  = exp(-(1/(2*sig2))*thmax^2);
M3 = fwig_ord/swig_ord;
ara = 1/(M3*sqrt(2*pi*sig2)); %analytical acceptance rate

keep=0; %running counter for kept draws
keepmat3=zeros(1,R);
i=0; %"i" will count total number of draws, accepted or not
    while keep < R;
        thcan=sqrt(sig2)*randn; %draw from the source density
        U=rand; %threshold uniform
        fwig = max([(1-abs(thcan)) 0]);%true kernel
        swig = exp(-(1/(2*sig2))*thcan^2);%source kernel
        if U <= fwig/(M3*swig)
            keep = keep+1;
            keepmat3(keep) = thcan;
        end;
        i=i+1;    
    end;
are=R/i;
fprintf(fid,'Analytical and empirical acceptance rate \n'); 
fprintf(fid,'ara\t\tare\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[ara are]); 
fprintf(fid,'\n');

% Let's see if all three source densities lead to draws from the target
% density
[f1,x1]=ksdensity(keepmat1,'kernel','epanechnikov','npoints',1000);
[f2,x2]=ksdensity(keepmat2,'kernel','epanechnikov','npoints',1000);
[f3,x3]=ksdensity(keepmat3,'kernel','epanechnikov','npoints',1000);


figure(1);
xl=-2;
xu=2;
yl=0;
yu=1.2;
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k');
title('AR for triangular (-1,1))');
legend('via U(-1,1)','via n(0,1)','via n(0,1/6)');
axis([xl xu yl yu]);

%Next, let's plot the true target density, plus all approximation densities for
%comparison (approximation density = source kernel * Mwiggle)
x=[-1:0.01:1]; %just some support points for the x-axis
ytr=1-abs(x); %target density - triangular
tun=M1*(1/2)*(x./x); %U(-1,1) approximation density - the last term seems unnecessary, 
% but the plot function below needs a continuous function based on the
% original "x" to behave correctly.
yn1=M2*exp((-1/2)*x.^2);%normal(0,1) approximation density
yn2=M3*exp(-3*x.^2);% normal(0,1/6) approximation density

figure(2);
xl=-1;
xu=1;
yl=0;
yu=1.5;
plot(x,ytr,'-k',x,tun,'--k',x,yn1,':k',x,yn2,'-.k');
title('AR for triangular (-1,1): target and approximations');
legend('target','via U(-1,1)','via n(0,1)','via n(0,1/6)');
axis([xl xu yl yu]);



finish = toc; 
fprintf(fid,'Time elapsed in seconds \n\n');
fprintf(fid,'%6.3f\n',finish);

st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;

    