% Script to illustrate importance sampling
% based on ex. 11.4 in KPT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% Monte Carlo Integration
%%%%%%%%%%%%%%%%%%%%%%%%%

% In this exercise we assume that we actually know the full posterior as a
% standard normal density.  Lets' first derive its mean and std via MCI,
% for comparison purposes

R=1000;

nvec=randn(R,1);

MCImean=mean(nvec);
MCIstd=sqrt(mean(nvec.^2)-MCImean^2);

fprintf(fid,'Mean and std via MCI \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[MCImean MCIstd]); 
fprintf(fid,'\n');

% Importance sampling
%%%%%%%%%%%%%%%%%%%%%%%

% Following KPT, we'll use a t(0,1,2.5), a t(0,1,100), and a t(3,1,2.5), 
% and a t(3,1,100) density as
% importance functions. A t-density with 2.5 degrees of freedom has very
% fat tails and will not fit the standard normal very well. With 100 DoF's,
% the t is virtually identical to the normal and should thus provide 
% for a very suitable importance function. The last 2 IS's are way off in
% their mean. In that case, the t with the fatter tails actually performs
% better.

% mu=0, v=2.5:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'IS with t(0,1,2.5) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

mu=0;
v=2.5;
tvec=trnd_klaus(mu,1,v,R,1);% take draws from the importance function
wvec=normpdf(tvec,0,1)./tpdf_klaus(tvec,mu,1,v);%compute importance weights
ISmean=sum(tvec.*wvec)/sum(wvec);%approximate posterior mean

int=sum(tvec.^2.*wvec)/sum(wvec);%approximate E(theta^2)
ISstd=sqrt(int-ISmean^2);

fprintf(fid,'Mean and std via IS, with t(0,1,2.5) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[ISmean ISstd]); 
fprintf(fid,'\n');

wmean=mean(wvec);
int=mean(wvec.^2);
wstd=sqrt(int-wmean^2);

fprintf(fid,'Mean and std of weights, with t(0,1,2.5) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[wmean wstd]); 
fprintf(fid,'\n');

%plots
[f1,x1]=ksdensity(nvec,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(tvec,'kernel','epanechnikov','npoints',300);

figure(1);
subplot(2,2,1);
xl=-4;
xu=6;
yl=0;
yu=0.5;
plot(x1,f1,'-k',x2,f2,'--k');
title('IS with t(0,1,2.5)');
legend('true density','importance function');
axis([xl xu yl yu]);


% mu=0, v=100:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'IS with t(0,1,100) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

mu=0;
v=100;
tvec=trnd_klaus(mu,1,v,R,1);% take draws from the importance function
wvec=normpdf(tvec,0,1)./tpdf_klaus(tvec,mu,1,v);%compute importance weights
ISmean=sum(tvec.*wvec)/sum(wvec);%approximate posterior mean

int=sum(tvec.^2.*wvec)/sum(wvec);%approximate E(theta^2)
ISstd=sqrt(int-ISmean^2);

fprintf(fid,'Mean and std via IS, with t(0,1,100) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[ISmean ISstd]); 
fprintf(fid,'\n');

wmean=mean(wvec);
int=mean(wvec.^2);
wstd=sqrt(int-wmean^2);

fprintf(fid,'Mean and std of weights, with t(0,1,100) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[wmean wstd]); 
fprintf(fid,'\n');

[f1,x1]=ksdensity(nvec,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(tvec,'kernel','epanechnikov','npoints',300);
subplot(2,2,2);
xl=-4;
xu=6;
yl=0;
yu=0.5;
plot(x1,f1,'-k',x2,f2,'--k');
title('IS with t(0,1,100)');
legend('true density','importance function');
axis([xl xu yl yu]);

% mu=3, v=2.5:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'IS with t(3,1,2.5) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

mu=3;
v=2.5;
tvec=trnd_klaus(mu,1,v,R,1);% take draws from the importance function
wvec=normpdf(tvec,0,1)./tpdf_klaus(tvec,mu,1,v);%compute importance weights
ISmean=sum(tvec.*wvec)/sum(wvec);%approximate posterior mean

int=sum(tvec.^2.*wvec)/sum(wvec);%approximate E(theta^2)
ISstd=sqrt(int-ISmean^2);

fprintf(fid,'Mean and std via IS, with t(3,1,2.5) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[ISmean ISstd]); 
fprintf(fid,'\n');

wmean=mean(wvec);
int=mean(wvec.^2);
wstd=sqrt(int-wmean^2);

fprintf(fid,'Mean and std of weights, with t(3,1,2.5) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[wmean wstd]); 
fprintf(fid,'\n');

[f1,x1]=ksdensity(nvec,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(tvec,'kernel','epanechnikov','npoints',300);
subplot(2,2,3);
xl=-4;
xu=6;
yl=0;
yu=0.5;
plot(x1,f1,'-k',x2,f2,'--k');
title('IS with t(3,1,2.5)');
legend('true density','importance function');
axis([xl xu yl yu]);

% mu=3, v=100:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'IS with t(3,1,100) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

mu=3;
v=100;
tvec=trnd_klaus(mu,1,v,R,1);% take draws from the importance function
wvec=normpdf(tvec,0,1)./tpdf_klaus(tvec,mu,1,v);%compute importance weights
ISmean=sum(tvec.*wvec)/sum(wvec);%approximate posterior mean

int=sum(tvec.^2.*wvec)/sum(wvec);%approximate E(theta^2)
ISstd=sqrt(int-ISmean^2);

fprintf(fid,'Mean and std via IS, with t(3,1,100) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[ISmean ISstd]); 
fprintf(fid,'\n');

wmean=mean(wvec);
int=mean(wvec.^2);
wstd=sqrt(int-wmean^2);

fprintf(fid,'Mean and std of weights, with t(3,1,100) \n'); 
fprintf(fid,'mean\t\tstd\n'); 
fprintf(fid,'%6.3f\t%6.3f\n',[wmean wstd]); 
fprintf(fid,'\n');

[f1,x1]=ksdensity(nvec,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(tvec,'kernel','epanechnikov','npoints',300);
subplot(2,2,4);
xl=-4;
xu=6;
yl=0;
yu=0.5;
plot(x1,f1,'-k',x2,f2,'--k');
title('IS with t(3,1,100)');
legend('true density','importance function');
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;









