%%  Sampling importance re-sampling 

cd('/Volumes/iamdavecampbell/stat890-2009')
load('stat_890_week1_number10.mat')

iter=iter-1;% because I stopped it the 'iter' value is incomplete so I discard it here


SIS_samples=randsample(1:iter,50000,true,exp(log_alpha_denom{end}(1:iter)));
%the index of values is 1:iter 
% I couldn't use the actual samples because they are multidimensional.
% They're defined by the index anyways...
% I want 50000 new samples (with replacement)
%true <- says that I am including weights
% the last entry are the weights


for lp=1:5;subplot(2,3,lp);plot(pars(lp,SIS_samples),'.');hold on, q=quantile(pars(lp,SIS_samples),[.025,.975]);hold on;plot([1,50000],[q(1),q(1)],'r',[1,50000],[q(2),q(2)],'r');hold off;
    end;subplot(2,3,6);q=quantile(log_alpha_denom(SIS_samples),.05);plot(1:50000,log_alpha_denom(SIS_samples),'.',[1,50000],[q,q],'r');end

figure;% make a 3D plot to see what the joint posterior looks like
plot3(pars(2,SIS_samples),pars(3,SIS_samples),log_alpha_denom(SIS_samples),'.')
xlabel('par 2')
ylabel('par 3')
zlabel('log posterior')
grid on





%% Old way that we ran m=10 completely parallel MCMC runs
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)


%% We can also do this with a parfor loop
% but we have to slice the variable beta in such a way that it can be
% passed to separate workers.
matlabpool open
clear
tic
RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(clock)))
niter=50000;
y=36;  N=5766;
stepvar=.004;
m=10;
betas=cell(m,1);
betas(:)={zeros(niter,1)};
parfor j=1:m
    betas{j}(1)=m/11;
    iter=1;
    log_alpha_bot=(y*log(betas{j}(iter))+(N-y)*log(1-betas{j}(iter))+...
        log(2-2*betas{j}(iter)));
    for  iter=2:niter
        X=unifrnd(betas{j}(iter-1)-stepvar,betas{j}(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))
            betas{j}(iter)=X;
            log_alpha_bot=log_alpha_top;
        else
            betas{j}(iter)=betas{j}(iter-1);
        end
    end
end
toc
betas = cell2mat(betas(1:2)');
plot(betas')



%%

%%  Parallel Tempering in Parallel on a simplified lipomodel with
%%  parameters gamma1 gamma2 and sigma2.  The initial conditions are
%%  assumed to be known L(0)=100 and B(0)=0

RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(clock)))

matlabpool open
cd('/Volumes/iamdavecampbell/stat890-2009')
clear

time=[.5,1,1.5,2,3,4,5,6,7,8,9,10]';
data=[46.1,25.9,17,12.10,7.22,4.51,3.19,2.4,1.82,1.41,1,.94]';

fn.ode=@lipomodelode_abs; % 2 ODE parameters
fn.like=@lipolike;
fn.prior=@lipoprior_abs;

