// P53LogReg.sce
// Classification with logistic regression
// ----------------------------------------
[u,t,n]=file();                         // find working directory
chdir(dirname(n(2)));                   // set working directory
clear("u","t","n")                      // clear auxiliary data
exec("ScIntro.sce",-1),mode(0)          // intro to sesion                 

nL=100;                       // number of data for learning
nT=100;                       // number of data for testing
sd=1;                         // magnitude of disturbance

x1=rand(1,nL+nT,'norm');      // 1. regression vector
x2=rand(1,nL+nT,'unif')<.3;   // 2. regression vector
x=[x1;x2]';                   // full regression vector
sL=1:nL;                      // inteval for learning
sT=1:nT;                      // interval for testing
yy=2*x1-3*x2+1+sd*rand(1,nL+nT,'norm');
y=double(yy>0);               // output (true classes)
yL=y(sL);  
xL=x(sL,:);

// Learning
// py     probabilities of y=0 and y=1
// yp     prediction of y  yp=py(:,2)
// yr     point prediction of y  yr=round(yp)
// b      model parameters
[py,yp,yr,b]=lrEst([],yL,xL); // log. reg. estimation

// Testing
ct=zeros(sT);
yT=y(1,nL+sT);                // true classes
xT=[ones(nT,1) x(nL+sT,:)];   // data for testing
for t=sT
  z=xT(t,:)*b;                // regression
  yp=exp(z)/(1+exp(z));       // logistic function
  ct(t)=round(yp);            // class point estimates
end

// RESULTS
s=1:nT;
set(scf(1),'position',[600 50 600 600])
plot(s,yT(s),'bo',s,ct(s),'r.')
set(gca(),"data_bounds",[min(s)-.1 max(s)+.1 -.1 1.1])
title('Logistic regression')
legend('output','point estimate');

disp('Logistic regression:')
printf(' Wrong %d from %d\n',sum(yT~=ct),nT)
