% Simulated data
%clear; 
%simsamplesize=300; % Sample size of simulation.
% true w is  [0.8301,0.0692,-0.5534,0]
%w_sim=[12,1,-8,0]'; w_sim=w_sim/norm(w_sim); 
% true w is [.6369,.0531,-.4246,.3715,-.4777,0.2123,0,0,0,0]
%w_sim=[12,1,-8,7,-9,4,0,0,0,0]'; w_sim=w_sim/norm(w_sim); 
% Generate data by the link function % f(si)=(si^2)exp(si), with single index "si"
%Xsim=rand(simsamplesize,length(w_sim)); 
% Standardize predictors.
%for j=1:length(w_sim);
%if j==1; StandCovs=[]; end;
%StandCovs=[StandCovs,(Xsim(:,j)-mean(Xsim(:,j)))/std(Xsim(:,j),1)];
%end;
%si_sim=sum(Xsim.*repmat(w_sim,1,simsamplesize)',2);
%
% Error distribution.
%errors=0; 
%errors=normrnd(0,.1,simsamplesize,1); % error distribution
%errors=normrnd(0,.2,simsamplesize,1); % error distribution
%errors=normrnd(0,.4,simsamplesize,1); % error distribution
%Linksim=(si_sim.^2).*exp(si_sim); 
%DATA=[Linksim+errors,Xsim]; Outcome=1; Covariates=2:(length(w_sim)+1);
%Labels={'X1','X2','X3','X4'};
%Labels={'X1','X2','X3','X4','X5','X6','X7','X8','X9','X10'};
%s=0; MaxKnots=15;  % Maximum number of knots.
%NumberSamples=75000; thin=5; % Number and thinning of MCMC samples.
%Heteroscedastic=1; %Indicate whether you want heteroscedastic model (0 for homoscedastic)
%delta1=1; delta2=.5; % gamma prior parameters for v


% Example 1: UNESCO data for years 2000-2004.
% Analysing relationship between 
% country educational spending and its GDP.
%for j=1:4;
%if j==1; StandCovs=[]; end;
%StandCovs=[StandCovs,(data(:,j)-mean(data(:,j)))/std(data(:,j),1)];
%end;
%DATA=[StandCovs,data(:,5)]; Covariates=1:4; Outcome=5;
%Labels=textdata(1,2:5);
%s=0; MaxKnots=15;  % Maximum number of knots.
%NumberSamples=75000; % Number and thinning of MCMC samples.
%Heteroscedastic=0; %Indicate whether you want heteroscedastic model (0 for homoscedastic)
%delta1=1; delta2=1.5; % gamma prior parameters for v

% Example 2: Air pollution data.
%for j=3:5;
%if j==3; StandCovs=[]; end;
%StandCovs=[StandCovs,(data(:,j)-mean(data(:,j)))/std(data(:,j),1)];
%end;
%DATA=[nthroot(data(:,2),3),StandCovs]; 
%Outcome=1; Covariates=2:4; 
%Labels=textdata(1,3:5);
%s=0; MaxKnots=15;  % Maximum number of knots.
%NumberSamples=75000; thin=5; % Number and thinning of MCMC samples.
%Heteroscedastic=0; %Indicate whether you want heteroscedastic model (0 for homoscedastic)
%delta1=1; delta2=.25; % gamma prior parameters for v


% Example 3: Rock Data (image processing).
%clear;
%for j=1:3;
%if j==1; StandCovs=[]; end;
%StandCovs=[StandCovs,(data(:,j)-mean(data(:,j)))/std(data(:,j),1)];
%end;
%DATA=[StandCovs,log(data(:,4))];
%clear StandCovs j;
%Covariates=1:3; Outcome=4;   % Outcome is log permeability.
%Labels=[colheaders(:,1:3)];
%s=0; MaxKnots=15;  % Maximum number of knots.
%NumberSamples=75000; % Number and thinning of MCMC samples.
%Heteroscedastic=1; %Indicate whether you want heteroscedastic model (0 for homoscedastic)
%delta1=1; delta2=.5; % gamma prior parameters for v
%plot((0:100)./10,gampdf((0:100)./10,delta1,delta2))

% CPS High school data 2006-07
% DV: Predomintely black school indicator
for j=[7,8,10,13];
if j==7; StandCovs=[]; end;
StandCovs=[StandCovs,(data(:,j)-mean(data(:,j)))/std(data(:,j),1)];
end;
clear j ans 
DATA=[StandCovs,data(:,5)];
Covariates=1:4; Outcome=5;
MaxCategories=repmat(1,length(DATA(:,5)),1);
Labels=textdata(1,[7,8,10,13]);
s=0; MaxKnots=25;  % Maximum number of knots.
NumberSamples=5000; % Number and thinning of MCMC samples.
delta1=1; delta2=1/2; % gamma prior parameters for v

%
%========================================================================================================
% MCMC SAMPLING LOOP.
%========================================================================================================
for s=1:NumberSamples
%========================================================================================================    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if s==1;  % SET UP ALL INITIAL VALUES FOR THE MCMC CHAIN.
format compact;
warning off all;
% Set up objects that describe basic characteristics of data.
ntot=length(DATA(:,Outcome)); 
n=sum(sum(DATA(:,[Covariates,Outcome]),2)>-Inf);
Y=DATA(:,Outcome); X=DATA(:,Covariates); 
p=length(Covariates); 
if (sum(strcmp(who,'MaxCategories'))==1); 
Ylow=NaN(length(Y),1); Ylow(Y==0)=-Inf; Ylow(Y>0 & Y<MaxCategories)=Y(Y>0 & Y<MaxCategories)-1; Ylow(Y==MaxCategories)=Y(Y==MaxCategories)-1; %Replaced -Inf with -100
Yhi=NaN(length(Y),1); Yhi(Y==0)=0; Yhi(Y>0 & Y<MaxCategories)=Y(Y>0 & Y<MaxCategories); Yhi(Y==MaxCategories)=Inf; %Replaced Inf with 100
end;
%
% Starting values for model parameters.
% Starting values of index vector using the estimator of Hristache Juditsky & Spokoiny (2001, Ann Stat).
C0=1; h_k=C0*(n^(-1/max(4,p))); h_max=C0;
a_rho=exp(-1/6); a_h=exp(1/(2*(max(4,p))));
rho_min=n^(-1/3); rho_k=1; betahat=zeros(p,1); k=1;
while rho_k>rho_min;
Sk=(eye(p)+((rho_k^-2)*(betahat*betahat')))^.5;
derivatives=zeros(p+1,1);
for i=1:n;
Term1=zeros(p+1,p+1); Term2=zeros(p+1,1);
    for j=1:n
        Xij=X(j,:)'-X(i,:)';
        t=(norm(Sk*Xij)^2)/(h_k^2);
        K=max(0,(1-(t^2))^2);  % Quartic kernel
        Term1=Term1+([1;Xij]*[1;Xij]'*K);
        Term2=Term2+(Y(j)*[1;Xij]*K);
    end;
derivatives=derivatives+(inv(Term1)*Term2);
end;
betahat=(1/n)*derivatives(2:p+1);
h_k=a_h*h_k; rho_k=a_rho*rho_k;
if rho_k>rho_min; k=k+1; end;
end;
w_hat=(betahat/norm(betahat)); 
w0=w_hat; % Starting value for w is the HJS estimate (Hristache, Jutisky, & Spokoiny, 2001, AS)
clear C0 h_k h_max a_rho a_h rho_min rho_k betahat k Sk derivatives i j Xij t K Term1 Term2
si0=X*w0; SI0=repmat(si0,1,MaxKnots); % Projected scores.
KNOTS0=repmat(quantile(si0,[0:(1/(MaxKnots)):((MaxKnots-1)/MaxKnots)]),n,1); % Knot matrix
B0=[ones(n,1),si0,(SI0>KNOTS0).*(SI0-KNOTS0)]; % Basis matrix
if Heteroscedastic==1; sigma0=unifrnd(.4,.5,n,1); else sigma0=.5; end;% regression variances
v0=.75; 
W0=v0*eye(MaxKnots+2); 
W0(1,1)=1; % starting value for variance of Betas.
%
% Starting value of Beta0, using a 
%Gibbs sampler (p. 330 of O'Hagan & Forster, 2004, Advanced Theory of Statistics, Second Ed.
Winv=inv(W0);
if (length(sigma0)>1); Cinv=inv(diag(sigma0)); else Cinv=inv(sigma0*eye(n)); end;
BCinvB=(B0'*Cinv*B0); BCinvY=B0'*Cinv*Y; 
Wstar=inv(Winv+BCinvB); Wstar=tril(Wstar)+tril(Wstar,-1)';
m_star=Wstar*BCinvY; Beta0=m_star;
clear Cinv BCinvB BCinvY Wstar m_star; 
%
% Set starting DP baseline parameters with Gibbs sampler.
if Heteroscedastic==1;
Nsigclus=length(unique(sigma0));  % Number of sigma clusters.
mu0=1/gamrnd(Nsigclus,1/sum(unique(sigma0))); %Starting value of exponential mean, using Gibbs sampler.
alpha0=1; % Starting value of alpha (precision parameter).
end;
t = cputime;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%===========================================================================================================================
% Gibbs sample latent variables for polychotomous model.
%===========================================================================================================================
if (sum(strcmp(who,'MaxCategories'))==1);
Yold=Y;
Y=norminv(normcdf(Ylow,B0*Beta0,sigma0)+(rand(n,1).*(normcdf(Yhi,B0*Beta0,sigma0)-normcdf(Ylow,B0*Beta0,sigma0))),B0*Beta0,sigma0);
Y(Y==-Inf)=Yold(Y==-Inf); Y(Y==Inf)=Yold(Y==Inf); Y(isnan(Y))=Yold(isnan(Y));
end;
%===========================================================================================================================
%
% Adaptive MH sample w, as vector;
if s==1;
        lambdaJ=10;  % parameter of J spherical (proposal) distribution for w.
        lambdaJ_samples=[lambdaJ;NaN*zeros(NumberSamples-1,1)];
        w_accept_prob_samples=[NaN*zeros(NumberSamples,1)];
        w_Samples=[NaN*zeros(NumberSamples,p)];
        EYsamples=[NaN*zeros(NumberSamples,n)]; siSamples=EYsamples;
end;
% Propose w and new basis matrix.
w1=normrnd(2*(lambdaJ^2)*w0,1,p,1); w1=w1/norm(w1); 
si1=X*w1; SI1=repmat(si1,1,MaxKnots);
KNOTS1=repmat(quantile(si1,[0:(1/(MaxKnots)):((MaxKnots-1)/MaxKnots)]),n,1);
B1=[ones(n,1),si1,(SI1>KNOTS1).*(SI1-KNOTS1)];
% Evaluate Likelihood x prior, for w0 and w1.
LP1=sum(log(normpdf(Y,B1*Beta0,sigma0.^.5)));
LP0=sum(log(normpdf(Y,B0*Beta0,sigma0.^.5)));
% Compute and store accpetance probabilities.
accept_prob_w=min(1,exp(LP1-LP0)); 
w_accept_prob_samples(s)=accept_prob_w;
% Update w and basis matrix.
if(rand<accept_prob_w); w0=w1; si0=si1; B0=B1; SI0=SI1; KNOTS0=KNOTS1; end;
w_Samples(s,:)=w0';
% Store adaptive MH parameters.
if s<=20000; lambdaJ=median([0,lambdaJ+((s^(-1/2))*(.234-accept_prob_w)),10^10]); end;
if s<NumberSamples; lambdaJ_samples(s+1)=lambdaJ; end;
clear LP1 LP0 accept_prob_w w1 si1 SI1 KNOTS1 B1;
%=========================================================================================================
%
%==========================================================================
% Gibbs sample Beta0;
% (p. 330 of O'Hagan & Forster, 2004, Advanced Theory of Statistics, Second Ed.)
%==========================================================================
if s==1; BetaSamples=[NaN*zeros(NumberSamples,length(Beta0))]; warning off all; end;
Winv=inv(W0); 
if (length(sigma0)>1); Cinv=inv(diag(sigma0)); else Cinv=inv(sigma0*eye(n)); end;
BCinvB=B0'*Cinv*B0; BCinvY=B0'*Cinv*Y; 
Wstar=inv(Winv+BCinvB); Wstar=tril(Wstar)+tril(Wstar,-1)'; m_star=Wstar*BCinvY; 
if isnan(sum(sum(Wstar)))==0; EIG=eig(Wstar); else EIG=eig(zeros(MaxKnots+2,MaxKnots+2)); end;
if sum(EIG>0)==(MaxKnots+2); BetaProp=mvnrnd(m_star,Wstar)'; else BetaProp=Beta0; end;
if min(normpdf(Y,B0*BetaProp,sigma0.^.5))>0; Beta0=BetaProp; end;
BetaSamples(s,:)=Beta0'; Means=B0*Beta0;
clear Winv Cinv BCinvB BCinvY BetaProp m_star W_star; 
%=========================================================================================================


%========================================================================================================
% Adaptive MH sample hypervariance (v) parameter of Beta
% Taken from Denison, Holmes, Mallick, & Smith, 2002, p.64, "Bayesian
% Methods for Nonlinear Classification and Regression".
%========================================================================================================
if s==1; % Initialization.
    pv_v=.2; vSamples=[NaN*zeros(NumberSamples,1)]; 
    pv_vSamples=vSamples; AcceptRate_vSamples=vSamples;
end;
%v0=gamrnd(delta1+((MaxKnots+2)/2),inv(delta2+((Beta0'*Beta0)/2)));
%vSamples(s)=v0;
%W0=v0*eye(MaxKnots+2); W0(1,1)=10^10;
% Log probability of v0
Like0=log(mvnpdf(Beta0,zeros(MaxKnots+2,1),W0))+log(gampdf(v0,delta1,delta2));
% Log probability of proposed v
Proposal=exp(normrnd(log(v0),sqrt(pv_v)));
W1=Proposal*eye(MaxKnots+2);  W1(1,1)=10;
Like1=log(mvnpdf(Beta0,zeros(MaxKnots+2,1),W1))+log(gampdf(Proposal,delta1,delta2));
if rand<min(1,exp(Like1-Like0)); v0=Proposal; W0=W1; end;
vSamples(s)=v0;
if s<=20000; pv_v=median([10^(-5),pv_v+((s^(-1/2))*(min(1,exp(Like1-Like0))-.44)),10^10]); end;
pv_vSamples(s)=pv_v;
AcceptRate_vSamples(s)=min(1,exp(Like1-Like0)); 
clear W1 Like0 Like1 Proposal
%========================================================================================================


if Heteroscedastic==1;
%========================================================================================================
% Gibbs sample DP precision (alpha) parameter for regression error variances (from Escobar & West,1995).
%========================================================================================================
if s==1; 
    alphaSamples=[NaN*zeros(NumberSamples,1)]; 
    muSamples=[NaN*zeros(NumberSamples,1)];
    a_alpha=.001; b_alpha=.001; % Noninformative shape and scale prior hyperparameters for alpha (Escobar & West, 1994, JASA, p.584)
end;
eta=betarnd(alpha0+1,ntot); Odds=(a_alpha+Nsigclus-1)/(ntot*(b_alpha-log(eta))); 
alpha0=gamrnd(a_alpha+Nsigclus-(rand>(Odds/(1+Odds))),b_alpha-log(eta));
alphaSamples(s)=alpha0; clear eta Odds; 
%========================================================================================================
% Gibbs Sample mean parameter (mu0) of Exponential baseline distribution,
% with mu0 under a reference prior.
%========================================================================================================
if s==1; muSamples=[NaN*zeros(NumberSamples,1)]; end;
mu0=1/gamrnd(Nsigclus,1/sum(unique(sigma0))); muSamples(s)=mu0;
%========================================================================================================
end;

%========================================================================================================
% Sample the regression variances under the DP model, using Algorithm 8 of Neal(2000);
%========================================================================================================
if s==1; %Initialize.
    if Heteroscedastic==1;  
        NsigclusSamples=[NaN*zeros(NumberSamples,1)]; 
        sigmaSamples=[NaN*zeros(NumberSamples,length(Y))];
        m=1;
    end;    
if Heteroscedastic==0; sigmaSamples=[NaN*zeros(NumberSamples,1)]; end;
AcceptRateSigSamples=[NaN*zeros(NumberSamples,1)]; pvSig=.2; pvSigSamples=AcceptRateSigSamples; 
end;
%
if Heteroscedastic==1;
%========================================================================================================
% First, sample the clusters of the sigmas;
tab=tabulate(sigma0);Counts=tab(:,2);COUNTS=repmat(Counts',ntot,1); clear tab Counts;
[SIGMAS,SIGMA0]=meshgrid(unique(sigma0)',sigma0');
IndSIGMA=SIGMA0==SIGMAS;COUNTS=COUNTS-IndSIGMA; clear SIGMA0 IndSIGMA; % Revise counts to n_-i,c.
uniq=sum((COUNTS==0)*1,2); uniq=sort(unique(uniq.*(1:ntot)')); uniq=uniq((2:length(uniq)),:);
RATIO1=COUNTS./(n-1+alpha0); clear COUNTS;
RATIO2=(alpha0/m)/(n-1+alpha0); RATIO2(1:ntot,1:m)=RATIO2;
RATIO=[RATIO1,RATIO2]; clear RATIO1 RATIO2;
SIGMAS=[SIGMAS,exprnd(mu0,ntot,m)];
SIGMAS(uniq,Nsigclus+m)=sigma0(uniq,:); clear uniq; % add unique sigmas as candidate sigmas.
MEANS=repmat(Means',Nsigclus+m,1)'; OUTCOME=repmat(Y',Nsigclus+m,1)';
LIKE=RATIO.*normpdf(OUTCOME,MEANS,SIGMAS.^.5); clear RATIO OUTCOME MEANS;
Norm=sum(LIKE,2); Norm=repmat(Norm',Nsigclus+m,1)'; 
LIKE=cumsum((LIKE./Norm),2); clear Norm;
sample=unifrnd(zeros(ntot,1),LIKE(:,Nsigclus+m));
[sample,cols]=meshgrid(sample,(1:(Nsigclus+m))); sample=sample'; cols=cols';
sample=sample>LIKE; clear LIKE; 
sample=max(sample.*cols,[],2)+1; sample=repmat(sample',Nsigclus+m,1)';
sample=(cols==sample); clear cols;
sigma0=sum(sample.*SIGMAS,2); clear sample SIGMAS;
Nsigclus=length(unique(sigma0)); 
NsigclusSamples(s,:)=Nsigclus;
% Given sample clusters, use adaptive MH step to generate new posterior sample of sigmas;
uniq=unique(sigma0); ClusterIDs=sum((repmat(uniq,1,ntot)'==repmat(sigma0,1,Nsigclus)).*repmat(1:Nsigclus,ntot,1),2);
LP0=log(normpdf(Y,Means,sigma0.^.5))+log(exppdf(sigma0,mu0)); 
LP0=accumarray(ClusterIDs,LP0);
Proposal=exp(normrnd(log(uniq),pvSig)); Proposal_n=Proposal(ClusterIDs,:);
LP1=log(normpdf(Y,Means,Proposal_n.^.5))+log(exppdf(Proposal_n,mu0));
LP1=accumarray(ClusterIDs,LP1);
AcceptProbs=min(1,exp(LP1-LP0)); clear Proposal_n LP1 LP0 i; 
Accept=rand(Nsigclus,1)<AcceptProbs; clear AcceptProbs;
Accepted=(Proposal.*Accept)+(uniq.*(1-Accept)); clear uniq Proposal;
sigma0=Accepted(ClusterIDs,:); AcceptRate=mean(Accept); clear ClusterIDs Accepted Accept;
if (AcceptRate>.441 && s<=20000); pvSig=exp(log(pvSig)+min(.01,s^(-.5))); end;
if (AcceptRate<.441 && s<=20000); pvSig=exp(log(pvSig)-min(.01,s^(-.5))); end;
sigmaSamples(s,:)=sigma0';
AcceptRateSigSamples(s)=AcceptRate; 
pvSigSamples(s)=pvSig; clear AcceptRate;
end;
% Use adaptive metropolis step for homoscedastic model.
if (Heteroscedastic==0 || Discrete==1);
if (Discrete==0); 
Like0=sum(log(normpdf(Y,Means,sigma0^.5)));
Proposal=exp(normrnd(log(sigma0),sqrt(pvSig)));
Like1=sum(log(normpdf(Y,Means,Proposal^.5)));
if rand<min(1,exp(Like1-Like0)); sigma0=Proposal; end;
else sigma0=1;
end;
sigmaSamples(s)=sigma0;
if (s<=20000); pvSig=median([10^(-5),pvSig+((s^(-1/2))*(min(1,exp(Like1-Like0))-.44)),10^10]); end;
pvSigSamples(s)=pvSig;
AcceptRateSigSamples(s)=min(1,exp(Like1-Like0)); 
clear Like0 Like1 Proposal
end;
%========================================================================================================
% Compute log-likelihood of model, given current MCMC state.
if s==1; LikeSamples=[NaN*zeros(NumberSamples,1)]; end;
LikeSamples(s)=sum(log(normpdf(Y,Means,sigma0.^.5)+(10^(-323))));
%
% Construct posterior predictive residuals
% and sample of link function.
if s==1; 
    PredictErrorSamples=[NaN*zeros(NumberSamples,n)]; 
    EYsamples=[NaN*zeros(NumberSamples,n)]; 
    siSamples=[NaN*zeros(NumberSamples,n)]; 
end;
siSamples(s,:)=si0'; EYsamples(s,:)=Means';
Predict=normrnd(Means,sigma0.^.5,n,1);
PredictErrorSamples(s,:)=(Y-Predict)';
clear Predict;
%
% This keeps track of time:
if s/500==round(s/500) 
    Iteration=s; Minutes=(cputime-t)/60; MinutesLeft=(Minutes*(NumberSamples/s))-Minutes;
    display(Iteration); display(Minutes); display(MinutesLeft);
    clear Iteration Minutes;
end;
%
end;
% end of single MCMC sampling loop;






% Convergence analysis from simulations.
thin=5; b=20000<1:75000; b=find(b,1,'first'):find(b,1,'last'); b=[min(b):thin:max(b)];
% Plot true link (black) against posterior mean link (red)
%figure('Color','white'); 
hold on;
subplot(4,2,7);
sorted=sortrows([si_sim';Linksim';1:simsamplesize]',1);
sorted_si_sim=sorted(:,1); sorted_Linksim=sorted(:,2); 
scatter(sorted_si_sim,Y(sorted(:,3)));
hold on;
plot(sorted_si_sim,sorted_Linksim,'Color','Black','LineWidth',2);
xlabel('Single Index','FontSize',11);
ylabel('Expected Response','FontSize',11);
xlim([-.5,1]); ylim([-.5,1.5]);
title('n=200, sigma=.2');
subplot(4,2,8);
plot(sorted_si_sim,sorted_Linksim,'Color','Black','LineWidth',2);
hold on;
sorted_est=sortrows([mean(siSamples(b,:),1); mean(EYsamples(b,:),1);1:simsamplesize]',1);
sorted_si_sim_est=sorted_est(:,1); sorted_Linksim_est=sorted_est(:,2);
plot(sorted_si_sim_est,sorted_Linksim_est,'Color','Red','LineWidth',2);
xlabel('Single Index','FontSize',11);
ylabel('','FontSize',11);
xlim([-.5,1]); ylim([-.5,1.5]);
title('n=200, sigma=.2');



% Comparing true link against posterior mean link:
for j=1:length(b);
    if j==1; LinkError=0; end;
    LinkError=LinkError+max(abs(Linksim'-EYsamples(b(j),:)));
end;
LinkError=LinkError/length(b);
% Comparing posterior mean w against simulated w
for j=1:length(b);
    if j==1; wError=0; end;
    wError=wError+max(abs(w_sim'-w_Samples(b(j),:)));
end;
wError=wError/length(b); % error of w posterior samples
%posterior mean of average sqrt(error variances)
meanSig=mean(sqrt(sigmaSamples(b)));
PostMeanw=mean(w_Samples(b,:)); PostMeanw=PostMeanw/norm(PostMeanw);
ErrorPostMeanw=max(abs(w_sim'-PostMeanw));
ErrorPostMeanLink=max(abs(Linksim'-mean(EYsamples(b,:),1)));
[wError,LinkError,ErrorPostMeanw,ErrorPostMeanLink,meanSig]



% CONVERGENCE PLOTS.
%Exclude burn-in samples in parameter summaries.
figure('Color','white');
% Plot samples of log-likelihood
subplot(3,2,1);
plot(lambdaJ_samples); xlabel('Iteration','FontSize',12); ylabel('Prop. Parameter {\it\lambda}','FontSize',12); set(findobj('Type','line'),'LineWidth',2);
subplot(3,2,2);
plot(cumsum(w_accept_prob_samples)'./[1:length(w_accept_prob_samples)]); 
xlabel('Iteration','FontSize',12); ylabel('Cumul.Acc.Rate of {\it\omega}','FontSize',11);
subplot(3,2,3);
plot(pv_vSamples); xlabel('Iteration','FontSize',12); ylabel('Prop. Var for {\it\nu}','FontSize',12); set(findobj('Type','line'),'LineWidth',2);
subplot(3,2,4);
plot(cumsum(AcceptRate_vSamples)'./[1:length(AcceptRate_vSamples)]); 
xlabel('Iteration','FontSize',12); ylabel('Cumul.Acc.Rate of {\it\nu}','FontSize',11);
subplot(3,2,5);
plot(pvSigSamples); xlabel('Iteration','FontSize',12); 
ylabel('Prop. Var for {\it\sigma^2}','FontSize',12); set(findobj('Type','line'),'LineWidth',2);
subplot(3,2,6);
plot(AcceptRateSigSamples); xlabel('','FontSize',12); ylabel('Accept. Rate','FontSize',12);
xlabel('Iteration','FontSize',12); ylabel('Acceptance Rate','FontSize',12); set(findobj('Type','line'),'LineWidth',2);


% Boxplot of w coefficients, and posterior mean of link function.
%Exclude burn-in samples in parameter summaries.
b=20000<1:s; thin=5; b=find(b,1,'first'):find(b,1,'last'); b=[min(b):thin:max(b)];
figure('Color','white'); 
subplot(4,1,1);
Covs=1:length(Beta0); % Choose covariates you wish to plot
Vals=w_Samples(b,:);
boxplot(Vals,'orientation','horizontal','symbol','');
xlabel('Index Coefficient, {\it\omega}','FontSize',12);
ylabel('Covariate','FontSize',12);
set(gca,'YTickLabel',Labels);
set(findobj('Type','line'),'LineWidth',2);
for j=1:length(Vals(1,:))
text(median(Vals(:,j)),j-.35,num2str(round(mean(Vals(:,j))*100)/100),'fontsize',8,'VerticalAlignment','middle');
text(quantile(Vals(:,j),.025),j+.35,num2str(round(quantile(Vals(:,j),.025)*100)/100),'fontsize',8,'VerticalAlignment','middle','HorizontalAlignment','right');
text(quantile(Vals(:,j),.975),j+.35,num2str(round(quantile(Vals(:,j),.975)*100)/100),'fontsize',8,'VerticalAlignment','middle','HorizontalAlignment','left');
end;
clear j Vals MinBeta;
% Plot posterior mean estimate of link function.
subplot(4,1,2);
sorted_est=sortrows([mean(siSamples(b,:),1); mean(EYsamples(b,:),1);1:n]',1);
sorted_si_sim_est=sorted_est(:,1); sorted_Linksim_est=sorted_est(:,2);
scatter(sorted_si_sim_est,Y(sorted_est(:,3)));
hold on;
plot(sorted_si_sim_est,sorted_Linksim_est,'Color','Red','LineWidth',3);
xlabel('Single Index, {\bfx\omega}','FontSize',11);
ylabel('Expected Response','FontSize',11);
% Plot residuals.
subplot(4,1,3);
for j=1:n
if ((j/2)==round(j/2)); plot([j,j],[quantile(PredictErrorSamples(b,j),.01),quantile(PredictErrorSamples(b,j),.99)],'Color','Blue'); end;
if ((j/2)~=round(j/2)); plot([j,j],[quantile(PredictErrorSamples(b,j),.01),quantile(PredictErrorSamples(b,j),.99)],'Color','Red'); end;
hold on;
end;
plot([1,n]',[0,0]')
xlabel('Individual {\iti}','FontSize',11);
ylabel('Residual (99%)','FontSize',11);
xlim([1,n]);
hold off;
subplot(4,1,4);
Vals=vSamples(b,:);
ksdensity(Vals,'support','positive','width',.20);
xlabel('Coefficient Variance Hyperparameter \it\nu','FontSize',11);
ylabel('Density','FontSize',11);
set(findobj('Type','line'),'LineWidth',2);
xlim([quantile(Vals,.01),quantile(Vals,.99)]);
text(mean(Vals),0,num2str(round(mean(Vals)*100)/100),'VerticalAlignment','bottom','HorizontalAlignment','center','fontsize',8)
text(min(xlim),max(ylim),num2str(round(quantile(Vals,.025)*100)/100),'VerticalAlignment','top','fontsize',8,'HorizontalAlignment','left')
text(max(xlim),max(ylim),num2str(round(quantile(Vals,.975)*100)/100),'VerticalAlignment','top','fontsize',8,'HorizontalAlignment','right')
hold off;

%IN A SINGLE FIGURE, PLOT THE ERROR VARIANCES AND THEIR DP BASELINE
%PARAMETERS.
%Exclude burn-in samples in parameter summaries.
b=20000<1:s; thin=5; b=find(b,1,'first'):find(b,1,'last'); b=[min(b):thin:max(b)];
figure('Color','white'); 
% Plot posterior mean of error variances;
subplot(2,1,1);
scatter(1:n,mean(sigmaSamples(b,:),1),50,'filled')
hold on;
xlabel('Individual {\iti}','FontSize',11);
ylabel('Error Variance {\it\sigma^2}  (Post.Mean)','FontSize',11);
xlim([1,length(sigmaSamples(1,:))]);
% Plot precision parameter;
subplot(2,3,4);
Vals=alphaSamples(b,:);
ksdensity(Vals,'support','positive','width',.05);
xlabel('DP precision \it\alpha','FontSize',11);
ylabel('Density','FontSize',11);
set(findobj('Type','line'),'LineWidth',2);
xlim([quantile(Vals,.01),quantile(Vals,.99)]);
text(mean(Vals),0,num2str(round(mean(Vals)*100)/100),'VerticalAlignment','bottom','HorizontalAlignment','center','fontsize',8)
text(min(xlim),max(ylim),num2str(round(quantile(Vals,.025)*100)/100),'VerticalAlignment','top','fontsize',8,'HorizontalAlignment','left')
text(max(xlim),max(ylim),num2str(round(quantile(Vals,.975)*100)/100),'VerticalAlignment','top','fontsize',8,'HorizontalAlignment','right')
% Plot mean parameters of regression variances;
subplot(2,3,5);
Vals=muSamples(b,:);
ksdensity(Vals,'support','positive','width',.11);
xlabel('Mean {\it\mu} of Error Variances','FontSize',11);
ylabel('','FontSize',11);
set(findobj('Type','line'),'LineWidth',2);
xlim([quantile(Vals,.01),quantile(Vals,.99)]);
text(mean(Vals),0,num2str(round(mean(Vals)*100)/100),'VerticalAlignment','bottom','HorizontalAlignment','center','fontsize',8)
text(min(xlim),max(ylim),num2str(round(quantile(Vals,.025)*100)/100),'VerticalAlignment','top','fontsize',8,'HorizontalAlignment','left')
text(max(xlim),max(ylim),num2str(round(quantile(Vals,.975)*100)/100),'VerticalAlignment','top','fontsize',8,'HorizontalAlignment','right')
% Plot histogram of number of clusters.
subplot(2,3,6);
Vals=NsigclusSamples(b,:);
[Count,bins]=hist(Vals,length(min(Vals):max(Vals)));
bar(round(bins),Count./sum(Count));
xlim([quantile(Vals,.01),quantile(Vals,.99)]);
xlabel('Number of {\it\sigma^2} Clusters','FontSize',11);
ylabel('Probability','FontSize',11);
set(gca,'XTick',min(round(bins))-1:7:(max(round(bins))+1));
text(mean(Vals),max(ylim),num2str(round(mean(Vals))),'VerticalAlignment','top','HorizontalAlignment','center','fontsize',8);
text(min(xlim),max(ylim),num2str(round(quantile(Vals,.025))),'VerticalAlignment','top','HorizontalAlignment','left','fontsize',8);
text(max(xlim),max(ylim),num2str(round(quantile(Vals,.975))),'VerticalAlignment','top','HorizontalAlignment','right','fontsize',8);
hold off;



% Table of posterior summaries.
b=20000<1:s; thin=5; b=find(b,1,'first'):find(b,1,'last'); b=[min(b):thin:max(b)];
wTable=[mean(w_Samples(b,:),1)./norm(mean(w_Samples(b,:),1));quantile(w_Samples(b,:),.025);quantile(w_Samples(b,:),.975)]';
vTable=[mean(vSamples(b,:)); quantile(vSamples(b,:),.025); quantile(vSamples(b,:),.975)]';
DpTable=[mean(alphaSamples(b,:));quantile(alphaSamples(b,:),.025);quantile(alphaSamples(b,:),.975)]';
DpTable=[DpTable;mean(muSamples(b,:)),quantile(muSamples(b,:),.025),quantile(muSamples(b,:),.975)];
DpTable=[DpTable;mean(NsigclusSamples(b,:)),quantile(NsigclusSamples(b,:),.025),quantile(NsigclusSamples(b,:),.975)];
Table=[wTable;vTable;DpTable];
Table
LabelTable=[Labels,'nu','alpha','mu','NumClusters']'



% Plot posterior mean link function.
figure('Color','white');
b=20000<1:s; thin=5; b=find(b,1,'first'):find(b,1,'last'); b=[min(b):thin:max(b)];
% Plot posterior mean estimate of link function.
sorted_est=sortrows([mean(siSamples(b,:),1); mean(EYsamples(b,:),1);1:n]',1);
sorted_si_sim_est=sorted_est(:,1); sorted_Linksim_est=sorted_est(:,2);
scatter(sorted_si_sim_est,Y(sorted_est(:,3)));
hold on;
plot(sorted_si_sim_est,sorted_Linksim_est,'Color','Red','LineWidth',3);
xlabel('Single Index, {\bfx\omega}','FontSize',11);
ylabel('Expected Response','FontSize',11);
hold off;


b=20000<1:s; thin=5; b=find(b,1,'first'):find(b,1,'last'); b=[min(b):thin:max(b)];
figure('Color','white'); 
% Plot residuals.
subplot(2,1,1);
for j=1:n
if ((j/2)==round(j/2)); plot([j,j],[quantile(PredictErrorSamples(b,j),.01),quantile(PredictErrorSamples(b,j),.99)],'Color','Blue'); end;
if ((j/2)~=round(j/2)); plot([j,j],[quantile(PredictErrorSamples(b,j),.01),quantile(PredictErrorSamples(b,j),.99)],'Color','Red'); end;
hold on;
end;
plot([1,n]',[0,0]')
xlabel('Individual {\iti}','FontSize',11);
ylabel('Residual (99%)','FontSize',11);
xlim([1,n]);
hold off;
% Plot posterior mean of error variances;
subplot(2,1,2);
scatter(1:n,mean(sigmaSamples(b,:),1),50,'filled')
hold on;
xlabel('Individual {\iti}','FontSize',11);
ylabel('Error Variance {\it\sigma^2}  (Post.Mean)','FontSize',11);
xlim([1,length(sigmaSamples(1,:))]);
hold off;


% Model Selection Analysis:
% Calculate posterior-mean log-likelihood
%Exclude burn-in samples in parameter summaries.
b=20000<1:s;  thin=5; b=find(b,1,'first'):find(b,1,'last'); b=[min(b):thin:max(b)];
LogLike=mean(LikeSamples(b));
% Compute L2 model selection criterion.
OBS=repmat(Y',length(b),1);
PREDICTIONS=OBS-PredictErrorSamples(b,:);
MeanPredictError=sum((Y'-mean(PREDICTIONS,1)).^2);
VarPredict=sum(var(PREDICTIONS,1,1));
L2=sum(MeanPredictError)+sum(VarPredict);
clear OBS;
ModelSelect=[MeanPredictError,VarPredict,L2]