function [ boot,stored_beta,accept,times_in_model,markers_in_model ] = human_MCMC_norm( runs, SNP, vfrac,scale)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

load Region18qnhnc2.txt
data= Region18qnhnc2;
data2=data;


data2(:,5)=[];
data2(:,4)=[];
data2(:,3)=[];
data2(:,2)=[];
data2(:,1)=[];
boot=zeros(2000,183);
beta=zeros(183,1);
beta(SNP)=vfrac;

n=1;
while n<1001
    
    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=1001;

while m<2001
          
          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=x_dim(1);
no_cases=0.5*no_patients;
%Starting Beta Values
current_beta=zeros(no_SNPs,1);

stored_beta=zeros(no_runs+1,no_SNPs);

perturbed_beta=0;

indicator_vector=-1*ones(no_SNPs,1);

perturbed_density=0;

current_density=0;

current_likelihood=0;

perturbed_likelihood=0;

% vector of indicators for whether assigned to zero or not


beta_intercept=0*ones(no_patients,1);

linear_predictor = zeros(no_patients,1);

prediction=ones(no_patients,1)./(1+exp(-(linear_predictor)));

Rmatrix=zeros(no_runs,no_SNPs);

Acceptancematrix=zeros(no_runs,no_SNPs);

mastore=zeros(no_SNPs,no_SNPs);
mastore2=zeros(no_SNPs,no_SNPs);



scale=sqrt(scale);
eta=1;
SNP_prior=1;
death_rate=0.3;
switch_rate=0.4;
switch_beta=0;
diff=0;
i=1;

norm_vec = eta*normrnd(0,1,no_runs,no_SNPs);


