function [ boot,accept,times_in_model] = human_MCMC_norm( runs, SNP, vfrac,scale)

  load Region18qnhnc2.txt %199p 71 snp
  data= Region18qnhnc2;
  data2=data;
   %Remove unneeded columns
    data2(:,5)=[];
    data2(:,4)=[];
    data2(:,3)=[];
    data2(:,2)=[];
    data2(:,1)=[];
   boot=zeros(100,183);
beta=zeros(183,1);
beta(SNP)=vfrac;
n=1;
while n<51
    
    y=floor((183-1+1)*rand+1);
    prob= exp((data2(y,:)*beta))/(1+exp((data2(y,:)*beta)));
    random=rand;
    if random<prob
        boot(n,:)=data2(y,:);
        n=n+1;
    end
end
m=51;

while m<101
          
          y=floor((183-1+1)*rand+1);
          q=exp((data2(y,:)*beta))/(1+exp((data2(y,:)*beta)));
    prob= 1-q; 
    random=rand;
    if random<prob
        boot(m,:)=data2(y,:);
        m=m+1;
    end
end

  %acceptance counter: matrix where columns are move , top row times chosen,
%bottom row times accepted. Flag 1 add in, 2 remove, 3 update, 4 mv

accept=zeros(2,5);

%genotype matrix
X=boot;
x_dim=size(X);
no_SNPs=x_dim(2);
no_runs=runs;
no_patients=100;
no_cases=0.5*no_patients;
%Starting Beta Values
    % Variables for qtl gibbs step
      a_0=0;
      b_0=1;
      V_0=0;
      m_0=0;
    % Run under prior ( 0 if yes, 1 if not)
      likelihood_flag=1;
    % Prior on number of markers in model
      SNP_prior=1;
      %death rate
      death_rate=0.3;
      switch_rate=0.4;
  %acceptance counter: matrix where columns are move , top row times chosen,
  %bottom row times accepted. Flag 1 add in, 2 remove, 3 update, 4 mv

    accept=zeros(2,5);
    
 % Memory Alocation
     current_beta=zeros(no_SNPs,1);
     perturbed_beta=zeros(no_SNPs,1);
     current_likelihood=0;
     perturbed_likelihood=0;
     beta_intercept=0*ones(no_patients,1);
     linear_predictor = zeros(no_patients,1);
     w = zeros(no_patients,1);
     prediction=ones(no_patients,1)./(1+exp(-(linear_predictor)));
     gamma=zeros(1,no_SNPs);
     beta_indicator=zeros(1,no_SNPs);
     X_gamma=0;
     V_gamma=0;
     beta_gamma=[];
     qvar=0;
     times_in_model=zeros(1,no_SNPs);
     indicator_vector=zeros(1,no_SNPs);
     burn=0.2*no_runs*no_SNPs;
% Variables for prior and noise updates
     prior_scale=(scale);%(2*qtlsim_p*(1-qtlsim_p)/v_fraction)^0.5;
     half=0.5*ones(no_patients,1);
     y=ones(no_patients,1);
     for ddd=51:100
         y(ddd)=0;
     end

% Model
     i=1;
     iteration=1;
     pk=10;
     
