logistic_Lhalf_gene_selection.m

clear
clc

load x
load y

%%%% input training dataset and testing dataset %%%%
%x = x;  
%y = y;
%test_data = x;
%test_label = y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
train_size=size(x,1);
test_size=size(test_data,1);
sample_size=train_size+test_size;
col=size(x,2);
row=size(x,1);

temp=sum(y)/row;
beta_zero=log(temp/(1-temp));    %intercept
beta=zeros(col,1);

%%%%%%%%%%%%%% compute lambda on the log scale %%%%%%%%%%%%%%%%%%%%%
eta = beta_zero + x * beta;
Pi=exp(eta)./(1+exp(eta));
W=diag(Pi.*(1-Pi));           %%%%%%%%% W is diagonal matrix%%%%%%%%%%
r=(W^-1)*(y-Pi);              % residual
S=(x'*W*r)/row;
lambda_max=(4/3*(max(S)))^(1.5);
lambda_min = lambda_max*0.001;
m = 10;

for i=1:m
    Lambda(i) = lambda_max*(lambda_min/lambda_max)^(i/m);
end
for i=1:m    
   lambda=Lambda(i);
   iter=0;
   maxiter=1000;
   beta_path(:,i)=beta(:,1);
   
   while iter<maxiter
    
    beta_temp=beta;
    beta_zero_temp=beta_zero;
    
    eta = beta_zero_temp + x * beta_temp;   
    Pi=exp(eta)./(1+exp(eta));
    W=diag(Pi.*(1-Pi));         
    [m,n]=size(W);
    r=(W^-1)*(y-Pi);      
    
    %%%%%%%%%%%%%%%%%%%% initializing beta_0 and beta_k %%%%%%%%%%%%%%%%%%%%%%
    beta_zero=sum(W*r)/sum(sum(W))+beta_zero_temp;    
    r=r-(beta_zero-beta_zero_temp);                
    
    for j=1:col
        v=x(:,j)'*W*x(:,j)/row;
        S=(x(:,j)'*W*r)/row+beta_temp*v;       
        
    %%%%%%%%%%%%%%%%%% Half Thresholding(Xu,et al 2010) %%%%%%%%%%%%%%%%%
        if abs(S(j)) > ((3/4)*(lambda^(2/3)))
           phi = acos(lambda/8*((abs(S(j))/3)^(-1.5)));
           beta(j) = real(2/3*S(j)*(1 + cos(2/3*(pi - phi))))/v;
        else
           beta(j) = 0;
        end    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%--update r---%%%%%%%%%%%%%%%%%%%%%%%%%%
        r=r-x(:,j)*(beta(j)-beta_temp(j));        
    end    
    if norm(beta_temp - beta) < (1E-5) 
        break;
    end    
    iter=iter+1
   end    
end
[opt,Mse]=CV_logistic(x,y,Lambda,beta_path);

beta=beta_path(:,opt);
l= intercept + x_test * beta;
prob=exp(l)./(1 + exp(l));
for i=1:test_size
    if prob(i)>0.5
        test_y(i)=1;
    else
        test_y(i)=0;
    end
end
error=test_y'-y_test;
count=find(error~=0)
fail=length(count)

beta_non_zero=find(beta~=0);

plot(beta_path','linewidth',1.5)
ax = axis;
line([opt opt], [ax(3) ax(4)], 'Color', 'b', 'LineStyle', '-.');
xlabel('Steps')
ylabel('Coefficeints')

figure;
hold on
plot(Mse,'linewidth',1.5);
ax = axis;
line([opt opt], [ax(3) ax(4)], 'Color', 'b', 'LineStyle', '-.');
xlabel('Steps')
ylabel('Misclassification Error')