while i<=no_runs
    j=1;
    while j<=no_SNPs
      variablenumber=rand;
      
      if(variablenumber<0.8)
 		
 			thresh_pi=rand;
            
 	

			noise=norm_vec(i,j);
            
            if ( current_beta(j)~=0 && thresh_pi < death_rate)
                perturbed_beta = 0;
                flag=2;
                else 
                perturbed_beta=current_beta(j)+noise;
                 if (current_beta(j)==0)
                    flag=1;
                else
                    flag=3;
                end
            end
            
            if (current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
                     diff=floor((10-1+1)*rand+1);
                     random_sign=rand;
                     if random_sign<0.5
                         diff=diff*-1;
                     end
                     while ((j+diff)>no_SNPs || (j+diff)<1 || current_beta(j+diff)==1)
                         diff=floor((10-1+1)*rand+1);
                     random_sign=rand;
                     if random_sign<0.5
                         diff=diff*-1;
                     end
                     end
                     
                    
                     perturbed_beta=0;
                     switch_beta=current_beta(j+diff)+noise;
                     flag=5;
            end
            
            accept(1,flag)= accept(1,flag)+1;
            
            %change_in_log_prior = -(1/scale)*(abs(perturbed_beta)-abs(current_beta(j)));
            change_in_log_prior = -(1/(2*scale^2))*(perturbed_beta^2-current_beta(j)^2);
            if(current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
            change_in_log_prior = -(1/(2*scale^2))*(perturbed_beta^2+switch_beta^2-current_beta(j)^2);
            end
            change_in_beta = -current_beta(j) + perturbed_beta;
            
            if(current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
            change_in_beta = -current_beta(j) + perturbed_beta;
            change_in_switch=-current_beta(j+diff)+switch_beta;
            end
            proposed_linear_predictor = linear_predictor + change_in_beta*X(:,j);
            
            if(current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
            proposed_linear_predictor = linear_predictor + change_in_beta*X(:,j)+change_in_switch*X(:,j+diff);
            end
            
            
            cs = [prediction(1:no_cases); 1-prediction(no_cases+1:no_patients)];
            current_log_likelihood = sum(log(cs));
           
            %current_log_likelihood=0;
            pb=current_beta;
            pb(j)=perturbed_beta;
            
           
            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));
		
            %perturbed_log_likelihood=0;
            
            %calculate ratio=prior raio * likelihood ratio * proposal ratio
			
            if (current_beta(j)==0)
                %q_beta = -sqrt(eta) times proposed value sqaured
                log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior - 0.5*log(2*pi)-0.5*(perturbed_beta)^2+log(0.2)+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
            
            elseif (current_beta(j)~=0 && thresh_pi < death_rate)
                log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior - log(0.2)+0.5*log(2*pi)+0.5*(current_beta(j))^2-log(1/(no_SNPs-SNP_prior))+log(1-(1/(no_SNPs-SNP_prior)));
            
                %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
            Rmatrix(i,j)=log_R;
            
            randomnumber=rand;

            % MONTIOR ACCEPTANCE RATES

            
            if (log_R > 0)
                accept(2,flag)= accept(2,flag)+1;
                if (current_beta(j)~=0 && thresh_pi < 0.2)
                    indicator_vector(j)=1;
                    
                elseif (current_beta(j)==0)
                    indicator_vector(j)=0;%-1;
                   
                else
                    
                    indicator_vector(j)=1;
                end
                Acceptancematrix(i,j)=1;
                stored_beta(i+1,j) = perturbed_beta;
                if(current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
                    current_beta(j+diff)=switch_beta;
                    stored_beta(j+diff)=switch_beta;
                end
                current_beta(j) = perturbed_beta;
                linear_predictor = proposed_linear_predictor; 
                prediction = proposed_prediction;
                %mastore=mastore+indicator_vector*indicator_vector';
              mastore2(j,:)=indicator_vector'*indicator_vector(j);
              mastore2(:,j)=indicator_vector*indicator_vector(j);
              mastore=mastore+mastore2;
                
            elseif (randomnumber < exp(log_R))
                accept(2,flag)= accept(2,flag)+1;
                if (current_beta(j)~=0 && thresh_pi < 0.2)
                    indicator_vector(j)=1;
                    
                elseif (current_beta(j)==0)
                    indicator_vector(j)=-0;%1;
                    
                else
                    indicator_vector(j)=1;
                    
                end
                Acceptancematrix(i,j)=1;
                stored_beta(i+1,j) = perturbed_beta;
                if(current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
                    current_beta(j+diff)=switch_beta;
                    stored_beta(j+diff)=switch_beta;
                end
                current_beta(j) = perturbed_beta;
                linear_predictor = proposed_linear_predictor; 
                prediction = proposed_prediction;
                %mastore=mastore+indicator_vector*indicator_vector';
              mastore2(j,:)=indicator_vector'*indicator_vector(j);
              mastore2(:,j)=indicator_vector*indicator_vector(j);
              mastore=mastore+mastore2;
            else
                stored_beta(i+1,j)=current_beta(j);
                mastore=mastore+mastore2;
            end

      else
          if (j==1)
              neighbour=2;
          elseif (j==no_SNPs)
              neighbour=no_SNPs-1;
          else
              random=rand;
              if(random<0.5)
                  neighbour=j-1;
              else
                  neighbour=j+1;
              end
          end

          mu=zeros(1,2);
	 	  alpha=0.6;
	 	  sigma=[1,-alpha;-alpha,1];
          flag=4;
          accept(1,flag)= accept(1,flag)+1;
	 	 
	 	  noise=mvnrnd(mu,sigma);
	 	  perturbed_beta=current_beta;

          
          if (j>neighbour)
              perturbed_beta(j)=current_beta(j)+noise(1);
              perturbed_beta(neighbour)=current_beta(neighbour)+noise(2);
          else
              perturbed_beta(j)=current_beta(j)+noise(2);
              perturbed_beta(neighbour)=current_beta(neighbour)+noise(1);
          end

p_current=sum(indicator_vector);
          p_proposed=p_current;
          if (current_beta(j)==0)
              p_proposed=p_proposed+1;
          end
            if (current_beta(neighbour)==0)
              p_proposed=p_proposed+1;
          end

          % use next line
		  
          change_in_log_prior = -(1/scale)*(abs(perturbed_beta(j))-abs(current_beta(j)) + abs(perturbed_beta(neighbour))-abs(current_beta(neighbour)));
          change_in_beta_j = current_beta(j) - perturbed_beta(j);
          change_in_beta_neighbour = current_beta(neighbour) - perturbed_beta(neighbour);
          proposed_linear_predictor = linear_predictor + change_in_beta_j*X(:,j)+change_in_beta_neighbour*X(:,neighbour);
          
          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)); 
          %perturbed_log_likelihood=0;

          
          %calculate ratio=prior raio * likelihood ratio * proposal ratio
		  log_R = (perturbed_log_likelihood-current_log_likelihood) + change_in_log_prior+(p_proposed-p_current)*log(1/(no_SNPs-SNP_prior))-(p_proposed-p_current)*log(1-(1/(no_SNPs-SNP_prior)));
          % proposal ratio= 1
          
          Rmatrix(i,j)=log_R;
          randomnumber=rand;
          
          if (log_R > 0)
             
              accept(2,flag)= accept(2,flag)+1;
              Acceptancematrix(i,j)=1;
              stored_beta(i+1,j) = perturbed_beta(j);
              stored_beta(i+1,neighbour) = perturbed_beta(neighbour);
              current_beta(j) = perturbed_beta(j);
              current_beta(neighbour) = perturbed_beta(neighbour);
              linear_predictor = proposed_linear_predictor; 
              prediction = proposed_prediction;
              indicator_vector(j)=1;
              indicator_vector(neighbour)=1;
              %mastore=mastore+indicator_vector*indicator_vector';
              mastore2(j,:)=indicator_vector'*indicator_vector(j);
              mastore2(:,j)=indicator_vector*indicator_vector(j);
              mastore2(neighbour,:)=indicator_vector'*indicator_vector(neighbour);
              mastore2(:,neighbour)=indicator_vector*indicator_vector(neighbour);
              mastore=mastore+mastore2;
          
          elseif (randomnumber < exp(log_R))
             
              accept(2,flag)= accept(2,flag)+1;
              Acceptancematrix(i,j)=1;
              stored_beta(i+1,j) = perturbed_beta(j);
              stored_beta(i+1,neighbour) = perturbed_beta(neighbour);
              current_beta(j) = perturbed_beta(j);
              current_beta(neighbour) = perturbed_beta(neighbour);
              linear_predictor = proposed_linear_predictor; 
              prediction = proposed_prediction;
              indicator_vector(j)=1;
              indicator_vector(neighbour)=1;
              %mastore=mastore+indicator_vector*indicator_vector';
              mastore2(j,:)=indicator_vector'*indicator_vector(j);
              mastore2(:,j)=indicator_vector*indicator_vector(j);
              mastore2(neighbour,:)=indicator_vector'*indicator_vector(neighbour);
              mastore2(:,neighbour)=indicator_vector*indicator_vector(neighbour);
              mastore=mastore+mastore2;
          
          else

              stored_beta(i+1,j) = current_beta(j);
              stored_beta(i+1,neighbour) = current_beta(neighbour);
              mastore=mastore+mastore2;
          end

      end
      
        

      j=j+1;
      
      
    end

  
 	 
         
    i=i+1;
end
markers_in_model=zeros(1,(no_runs+1));
     times_in_model=zeros(1,no_SNPs);

for i=1:no_SNPs
    for j=1:(no_runs+1)
	
        if(stored_beta(j,i)~=0)
		times_in_model(i)=times_in_model(i)+1;
        end
    end
end

for i=1:(no_runs+1)
    for j=1:no_SNPs
        if(stored_beta(i,j)~=0)
		markers_in_model(i)=markers_in_model(i)+1;
        end
    end
end

end

