RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(clock)))
clear
niter=500000;
beta=zeros(niter,1);
y=36;
N=5766;
stepvar=.004;
accepts=0;
beta(1)=.89;% although perhaps 36/5766 might be better;
iter=1
tic
log_alpha_bot=(y*log(beta(iter))+(N-y)*log(1-beta(iter))+...
    log(2-2*beta(iter)));
for  iter=2:niter
    X=unifrnd(beta(iter-1)-stepvar,beta(iter-1)+stepvar);
    log_alpha_top= y*log(X)+(N-y)*log(1-X)+log(2-2*X);
    if(rand<exp(log_alpha_top-log_alpha_bot))
        beta(iter)=X;
        log_alpha_bot=log_alpha_top;
        accepts=accepts+1;
    else
        beta(iter)=beta(iter-1);
    end
end
acceptance_rate = accepts/niter
toc
%%
plot(beta(1:10000))


%%
addpath('/Volumes/iamdavecampbell/coda Matlab')
quantile(beta(end/2:end),.975)
q=.975;r=.005;s=.95;
raftery(beta,q,r,s)


%output is 
% nburn - the burn in 
% nprec - the number of iterations that it takes to reach the desired precision
% nmin  - the number of iterations it would take if the samples were
%       actually from an independent distribution and not from a 
%       Markov Chain with time dependence
% n = nburn+nprec the number of iterations that we really have to run
%       because we don't have independent samples
% irl =n/nmin  the impact of not having independent samples
%%

quantile(beta(end/2:end),.025)
q=.025;r=.005;s=.95;
raftery(beta,q,r,s)



% now let's re-run this diagnostic after discarding some burn-in
raftery(beta(1000:end),.025,r,s)
raftery(beta(1000:end),.975,r,s)


quantile(beta(1000:end),[.025,.5,.975])
%% Geweke and Raftery-Lewis
coda(beta)


vnames='beta the variable';
info.q=.025;
info.r=.005;
info.s=.95;
info.p1=.2;
info.p2=.5;

 coda(beta,vnames,info)
 coda(beta(1:10000),vnames,info)
 coda(beta(10000:end),vnames,info)

