function [log_predictive,ESS,var_coef]=Standard(Dataset,valid_set,c_0,constr_c,K,prior_intrcp,c1,thita,N,P,num_iter,num_burn,num_thin,varargin)

%STANDARD           The K-fold   distribution \pi(\theta_\gamma,\gamma|y_{-valid_set},c_0) is used as the importance sampling density.  The Holmes 
%                   and Held algorithm described in the article "Bayesian auxiliary  variable models for binary and multinomial regression"
%                   ( Bayesian Analysis (2006)) is used to sample from the importance density.


%Last Modification: 29 March 2009
%%%%%%%%%%%%%%%%
%INPUT variables
%%%%%%%%%%%%%%%%
%Dataset:two matrix
%       n x 1 is the response variable
%       n x p contains the expession levels of all the genes.This matrix
%       has been centered.
%valid_set : cell array
%           Determine the cross-validation set.
%c_0 : matrix 1 x 1
%    The default value of c_0 that resuls in the K-fold standard importance sampler.
%constr_c : matrix 1 x 2
%           The lower and upper bound of the prior variance of the regression coefficients.
%K : matrix 1 x 1
%   Number of different values of c.
%Prior_intrcp : character
%                'Proper' indicates to use a normal prior for the inercept.
%                'Improper' indicates to use the non-informative improper
%                prior for a location parameter with constant density.
%c1 : matrix 1 x 1
%   Prior variance of the constant term.
%thita : matrix 1 x 1
%      The prior mean of the model size.
%N : matrix 1 x 1
%  Maximum number of proposed new variables.
%P : matrix 1 x 1
%  The probability parameter of the binomial distribution.
%num_iter : matrix 1 x 1
%         Number of iteration of the MCMC.
%num_burn : matrix 1 x 1
%         Burn in period of the MCMC.
%num_thin : matrix 1 x 1
%         Thining of the MCMC.
%Optional Input:
%     numGenes : matrix 1 x 1.
%     Number of genes included in the  design matrix.
%
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
%OUTPUT variables
%%%%%%%%%%%%%%%%%
%log_predictive : matrix m x K
%               Log of Cross-validation estimates of the validation set at
%               different values of c. 
%ESS : matrix 1 x K
%    Effective sample size of the importance sampling method.
%var_coef : matrix 1 x K
%          The different values of c (the prior variance of the regression
%          coeficients).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


target=Dataset(:,1);
X=Dataset;X(:,1)=[];

y_vs=target(valid_set);
X_vs=X(valid_set,:);

target(valid_set)=[];
X(valid_set,:)=[];

[n,p]=size(X);


if (nargin==14 & p~=varargin{1})
    numGenes=varargin{1};
    numClass=max(target)+1;
    method1='Dudoit';
    newX=select_genes(X,target,numClass,numGenes,method1);
    X=newX.GENES;
    [n,p]=size(X);

end

c_width=repmat((log(constr_c(2))-log(constr_c(1)))/(K-1),1,K);c_width(1)=0;
var_coef=exp(log(constr_c(1))+cumsum(c_width));%Log-predictive score is estimated at this 12 points.
c=c_0;

%Initial model.
thita=thita/p;
model=binornd(1,2*thita,1,p);
l=find(model);
k=length(l);
t=ones(n,1);

%Model Probability
prob=k*log(thita)+(p-k)*log(1-thita);

%Initial values of the auxiliary variable.
Y(find(target==1),1)=2;Y(find(target==0),1)=-2; % Initial value of Y.

%Initialise the vectors. 
num=(num_iter-num_burn)/num_thin;
numerator=zeros(length(valid_set),length(var_coef));
denom=zeros(length(valid_set),length(var_coef));

%Initial values of intercept and regression coefficients.

if isequal(prior_intrcp,'Proper')

    data=[t X(:,l)];
    H=diag([c1 c*ones(1,k)]);
    K=data'*data+diag([1/c1 1/c*ones(1,k)]);%Precision matrix.
    var=inv(K);%Covariance matrix.

    mu=var*data'*Y;
    %beta=mvnrnd(mu,var)';% Initial sample  of regression coefficients.
    beta=mu+chol(var)'*randn(k+1,1);


else
    %For improper prior.
    data=X(:,l);
    H=diag(c*ones(1,k));
    K=data'*data+diag(1/c*ones(1,k));%Precision matrix.
    var=inv(K);%Covariance matrix.

    mu=var*data'*Y;
    beta=mu+chol(var)'*randn(k,1);
    intrcp=mean(Y)+1/(sqrt(n))*randn;
    beta=[intrcp;beta];

end

%Start the Reversible Jump MCMC algorithm.

