logistic_MCP_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;       
      %%%%%%%%%%%%%%%%%% MCP Thresholding(Zhang,et al 2010) %%%%%%%%%%%%%%%%%
      gamma=19;%20;
      if abs(S(j))<= gamma*lambda
          if S(j) > lambda
              beta(j) = (S(j) - lambda)/((1-1/gamma));
          elseif S(j) < -lambda
              beta(j) = (S(j) + lambda)/((1-1/gamma));
          elseif abs(S(j)) <= lambda
              beta(j) = 0;
          end
      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_MCP_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')