function [qinfo,stored_beta,accept,qvarmonitor,log_R_track, times_in_model,markers_in_model,frequentist_log_likelihood ] = mouse_MCMC_function_new( runs, SNP, vfrac,scale )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

% Load data
  load crlCFW-ALP.txt %199p 71 snp
  data= crlCFW_ALP;
  data2=data;
  

 
  %Remove unneeded columns
    data2(:,5)=[];
    data2(:,4)=[];
    data2(:,3)=[];
    data2(:,2)=[];
    data2(:,1)=[];
    
  % variables for simulation
    % Genotype matrix
      X=data2;
%       for (ci=1:199)
%           for(cj=1:71)
%               if(data2(ci,cj)==0)
%                   X(ci,cj)=2;
%               elseif(data2(ci,cj)==2)
%                   X(ci,cj)=0;
%               end
%           end
%       end
    % SNP for QTL to be simed around
      causal_snp=SNP;
    % Variance explained
      v_fraction=vfrac;
      
  [qinfo,qtlsim_p,qtlsim_a] = qtlsim(X,causal_snp,v_fraction);
    
  % variables for MCMC
    % Number of runs
      no_runs=runs;
    % Number of patients and SNPs
      x_dim=size(X);
      no_SNPs=x_dim(2);
      no_patients=x_dim(1);
    % 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
     stored_beta=zeros(no_SNPs*no_runs+1,no_SNPs);
     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);
     qvarmonitor=zeros(no_SNPs*no_runs+1,1);
     gamma=zeros(1,no_SNPs);
     beta_indicator=zeros(1,no_SNPs);
     X_gamma=0;
     V_gamma=0;
     beta_gamma=[];
     log_R_track=zeros(no_runs,no_SNPs);
     qvar=0;
     
     
 
% Variables for prior and noise updates
     prior_scale=sqrt(scale);%(2*qtlsim_p*(1-qtlsim_p)/v_fraction)^0.5;
     
     

% Model
     i=1;
     iteration=1;
     
while i<=no_runs
    j=1;
    while j<=no_SNPs
        iteration=iteration+1;
    %Variance calculation
    sgamma=sum(gamma);
   if sgamma==0
       
       a_1=a_0+no_patients/2;
       V_1=0;
       m_1=0;
       b_1=b_0;

       tau=random('gam', a_1,(1/(b_1)));
       qvar=tau^-1;
       qvarmonitor(iteration)=qvar;
       
   else
       
   
    
            a_1=a_0+no_patients/2;
            V_1=(X_gamma'*X_gamma+V_gamma^-1)^-1;
            m_1=V_1*(X_gamma'*qinfo+V_gamma^-1*m_0);
            b_1=b_0+0.5*(qinfo'*qinfo+m_0'*V_gamma^-1*m_0-m_1'*V_1^-1*m_1);
            tau=random('gam', a_1,(1/(b_1)));
            qvar=tau^-1;
           
            qvarmonitor(iteration)=qvar;
            beta_gamma=mvnrnd(m_1,(prior_scale*qvar*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;
%                 elseif (current_beta(j)~=0 && thresh_pi > death_rate && thresh_pi<(death_rate+switch_rate))
%                      diff=floor((no_SNPs-1+1)*rand+1);
%                      random_sign=rand;
%                      if random_sign<0.5
%                          diff=diff*-1;
%                      end
%                      while j==diff || (diff)<1 || gamma(diff)==1
%                          diff=floor((no_SNPs-1+1)*rand+1);
%                      random_sign=rand;
%                      if random_sign<0.5
%                          diff=diff*-1;
%                      end
%                      end
%                      
%                     
%                      gamma(j)=0;
%                      gamma(diff)=1;
%                      flag=5;
%                      
                         
                     
                         
            else 
                gamma(j) = 1;
                 if (current_beta(j)==0)
                    flag=1;
                else
                    flag=3;
                end
            end
            
            X_gamma_1=[];
            for col=1:no_SNPs
                if gamma(col)==1
                    X_gamma_1=[X_gamma_1 X(:,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=0;
                beta_gamma_1=0;
                
            else
                
           
            V_gamma_1=eye(xgamma_dim(2));
            m_0_gamma=zeros(xgamma_dim(2),1);
            V_1_gamma=(X_gamma_1'*X_gamma_1+V_gamma_1^-1)^-1;
            m_1_gamma=V_1_gamma*(X_gamma_1'*qinfo+V_gamma_1^-1*m_0_gamma);
            
            beta_gamma_1=mvnrnd(m_1_gamma,prior_scale*qvar*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/prior_scale)*sum((abs(perturbed_beta)-abs(current_beta)));
            linear_predictor = X*current_beta;
            
            proposed_linear_predictor = X*perturbed_beta;
            
            if likelihood_flag==1
            current_log_likelihood=-0.5*no_patients*log(qvar) -0.5*sum((qinfo-linear_predictor).^2/(qvar));
            perturbed_log_likelihood=-0.5*no_patients*log(qvar) -0.5*sum((qinfo-proposed_linear_predictor).^2/(qvar));
            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,prior_scale*qvar*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,prior_scale*qvar*V_1_gamma))+log(death_rate)+log(mvnpdf(beta_gamma',m_1,prior_scale*qvar*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,prior_scale*qvar*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,prior_scale*qvar*V_1_gamma))+log(mvnpdf(beta_gamma',m_1,prior_scale*qvar*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;
            log_R_track(i,j)=log_R;

            % 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;
                stored_beta(iteration,:)=perturbed_beta;
                current_beta=perturbed_beta;
                linear_predictor = proposed_linear_predictor; 
                
            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;
                stored_beta(iteration,:)=perturbed_beta;
                current_beta=perturbed_beta; 
                linear_predictor = proposed_linear_predictor; 
                
            else
                gamma=beta_indicator;
                stored_beta(iteration,:)=current_beta;
                
                
            end
            
            
           lpa=1;
           lpb=1;
           prior_scale_in=random('gam',(lpa+sum(beta_indicator)),(1/(lpb+sum(abs(current_beta)))));
           prior_scale=prior_scale_in^-1;
%            tau=random('gam', a_1,(1/(b_1)));
%             qvar=tau^-1;

      else

          
        end
      
       
      j=j+1;
      
      
    end

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

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

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

for j=1:71
    tlinear_predictor =zeros(no_patients,1);
      
    for i=1:no_patients
    if data2(i,j)==0
        
        tlinear_predictor(i)=-qtlsim_a;
    end
    
    if data2(i,j)==1
     
        tlinear_predictor(i)=0;
    end
    
    if data2(i,j)==2
        
        tlinear_predictor(i)=qtlsim_a;
    end
    end
    
    beta_mle= (sum(qinfo)*sum(tlinear_predictor)-no_patients*sum(qinfo.*tlinear_predictor))/(sum(tlinear_predictor)*sum(tlinear_predictor)-no_patients*sum(tlinear_predictor.*tlinear_predictor));
    mu_mle= (1/no_patients)*(sum(qinfo)-beta_mle*sum(tlinear_predictor));
    sigma_mle=(1/no_patients)*(sum((qinfo-beta_mle*tlinear_predictor-mu_mle*(ones(no_patients,1))).^2));

    %slp(j)=sum(linear_predictor);
    frequentist_log_likelihood(j)=-0.5*no_patients*log(2*pi*sigma_mle) -0.5*sum((qinfo-mu_mle*(ones(no_patients,1))-beta_mle*tlinear_predictor).^2/(sigma_mle));
    
 
end

end