%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;