while i<=no_runs
    j=1;
    while j<=no_SNPs
        iteration=iteration+1;
    %Variance calculation
    sgamma=sum(gamma);
   if sgamma==0
      V_1=1;
       m_1=0;
       w = zeros(no_patients,1);
       
       for q=1:no_patients
           
           for kk=1:pk
               C=0;
               randgam=exprnd(1);
           w(q)=w(q)+(randgam/((kk-0.5)^2+C^2/(4*pi^2)))/(2*pi);
           end
       end
       z=(y-half)./w;
       omega=diag(w);
   else
       w = zeros(no_patients,1);
       for q=1:no_patients
           for kk=1:pk
               C=X_gamma(q,:)*beta_gamma';
               randgam=exprnd(1);
           w(q)=w(q)+(randgam/((kk-0.5)^2+C^2/(4*pi^2)))/(2*pi^2);
           end
       end
       
   
       z=(y-half)./w;
       omega=diag(w);
            V_1=(X_gamma'*omega*X_gamma+V_gamma^-1)^-1;
            m_1=V_1*(X_gamma'*omega*z);
            
           

            beta_gamma=mvnrnd(m_1,V_1);
             ww=1;
                for w=1:no_SNPs
                    
                    if gamma(w)==0
                        current_beta(w)=0;
                    else
                        current_beta(w)=beta_gamma(ww);
                        ww=ww+1;
                    end
                end
   end
   
            
    
        
        variablenumber=rand;
            
        if(variablenumber<1)
 		
 			thresh_pi=rand;
 	

            
            if ( current_beta(j)~=0 && thresh_pi < death_rate)
                gamma(j) = 0;
                flag=2;

                         
                     
                         
            else 
                gamma(j) = 1;
                 if (current_beta(j)==0)
                    flag=1;
                else
                    flag=3;
                end
            end
            
            X_gamma_1=[];
            beta_gamma_1=[];
            for col=1:no_SNPs
                if gamma(col)==1
                    X_gamma_1=[X_gamma_1 X(:,col)];
                    beta_gamma_1=[beta_gamma_1 current_beta(col)];
                end
            end
            
            %X_gamma_1=X(:,[gamma]);
            xgamma_dim=size(X_gamma_1);
            
            if sum(xgamma_dim)==0
                m_1_gamma=0;
                V_1_gamma=1;
                beta_gamma_1=0;
                w = zeros(no_patients,1);
                for q=1:no_patients
                for kk=1:pk
               C=0;
               randgam=exprnd(1);
                w(q)=w(q)+(randgam/((kk-0.5)^2+C^2/(4*pi^2)))/(2*pi^2);
                end
                end
                z_1=(y-half)./w;
                omega_1=diag(w);
            else
                w = zeros(no_patients,1);
                for q=1:no_patients
           for kk=1:pk
               C=X_gamma_1(q,:)*beta_gamma_1';
               randgam=exprnd(1);
           w(q)=w(q)+(randgam/((kk-0.5)^2+C^2/(4*pi^2)))/(2*pi);
           end
       end
       
   
       z_1=(y-half)./w;
       omega_1=diag(w);
            
           
            V_gamma_1=eye(xgamma_dim(2));
            m_0_gamma=zeros(xgamma_dim(2),1);
            V_1_gamma=(X_gamma_1'*omega_1*X_gamma_1+V_gamma_1^-1)^-1;
            m_1_gamma=V_1_gamma*(X_gamma_1'*omega_1*z_1);
            
            beta_gamma_1=mvnrnd(m_1_gamma,V_1_gamma);
            end
            
            accept(1,flag)= accept(1,flag)+1;
            ww=1;
                for w=1:no_SNPs
                    
                    if gamma(w)==0
                        perturbed_beta(w)=0;
                    else
                        perturbed_beta(w)=beta_gamma_1(ww);
                        ww=ww+1;
                    end
                end
            change_in_log_prior = -(1/(2*prior_scale^2))*sum((perturbed_beta-current_beta).^2);
            %change_in_log_prior = -(1/10)*sum((abs(perturbed_beta)-abs(current_beta)));
            linear_predictor = X*current_beta;
            
            proposed_linear_predictor = X*perturbed_beta;
            
            if likelihood_flag==1
            cs = [prediction(1:no_cases); 1-prediction(no_cases+1:no_patients)];
            current_log_likelihood = sum(log(cs));
           
            %current_log_likelihood=0;
            
            
           
            proposed_prediction=ones(no_patients,1)./(1+exp(-(proposed_linear_predictor)));
            
            ps = [proposed_prediction(1:no_cases); 1-proposed_prediction(no_cases+1:no_patients)];
            perturbed_log_likelihood = sum(log(ps));
            end
            
            if likelihood_flag==0
            current_log_likelihood=0;
            perturbed_log_likelihood=0;
            end
            
		
            
            %calculate ratio=prior raio * likelihood ratio * proposal ratio
                
                if (current_beta(j)==0)
                %q_beta = -sqrt(eta) times proposed value sqaured
                if (sgamma==0)
                    
                    log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior - log(mvnpdf(beta_gamma_1',m_1_gamma,V_1_gamma))+log(death_rate)+0+log(1/(no_SNPs-SNP_prior))-log(1-(1/(no_SNPs-SNP_prior)));
                else
                    log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior - log(mvnpdf(beta_gamma_1',m_1_gamma,V_1_gamma))+log(death_rate)+log(mvnpdf(beta_gamma',m_1,V_1))+log(1/(no_SNPs-SNP_prior))-log(1-(1/(no_SNPs-SNP_prior)));
		 		%proposal ratio= prob.birth*prob.that value of beta/prob.death log(1/182) chance of inclusion
                end
                
                elseif (current_beta(j)~=0 && thresh_pi < death_rate)
                    sgamma=sum(gamma);
                    if (sgamma==0)
                     log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior - log(death_rate)- 0+log(mvnpdf(beta_gamma',m_1,V_1))-log(1/(no_SNPs-SNP_prior))+log(1-(1/(no_SNPs-SNP_prior)));
                     else
                    log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior - log(death_rate)- log(mvnpdf(beta_gamma_1',m_1_gamma,V_1_gamma))+log(mvnpdf(beta_gamma',m_1,V_1))-log(1/(no_SNPs-SNP_prior))+log(1-(1/(no_SNPs-SNP_prior)));
                    end
                %proposal ratio=1/birth proposal ratio
%                 elseif (current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
%                     log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior;
%                 
                else    		 		  
                log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior;
                %proposal ratio=1
              
                end
           
            randomnumber=rand;


            % MONTIOR ACCEPTANCE RATES

            
            if (log_R > 0)
                accept(2,flag)= accept(2,flag)+1;
                X_gamma=X_gamma_1;
                V_gamma=V_gamma_1;
                m_0=m_0_gamma;
                beta_gamma=beta_gamma_1;   
                beta_indicator=gamma;
                current_beta=perturbed_beta;
                linear_predictor = proposed_linear_predictor; 
                prediction = proposed_prediction;
                
            elseif (randomnumber < exp(log_R))
                
                accept(2,flag)= accept(2,flag)+1;
                X_gamma=X_gamma_1;
                V_gamma=V_gamma_1;
                m_0=m_0_gamma;
                beta_gamma=beta_gamma_1;   
                beta_indicator=gamma;
                current_beta=perturbed_beta; 
                linear_predictor = proposed_linear_predictor; 
                prediction = proposed_prediction;
                
            else
                gamma=beta_indicator;

                
                
            end
            
            if (current_beta(j)~=0)
                    indicator_vector(j)=1;
            else 
                    indicator_vector(j)=0;
            end
            
            if (iteration>burn)
                times_in_model=times_in_model+indicator_vector;
            end
           

      else

          
        end
      
       
      j=j+1;
      
      
    end

         
    i=i+1;
end
end
