%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   M. Welling %%%%%%% March 2000 %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Logistic Regression %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
colordef black;

N_0 = input('Number of points you may want to select for 0 ');
N_1 = input('Number of points you may want to select for 1 ');

x_0=zeros(1,N_0);
y_0=zeros(1,N_0);

figure(1);clf;
axis([-10 10 -0.2 1.2]);
hold on; 

%%%%%%%% acquire points %%%%%%%%%%%%%%%%%%%
for i=1:N_0,

    fprintf(1,'%d - Select a point labelled 0 with the mouse:\n ',i);
    [x,dm] = ginput2(1);
    x_0(i) = x;
    plot(x_0(i),y_0(i),'r+');
    
end;


x_1=zeros(1,N_1);
y_1=ones(1,N_1);

for i=1:N_1,

    fprintf(1,'%d - Select a point labelled 1 with the mouse: \n',i);
    [x,dm] = ginput2(1);
    x_1(i) = x;
    plot(x_1(i),y_1(i),'r+');
    
end;


X = [x_0 x_1;ones(1,N_0+N_1)]';
Y = [y_0 y_1]';

eta = 0.03;
th = [0 0]';
th2 = [0 0]';

t = -10:0.01:10;
T = [t;ones(1,length(t))]';

sg_t = 1./(1+exp(-T*th));
sg_t2 = 1./(1+exp(-T*th2));

figure(2);clf;
for i=1:200
    
    %% gradient =descend update %%%%%%%
    
    sg = 1./(1+exp(-X*th));
    th = th + eta * X'*(Y - sg);
    
    tht(i,:) = th';
    
    %% IRLS update %%%%%%%%%%%%%%%%%%%%
    
    sg2 = 1./(1+exp(-X*th2));
    W = diag(sg2 .* (1-sg2));
    z = X*th2 + inv(W)*(Y-sg2);
    
    th2 = inv(X'*W*X) *X'*W*z; 
    
    tht2(i,:) = th2';

        
 %%%%%%%%%% plotting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 
    sg_t = 1./(1+exp(-T*th));
    sg_t2 = 1./(1+exp(-T*th2));
    
    figure(1);clf
    plot(x_0,y_0,'r+');hold on;
    plot(x_1,y_1,'r+');
    plot(t,sg_t,'g-','linewidth',2);
    plot(t,sg_t2,'c-','linewidth',2);
    axis([-10 10 -0.2 1.2]);  
    drawnow;
    
    figure(2);clf;
    plot(tht(:,1),'g-','linewidth',2);hold on;
    plot(tht(:,2),'g--','linewidth',2);
    plot(tht2(:,1),'c-','linewidth',2);
    plot(tht2(:,2),'c--','linewidth',2);
    drawnow;
    
    if i <=20
    pause
    end

end
%plot([-th(2)/th(1) -th(2)/th(1)],[-0.2 1.2],'m-','linewidth',2);