%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%        RUN CODE OF THE PAPER                        %%%%%%%%%%%%%
%%%%%%%%% "Cross-validation prior choice in Bayesian probit   %%%%%%%%%%%%%
%%%%%%%%%          regression with many cocariates"           %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Last Modification: 29 March  2009


%Change the directory according to the folder you have saved the Suplementary_Material.
%In the case you want to run the code from another directory use the
%command addpath.
cd 'C:\Supplementary_Material'

cd 'Importance_Samplers'
addpath '..\Matlib'

%Load the Dataset in the case you are using one of the provided datasets.
Dataset='Arthritis';
load (['..\Data\',Dataset,'.mat'])


%Preprocess the design matrix using the PreProcess_Genes function in the case you are using your own dataset.
%The following specific format of the function is only centering the columns of the design matrix.
%X=PreProcess_Genes(1:length(target),X,target,[0 0 0 0 1]);

%Create randomly the K-fold cross-validation sets using the cross_valid_set
%function in the case you are using your own dataset.
%valid_set=cross_valid_set(length(target),k);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                                                                  %%%%%
%%%%             K-FOLD STANDARD IMPORTANCE SAMPLER                   %%%%%
%%%%                                                                  %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Input the parameters of the standard importance sampler.
c_0=100;constr_c=[0.1 1000];K=12;
prior_intrcp='Proper';
h=100;w=5;
N=0;P=0;
num_iter=2000000;num_burn=100000;num_thin=10;num=(num_iter-num_burn)/num_thin;
index=0;
S=[];
%Run the standard importance sampler.
for j=1:length(valid_set)

    [log_predictive,ESS_1,var_coef]=Standard([target X],valid_set{j},c_0,constr_c,K,prior_intrcp,h,w,N,P,num_iter,num_burn,num_thin);

    S((index+1):(index+length(valid_set{j})),:)=log_predictive;
    ESS(j,:)=ESS_1;
    index=index+length(valid_set{j});

end

Score=-sum(S)/size(X,1);
Effective_SS=mean(ESS);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%           PLOT THE STANDARD IMPORTANCE ESTIMATES OF             %%%%%%
%%%%                    THE LOG-PREDICTIVE SCORE                     %%%%%%
%%%%                                                                 %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pp=csaps(log(var_coef),Score')

figure(1)
fnplt(pp,'-b')
hold on

figure(1)
xlabel('log(c)','Fontsize',12)
ylabel('S(c)','Fontsize',12)
title (['Standard importance sampler (c_0=',num2str(c_0),')'],'Fontsize',12)


%Plot the Effective Sample Size of the importance sampler.
pp=csaps(log(var_coef),log(Effective_SS'))

figure(2)
fnplt(pp,'-b')
hold on

figure(2)
xlabel('log(c)','Fontsize',12)
ylabel('Log ESS','Fontsize',12)
title (['Standard importance sampler (c_0=',num2str(c_0),')'],'Fontsize',12)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                                                                  %%%%%
%%%%             K-FOLD MIXTURE IMPORTANCE SAMPLER                     %%%%%
%%%%                                                                  %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Input parameters of the mixture importance sampler.
constr_c=[0.1 1000];K=[12 20];
prior_intrcp='Proper';
h=100;w=5;
N=0;P=0;
num_iter=200000;num_burn=50000;num_thin=5;num=(num_iter-num_burn)/num_thin;
S=cell(1,length(valid_set));

%Run the standard importance sampler.
for j=1:length(valid_set)

    [log_predictive,var_coef,dinst_c0]=Mixture([target X],valid_set{j},constr_c,K,prior_intrcp,h,w,N,P,num_iter,num_burn,num_thin);
    S{j}=log_predictive;
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%           PLOT THE MIXTURE IMPORTANCE ESTIMATES OF              %%%%%%
%%%%                    THE LOG-PREDICTIVE SCORE                     %%%%%%
%%%%                                                                 %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear log_predictive
lambda=0.5;
difference=normpdf(abs(repmat(log(var_coef)',1,K(2))-repmat(log(dinst_c0),K(1),1))./lambda);
s=size(X);
index=0;

for i=1:length(valid_set)

    if length(valid_set)==s(1)

        log_predictive(i,:)=sum(S{i}.*difference,2)'./sum(difference,2)';

    else

        for k=1:length(valid_set{j})


            log_predictive(index+k,:)=sum(S{i}(:,:,k).*difference,2)'./sum(difference,2)';

        end

        index=index+length(valid_set{i});
        
    end
end


Score=-sum(log_predictive)/size(X,1);

pp=csaps(log(var_coef),Score')

figure(1)
fnplt(pp,'-b')

xlabel('log(c)','FontSize',12)
ylabel('S(c)','FontSize',12)
title('Mixture importance sampler','FontSize',12)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                                                                  %%%%%
%%%%             K-FOLD AUXILIARY IMPORTANCE SAMPLER                  %%%%%
%%%%                                                                  %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Imput parameters of the auxiliary importance sampler.
aux_variable='IGamma';
kappa=0.001;thu=50;delta(1)=kappa;delta(2)=1./thu;delta(3)=0;
constr_c=[0.1 1000];K=12;
prior_intrcp='Proper';
h=100;w=5;
N=0;P=0;
num_iter=20000;num_burn=1000;num_thin=5;num=(num_iter-num_burn)/num_thin;
S=[];
%Run the auxiliary importance sampler.
index=0;

for j=1:length(valid_set)

    [log_predictive,ESS_1,var_coef]=Auxiliary([target X],valid_set{j},aux_variable,delta,constr_c,K,prior_intrcp,h,w,N,P,num_iter,num_burn,num_thin);

    S((index+1):(index+length(valid_set{j})),:)=log_predictive;
    ESS(j,:)=ESS_1;
    index=index+length(valid_set{j});


end

Score=-sum(S)/size(X,1);
Effective_SS=mean(ESS);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%           PLOT THE AUXILIARY IMPORTANCE ESTIMATES OF            %%%%%%
%%%%                    THE LOG-PREDICTIVE SCORE                     %%%%%%
%%%%                                                                 %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pp=csaps(log(var_coef),Score')

figure(1)
fnplt(pp,'-b')
hold on

figure(1)
xlabel('log(c)','Fontsize',12)
ylabel('S(c)','Fontsize',12)
title (['IG(a=',num2str(delta(1)),', b=',num2str(delta(2)),') auxiliary distribution on c'],'Fontsize',12)


%Plot the Effective Sample Size of the importance sampler.
pp=csaps(log(var_coef),log(Effective_SS'))

figure(2)
fnplt(pp,'-b')
hold on

figure(2)
xlabel('log(c)','Fontsize',12)
ylabel('Log ESS','Fontsize',12)
title (['IG(a=',num2str(delta(1)),', b=',num2str(delta(2)),') auxiliary distribution on c'],'Fontsize',12)
