% Let's take a closer look at the prior and posterior distributions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

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

R=10000;

load c:\klaus\AAEC6564\mlab\worksp\mod2_wetlands;

k=size(betamat,1);
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta

bmatprior=zeros(k,R);
for i=1:R
    bmatprior(:,i)=mvnrnd(mu0,V0)';
end

bprior1=bmatprior;
bpost1=betamat;


load c:\klaus\AAEC6564\mlab\worksp\mod2_wetlands2;

k=size(betamat,1);
mu0=[ 0     0.9413       0   0.0151]'; 
int=[100     0.4046^2    100   0.0065^2];
V0=diag(int);

bmatprior=zeros(k,R);
for i=1:R
    bmatprior(:,i)=mvnrnd(mu0,V0)';
end

bprior2=bmatprior;
bpost2=betamat;

%Compare priors and posteriors for log(income) and users
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

incprior1=bprior1(2,:);
incpost1=bpost1(2,:);
incprior2=bprior2(2,:);
incpost2=bpost2(2,:);

useprior1=bprior1(3,:);
usepost1=bpost1(3,:);
useprior2=bprior2(3,:);
usepost2=bpost2(3,:);


[f1,x1]=ksdensity(incprior1,'kernel','epanechnikov','npoints',1000);
[f2,x2]=ksdensity(incpost1,'kernel','epanechnikov','npoints',1000);
[f3,x3]=ksdensity(incprior2,'kernel','epanechnikov','npoints',1000);
[f4,x4]=ksdensity(incpost2,'kernel','epanechnikov','npoints',1000);

[f5,x5]=ksdensity(useprior1,'kernel','epanechnikov','npoints',1000);
[f6,x6]=ksdensity(usepost1,'kernel','epanechnikov','npoints',1000);
[f7,x7]=ksdensity(useprior2,'kernel','epanechnikov','npoints',1000);
[f8,x8]=ksdensity(usepost2,'kernel','epanechnikov','npoints',1000);


figure(1);
subplot(2,1,1);
xl=-2;
xu=4;
yl=0;
yu=1.4;
plot(x1,f1,':k',x2,f2,'-k');
title('prior and posterior for log(inc), original model');
legend('prior','posterior');
axis([xl xu yl yu]);

subplot(2,1,2);
xl=-2;
xu=4;
yl=0;
yu=1.4;
plot(x3,f3,':k',x4,f4,'-k');
title('prior and posterior for log(inc), refined model');
legend('prior','posterior');
axis([xl xu yl yu]);

figure(2);
subplot(2,1,1);
xl=-0.4;
xu=0.4;
yl=0;
yu=20;
plot(x5,f5,':k',x6,f6,'-k');
title('prior and posterior for users, original model');
legend('prior','posterior');
axis([xl xu yl yu]);

subplot(2,1,2);
xl=-0.4;
xu=0.4;
yl=0;
yu=20;
plot(x7,f7,':k',x8,f8,'-k');
title('prior and posterior for users, refined model');
legend('prior','posterior');
axis([xl xu yl yu]);

