logistic_SCAD_simulation.m

clear
clc
%rng(1);

%%%%%%%%%% generate benchmark dataset %%%%%%%%%%%%%%%%%%%%%%%%
beta=zeros(1,1000);
beta(1)=5;
beta(2)=-5;
beta(3)=5;
beta(4)=-5;
beta(5)=5;
beta(6)=-5;
beta(7)=5;
beta(8)=-5;
beta(9)=5;
beta(10)=-5;

% beta(1)=1.2;
% beta(4)=1.6;
% beta(7)=0.9;
% beta(15)=0.6;
% beta(19)=0.5;
% beta(23)=-1.2;
% beta(26)=1;
% beta(30)=-0.5;
% beta(35)=1.3;
% beta(36)=0.8;

actual_beta=beta; 
train_size=100;
test_size=50;
sample_size=train_size+test_size;
intercept=0.0;
X = normrnd(0, 1, sample_size, size(beta,2)+1);
[n,p]=size(X);
cor=0.0;
for i=1:n
    for j=1:p-1
        x(i,j)=X(i,j+1)*sqrt(1-cor)+X(i,1)*sqrt(cor);
    end
end

l=intercept+(x*beta'+0.2*normrnd(0, 1, n, 1));
prob=exp(l)./(1 + exp(l));
U=rand(1,sample_size);

for i=1:sample_size
    if prob(i)>0.5
        y(i)=1;
    else
        y(i)=0;
    end
end
y=y';

x_test=x(train_size+1:sample_size,:);
x=x(1:train_size,:);
y_test=y(train_size+1:sample_size,:);
y=y(1:train_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;   %%%%% eta= intercept + X*beta;
       Pi=exp(eta)./(1+exp(eta));
       W=diag(Pi.*(1-Pi));           %%%%%%%%% W is diagonal matrix%%%%%%%%%%
       [m,n]=size(W);
       r=(W^-1)*(y-Pi);            %residual= (w^-1)*(y-pi)
    
    %%%%%%%%%%%%%%%%%%%% intercept%%%%%%%%%%%%%%%%%%%%%%
        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;       
        %%%%%%%%%%%%%%%%%%%% SCAD Thresholding (FAN & Li, 2001)%%%%%%%%%%%%%%%%
        alpha=3.7; 
        s=0;
        if S(j)>0
            s=1;
        elseif S(j)<0
            s=-1;       
        end
        if abs(S(j)) <=lambda 
            beta(j)=0;
        elseif abs(S(j)) <= 2*lambda 
            beta(j)=s*(abs(S(j))-lambda);    
        elseif abs(S(j)) <= alpha*lambda
            beta(j)=s*(abs(S(j))-alpha*lambda/(alpha-1))/(1-1/(alpha-1));
        else
            beta(j)=S(j);
        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_SCAD_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')