
clear
tic
RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(clock)))
niter =50000;
y=36;  N=5766;
stepvar=.004;
m=10;
betas=zeros(niter,m);
betas(1,:)=(1:10)/11;
for j=1:m
    iter=1;
    log_alpha_bot=(y*log(betas(iter,j))+(N-y)*log(1-betas(iter,j))+...
        log(2-2*betas(iter,j)));
    for  iter=2:niter
        X=unifrnd(betas(iter-1,j)-stepvar,betas(iter-1,j)+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))
            betas(iter,j)=X;
            log_alpha_bot=log_alpha_top;
        else
            betas(iter,j)=betas(iter-1,j);
        end
    end
end
toc
plot(betas)
 
 
%%
addpath('/Volumes/iamdavecampbell/CompStatsToolboxV2')
open '/Volumes/iamdavecampbell/CompStatsToolboxV2/Contents'

size(betas)  
R=csgelrub(betas')
% using the above parfor loop we need to do csgelrub(cell2mat(beta')')