for i=1:num_iter


    %Between model move using Reversible Jump and Holmes and Held method.
    %Propose new gamma.

    proposemodel=model;
    w=1+binornd(N,P); %Block size.
    bin_prob=binopdf(w-1,N,P);

    if k<w
        move=2;
    elseif k==w
        move=1+unidrnd(2);
    else
        move=unidrnd(3);
    end

    if move==1
        %Deletion move

        test=0;

        while test~=w %Select w dinstinct predictors.
            r=unidrnd(k,1,w);
            test=length(unique(r));
        end

        proposemodel(l(r))=0;
        newl=find(proposemodel);
        newk=k-w;
        newprob=newk*log(thita)+(p-newk)*log(1-thita);

        %Probability to move from model to proposemodel.
        r12=1/(3*combination(k,w))*bin_prob;

        %Probability to move from proposemodel to model.
        if newk==w

            r21=1/(2*combination(p-(k-w),w))*bin_prob;

        elseif newk<w

            r21=1/(combination(p-(k-w),w))*bin_prob;

        else

            r21=1/(3*combination(p-(k-w),w))*bin_prob;

        end


        if isequal(prior_intrcp,'Proper')

            newdata=[t X(:,newl)];
            newH=diag([c1 c*ones(1,newk)]);
            K1=newdata'*newdata+diag([1/c1 1/c*ones(1,newk)]);%Precision matrix.
            newvar=inv(K1);%Covariance matrix.

        else
            %For improper prior.
            newdata=X(:,newl);
            newH=diag(c*ones(1,newk));
            K1=newdata'*newdata+diag(1/c*ones(1,newk));%Precision matrix.
            newvar=inv(K1);%Covariance matrix.

        end


        Qmodel=Y'*(data*var*data')*Y;
        Qpropose=Y'*(newdata*newvar*newdata')*Y;

        accept=newprob-prob+0.5*(log(det(eye(n)+data*H*data'))-log(det(eye(n)+newdata*newH*newdata')))+0.5*(Qpropose-Qmodel)+log(r21)-log(r12);


    elseif move==2

        %Adittion move

        covariate=find(model==0);
        k1=length(covariate);
        test=0;

        while test~=w %Select w dinstinct predictors.
            r=unidrnd(k1,1,w);
            test=length(unique(r));
        end

        proposemodel(covariate(r))=1;

        newl=find(proposemodel);
        newk=k+w;
        newprob=newk*log(thita)+(p-newk)*log(1-thita);

        if k<w
            %Probability to move from model to proposemodel.
            r12 = 1/(combination(p-k,w))*bin_prob;


        elseif k==w

            %Proobability to move from model to proposemodel.
            r12 = 1/(2*combination(p-k,w))*bin_prob;

        else
            %Probability to move from model to proposemodel.
            r12 = 1/(3*combination(p-k,w))*bin_prob;

        end

        %Probability to move from proposemodel to model.
        r21=1/(3*combination(k+w,w))*bin_prob;

        if isequal(prior_intrcp,'Proper')

            newdata=[t X(:,newl)];
            newH=diag([c1 c*ones(1,newk)]);
            K1=newdata'*newdata+diag([1/c1 1/c*ones(1,newk)]);%Precision matrix.
            newvar=inv(K1);%Covariance matrix.

        else
            %For improper prior.
            newdata=X(:,newl);
            newH=diag(c*ones(1,newk));
            K1=newdata'*newdata+diag(1/c*ones(1,newk));%Precision matrix.
            newvar=inv(K1);%Covariance matrix.

        end


        Qmodel=Y'*(data*var*data')*Y;
        Qpropose=Y'*(newdata*newvar*newdata')*Y;

        accept=newprob-prob+0.5*(log(det(eye(n)+data*H*data'))-log(det(eye(n)+newdata*newH*newdata')))+0.5*(Qpropose-Qmodel)+log(r21)-log(r12);

    else
        %Swaption move
        test=0;

        while test~=w %Select w dinstinct predictors.
            r=unidrnd(k,1,w);
            test=length(unique(r));
        end

        proposemodel(l(r))=0;

        covariate=find(model==0);
        k1=length(covariate);
        test=0;

        while test~=w %Select w dinstinct predictors.
            rr=unidrnd(k1,1,w);
            test=length(unique(rr));
        end

        proposemodel(covariate(rr))=1;

        newl=find(proposemodel);
        newk=k;
        newprob=prob;

        if isequal(prior_intrcp,'Proper')

            newdata=[t X(:,newl)];
            newH=diag([c1 c*ones(1,newk)]);
            K1=newdata'*newdata+diag([1/c1 1/c*ones(1,newk)]);%Precision matrix.
            newvar=inv(K1);%Covariance matrix.

        else
            %For improper prior.
            newdata=X(:,newl);
            newH=diag(c*ones(1,newk));
            K1=newdata'*newdata+diag(1/c*ones(1,newk));%Precision matrix.
            newvar=inv(K1);%Covariance matrix.

        end


        Qmodel=Y'*(data*var*data')*Y;
        Qpropose=Y'*(newdata*newvar*newdata')*Y;

        accept=+0.5*(log(det(eye(n)+data*H*data'))-log(det(eye(n)+newdata*newH*newdata')))+0.5*(Qpropose-Qmodel);


    end

    accept=exp(accept);
    u=rand;

    if u<=accept



        model=proposemodel;
        data=newdata;
        prob=newprob;
        k=newk;
        l=newl;

        H=newH;
        var=newvar;%Covariance matrix.


        %Simulate the new betas from its full conditional given a
        %change to the covariate set.
        if isequal(prior_intrcp,'Proper')

            mu=var*data'*Y;
            beta=mu+chol(var)'*randn(k+1,1);

        else

            %For improper prior.
            mu=var*data'*Y;
            beta=mu+chol(var)'*randn(k,1);
            intrcp=mean(Y)+1/(sqrt(n))*randn;
            beta=[intrcp;beta];

        end

    end


    %Impute  the new value of Y from its full conditional given the model
    %and its parameters.
    if isequal(prior_intrcp,'Proper')

        m=data*beta;

    else

        m=[t data]*beta;

    end

    for j=1:n
        if target(j)==1
            Y(j)=test_truncnorm1(m(j),1,0,+inf);
        else
            Y(j)=test_truncnorm1(m(j),1,-inf,0);
        end
    end

    %Within model move.

    if isequal(prior_intrcp,'Proper')

        mu=var*data'*Y;
        beta=mu+chol(var)'*randn(k+1,1);

    else

        %For improper prior.
        mu=var*data'*Y;
        beta=mu+chol(var)'*randn(k,1);
        intrcp=mean(Y)+1/(sqrt(n))*randn;
        beta=[intrcp;beta];

    end


 

    if i<num_burn,continue,end

    if i==num_burn

        beta1=beta(2:(k+1));

        Weight=(c./var_coef).^(k/2).*exp(-0.5*(beta1'*beta1)*(1./var_coef-1/c));

        if length(valid_set)==1

            numerator=numerator+normcdf([1 X_vs(l)]*beta)^y_vs*(normcdf(-([1 X_vs(l)]*beta)))^(1-y_vs)*Weight;

        else

            numerator=numerator+repmat((normcdf([ones(length(valid_set),1) X_vs(:,l)]*beta).^y_vs.*(normcdf(-([ones(length(valid_set),1) X_vs(:,l)]*beta))).^(1-y_vs)),1,length(var_coef)).*repmat(Weight,length(valid_set),1);

        end

        denom=denom+repmat(Weight,length(valid_set),1);
        dlmwrite('Weight.txt',Weight,'delimiter','\t')


    elseif (mod(i,num_thin)==0)

        beta1=beta(2:(k+1));
        Weight=(c./var_coef).^(k/2).*exp(-0.5*(beta1'*beta1)*(1./var_coef-1/c));

        if (all(isnan(Weight)==0)==1)

            if length(valid_set)==1

                numerator=numerator+normcdf([1 X_vs(l)]*beta)^y_vs*(normcdf(-([1 X_vs(l)]*beta)))^(1-y_vs)*Weight;

            else

                numerator=numerator+repmat((normcdf([ones(length(valid_set),1) X_vs(:,l)]*beta).^y_vs.*(normcdf(-([ones(length(valid_set),1) X_vs(:,l)]*beta))).^(1-y_vs)),1,length(var_coef)).*repmat(Weight,length(valid_set),1);

            end

            denom=denom+repmat(Weight,length(valid_set),1);

            dlmwrite('Weight.txt',Weight,'delimiter','\t','-append')


        end

    end
end



log_predictive=zeros(length(valid_set),length(var_coef));

for j=1:length(valid_set)

    K=max([sum(isinf(log(numerator(j,:)))) sum(isinf(log(denom(j,:))))]);

    if (K==0)

        log_predictive(j,:)=log(numerator(j,:))-log(denom(j,:));

    else

        log_predictive(j,(1:K))=log(numerator(j,K+1))-log(denom(j,K+1));
        log_predictive(j,(K+1):length(var_coef))=log(numerator(j,(K+1):length(var_coef)))-log(denom(j,(K+1):length(var_coef)));

    end

end

Weight=dlmread('Weight.txt');

bar_weight=mean(Weight);
std_weight=std(Weight);

K=sum((bar_weight.^2)==0);
bar_weight(1:K)=bar_weight(K+1);
std_weight(1:K)=std_weight(K+1);

coef_var=((std_weight).^2)./((bar_weight).^2);
ESS=num./(1+coef_var);


