%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example in Hoff 1.2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This simply works through the example in Hoff (2009), section 1.2.1
% We are interested in a single parameter: theta, the fraction of
% individuals in a city population with with a specific infectious disease.
%
% We have a small sample of 20 individuals to provide empirical evidence
% for our estimation of theta.


%set random number seed. 
%This will allow you to EXACTLY replicate your results in re-runs.
%Use "37" to get my results exactly.

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

tic;  % start stop watch - always fun to know your run time.

% Open log file - for your text output
[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod1_1a.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sampling distribution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%If we knew theta, than a reasonable sampling model (and likelihood
%function) for the total number of infected citizens would be the binomial
%(20,theta) density.

%Let's plot the bin(20,theta) for a few values of theta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=20; %our assumed sample size

th1=0.05;
th2=0.1;
th3=0.2;
th4=0.5;

R=10000; %number of draws - your choice

y1=binornd(20,th1,R,1); %this gives us an R by 1 vector of proportions.
y2=binornd(20,th2,R,1); 
y3=binornd(20,th3,R,1); 
y4=binornd(20,th4,R,1); 

%Plot the results in a single figure, but with separate windows
figure(1);
subplot(2,2,1); %call first window of a 2 x 2 grid
seq=0:1:20; %determine x-values over which to plot
[counts,bincenters]=hist(y1,seq); 
%pcaptures # of obs. for each x-value in "counts", and captures the integers
%over which we wish to plot in "bincenters" (so bincenters = seq);
bar(bincenters,counts/R); %generates actual plot
title('$\theta=0.05$','interpreter','latex','fontsize',15);
xlabel('# of successes (infected persons)');
ylabel('probability');
set(gca,'Xlim',[-1 20]);
set(gca,'Ylim',[0 0.5]);
set(gca,'XTick',0:2:20);

subplot(2,2,2); %fill second window
seq=0:1:20; 
[counts,bincenters]=hist(y2,seq); 
bar(bincenters,counts/R); 
title('$\theta=0.1$','interpreter','latex','fontsize',15);
xlabel('# of successes (infected persons)');
ylabel('probability');
set(gca,'Xlim',[-1 20]);
set(gca,'Ylim',[0 0.5]);
set(gca,'XTick',0:2:20);

subplot(2,2,3); %fill third window
seq=0:1:20; 
[counts,bincenters]=hist(y3,seq); 
bar(bincenters,counts/R); 
title('$\theta=0.2$','interpreter','latex','fontsize',15);
xlabel('# of successes (infected persons)');
ylabel('probability');
set(gca,'Xlim',[-1 20]);
set(gca,'Ylim',[0 0.5]);
set(gca,'XTick',0:2:20);

subplot(2,2,4); %fill fourth window
seq=0:1:20;
[counts,bincenters]=hist(y4,seq); 
bar(bincenters,counts/R); 
title('$\theta=0.5$','interpreter','latex','fontsize',15);
xlabel('# of successes (infected persons)');
ylabel('probability');
set(gca,'Xlim',[-1 20]);
set(gca,'Ylim',[0 0.5]);
set(gca,'XTick',0:2:20);

%Save figure in Matlab format
to = 'c:\klaus\AAEC6564\mlab\figures\BinomExample'; %designate destination
saveas(gcf,to);%"gcf" means "get current figure"

% Export figure so it can be imported into latex
to = 'c:\klaus\AAEC6564\module1\BinomExample'; %designate destination
print('-depsc',to);%this automatically selects the current figure
% depsc stands for "encapsulated color PostScript" - something latex can
% handle well

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prior
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Infection rates in other cities range from 0.05 to 0.2, with an average
%prevalence of 0.1. Thus, a reasonable prior distribution must satisfy the
%following conditions:

% bound between 0 and 1
% continuous between bounds
% mean near 0.1

% Even with these conditions, there are many priors to choose from. In such
% cases, we add a 4th condition: computational convenience. This leads to a
% beta prior. The beta(a,b) has mean a/(a+b) and mode (a-1)/(a-1+b-1).

a0=2;
b0=20;
pr=betarnd(a0,b0,R,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Posterior
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Assume we observe y=0 (i.e. nobody of the 20 individuals is infected);
%this implies the posterior distribution for theta is beta(2,40)
y=0;
a1=a0+y;
b1=b0+n-y;
po=betarnd(a1,b1,R,1);
po=betarnd(2,40,R,1);

%plot prior vs. posterior
[f1,x1]=ksdensity(pr,'kernel','epanechnikov','npoints',200,'support',[0 1]);
[f2,x2]=ksdensity(po,'kernel','epanechnikov','npoints',200,'support',[0 1]);

figure(2)
plot(x1,f1,'-b',x2,f2,'-k','LineWidth',1);
title('prior and posterior distributions');
xlabel('$\theta$','interpreter','latex','fontsize',12);
ylabel('density');
h=legend('$p\left(\theta\right)$','$p\left(\theta|y\right)$');
set(h,'interpreter','latex','fontsize',12);

%Send prior and posterior moments to output file
fprintf(fid,'sample size =\t%6.0f \n',n);
fprintf(fid,'Observed number of infections =\t%6.0f \n',y);
fprintf(fid,'\n');

out=[mean(pr) std(pr) min(pr) max(pr)];
fprintf(fid,'Prior moments \n\n');
fprintf(fid,'mean\t\tstd\t\tmin\t\tmax\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out');
fprintf(fid,'\n');

out=[mean(po) std(po) min(po) max(po)];
fprintf(fid,'Posterior moments \n\n');
fprintf(fid,'mean\t\tstd\t\tmin\t\tmax\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out');
fprintf(fid,'\n');


%stop timer
finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);

%save your numerical output
save c:\klaus\AAEC6564\mlab\worksp\mod1_1a;

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




