%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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_AR2.txt','w');
if fid==-1;
   error('File could not be opened');
  else;
   disp('File opened successfully');
end;

R=30000; %desired number of accepted draws

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ex1: Approximating the truncated normal(0,4) over [-2,2] via a U[-2,2]
M=4;
ara = (normcdf(1)-normcdf(-1))*sqrt(8*pi)/4; %analytical acceptance rate

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=-2+4*rand; %draw from the source density
        U=rand; %threshold uniform
        fwig = exp(-(1/8)*thcan^2);%true kernel
        swig = 1/4;%source kernel
        if U <= fwig/(M*swig)
            keep = keep+1;
            keepmat1(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');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ex2: Approximating the truncated normal(0,4) over [-2,2] via a normal
% with moments corresponding to the truncated moments
ut=2*(normpdf(-1)-normpdf(1))/(normcdf(1)-normcdf(-1));
%truncated mean, will be zero as well given that the truncation is
%symmetric
Vt=4*(1+((-normpdf(-1)-normpdf(1))/(normcdf(1)-normcdf(-1)))-...
    ((normpdf(-1)-normpdf(1))/(normcdf(1)-normcdf(-1)))^2);%truncated variance

M=exp((-(1/8)+(1/(2*Vt)))*4);
ara = ((normcdf(1)-normcdf(-1))*sqrt(8*pi))/(M*(sqrt(2*Vt*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=sqrt(Vt)*randn; %draw from the source density
        U=rand; %threshold uniform
        if abs(thcan)<2;
            fwig = exp(-(1/8)*thcan^2);%true kernel
        else
            fwig = 0;
        end
        swig = exp(-(1/(2*Vt))*thcan^2);%source kernel
        if U <= fwig/(M*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');

% Let's see if the source densities lead to draws from the target
% density
x1=[-2:0.01:2];
f1=((normcdf(1)-normcdf(-1))*sqrt(8*pi))^(-1)*exp(-(1/8)*x1.^2); %true density

[f2,x2]=ksdensity(keepmat1,'kernel','epanechnikov','npoints',1000);
[f3,x3]=ksdensity(keepmat2,'kernel','epanechnikov','npoints',1000);

figure(1);
xl=-2;
xu=2;
yl=0;
yu=0.4;
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k');
title('AR for truncated normal');
legend('target','via U(-2,2)','via n(ut,Vt)');
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;

    