niter=100000;
m=6;
pars=cell(m,1);  % sliced variable is a cell where each entry is a matrix
pars(:)={[[1, .1,  1 ]',zeros(3,niter)]}; % gamma1, gamma2 and sigma^2

% set up the prior parameters in a matrix called prior_pars;
% first row is for the function space,
% second row is for sigma
prior_pars=[0,100;1,1];

accepts=cell(m,1);
accepts(:)={0}; % fill in each entry of the cell array with a 1x1 matrix of zeros

tau=(1:m)/m; % defines the temperatures in our parallel tempering routine, some adjusting to get proper acceptance between chains will need to be done.

stepvar=cell(m,1);
log_alpha_denom=cell(m,1);
log_alpha_denom(:)={zeros(niter,1)};
log_alpha_denom_temp=cell(m,1);
[loglike,path]=fn.like(data,time,pars{1}(:,1),fn);  % we only need to do this once since initially all chains start at the same place
parfor lp=1:m
   stepvar{lp}=[.1,  .5, .5]'/20/tau(lp);
    log_alpha_denom{lp}(1) = (loglike)*tau(lp)+fn.prior(pars{lp}(:,1),path,prior_pars);
end

swapproposed{1}=zeros(m,m);
swapaccepts{1}=zeros(m,m);
parstemp=cell(m,1);
rate=cell(m,1);
iter=3;


tic
a=iter-1;% I do this so that if I want to stop it partway and check on the progress I can then continue it along by defining 'a' as the start of the loop iterates
for  iter=a:niter
    % propose a value from an easy distribution
    chain1=floor(rand*m)+1; % discrete values between 1 and 10 with equal probability
    chain2=floor(rand*m)+1;
    if chain1==chain2  %then don't propose an exchange move
        chain1=0;
        chain2=0;
    else
        ind=sort([chain1,chain2]);
        chain1=ind(1);
        chain2=ind(2);
    end
    
    
    parfor lp=1:m%update all but the ones involved in the swap proposal
        if lp~= chain1 && lp~=chain2
            X=normrnd(pars{lp}(:,iter-1),stepvar{lp});
            [loglike,path]=fn.like(data,time,X,fn);
            log_alpha_numer = loglike*tau(lp)+fn.prior(X,path,prior_pars);
            % make a decision
            if(rand<exp(log_alpha_numer - log_alpha_denom{lp}(iter-1)))
                log_alpha_denom_temp{lp}=log_alpha_numer;
                accepts{lp}=accepts{lp}+1;
                parstemp{lp}=X;
            else
                parstemp{lp}=pars{lp}(:,iter-1);
                log_alpha_denom_temp{lp}=log_alpha_denom{lp}(iter-1);
            end
        end
    end
    
    
    if chain2~=0 % propose a swap
        swapproposed{end}(chain1,chain2)=swapproposed{end}(chain1,chain2)+1;
        
        [loglike1,path1]=fn.like(data,time,pars{chain2}(:,iter-1),fn);
        [loglike2,path2]=fn.like(data,time,pars{chain1}(:,iter-1),fn);
        log_alpha_numer1 = loglike1*tau(chain1)+...
                    fn.prior(pars{chain2}(:,iter-1),path1,prior_pars);
        log_alpha_numer2 = loglike2*tau(chain2)+...
                    fn.prior(pars{chain1}(:,iter-1),path2,prior_pars);
        
        
        if(rand<exp((log_alpha_numer1+log_alpha_numer2) -...
                (log_alpha_denom{chain1}(iter-1)+log_alpha_denom{chain2}(iter-1))));
            log_alpha_denom_temp{chain1}=log_alpha_numer1;
            log_alpha_denom_temp{chain2}=log_alpha_numer2;
            parstemp{chain1}=pars{chain2}(:,iter-1);
            parstemp{chain2}=pars{chain1}(:,iter-1);
            swapaccepts{end}(chain1,chain2)=swapaccepts{end}(chain1,chain2)+1;
        else
            log_alpha_denom_temp{chain1}=log_alpha_denom{chain1}(iter-1);
            log_alpha_denom_temp{chain2}=log_alpha_denom{chain2}(iter-1);
            parstemp{chain1}=pars{chain1}(:,iter-1);
            parstemp{chain2}=pars{chain2}(:,iter-1);
        end
    end
    
    for lp=1:m
        log_alpha_denom{lp}(iter)=log_alpha_denom_temp{lp};
        pars{lp}(:,iter)=parstemp{lp};
    end
    
    if(rem(iter,5000)==0)
        save('PT lipoprotein_abs')
        for lp=1:m
            rate{lp}=[rate{lp};accepts{lp}/...
                  (5000-sum(swapproposed{end}(lp,:)-sum(swapproposed{end}(:,lp))))];
            accepts{lp}=0;
            if(rate{lp}(end)<.18||rate{lp}(end)>.28 && lp<5)  % less than 5 since I know that the stepvar is good for the target posterior near the mode, I'm counting on the other chains to help move around away from that mode
                stepvar{lp}=stepvar{lp}*rate{lp}(end)/.23;
            end
        end
        swapproposed{end+1}=zeros(m,m);
        swapaccepts{end+1}=zeros(m,m);
        disp(iter)
    end
end
toc
% swapaccepts./swapproposed

for kp=1:m;figure(kp);for lp=1:3;subplot(1,3,lp);plot(pars{kp}(lp,1000:iter-1),'.');hold on, q=quantile(pars{kp}(lp,1000:iter-1),[.025,.975]);hold on;plot([1000,iter],[q(1),q(1)],'r',[1000,iter],[q(2),q(2)],'r');hold off;
    end;end
    

figure;% make a 3D plot to see what the joint posterior looks like
plot3(pars{kp}(2,1000:iter-1),pars{kp}(3,1000:iter-1),log_alpha_denom{kp}(1000:iter-1),'.')
xlabel('par 2')
ylabel('par 3')
zlabel('log posterior')
grid on


% This is certainly mixing a lot better but it's still not clear what
% should be included in interval estimates.  


