wolffd@0: function [beta,p,lli] = logist2(y,x,w)
wolffd@0: % [beta,p,lli] = logist2(y,x) 
wolffd@0: %
wolffd@0: % 2-class logistic regression.  
wolffd@0: %
wolffd@0: % INPUT
wolffd@0: % 	y 	Nx1 colum vector of 0|1 class assignments
wolffd@0: % 	x 	NxK matrix of input vectors as rows
wolffd@0: % 	[w]	Nx1 vector of sample weights 
wolffd@0: %
wolffd@0: % OUTPUT
wolffd@0: % 	beta 	Kx1 column vector of model coefficients
wolffd@0: % 	p 	Nx1 column vector of fitted class 1 posteriors
wolffd@0: % 	lli 	log likelihood
wolffd@0: %
wolffd@0: % Class 1 posterior is 1 / (1 + exp(-x*beta))
wolffd@0: %
wolffd@0: % David Martin <dmartin@eecs.berkeley.edu> 
wolffd@0: % April 16, 2002
wolffd@0: 
wolffd@0: % Copyright (C) 2002 David R. Martin <dmartin@eecs.berkeley.edu>
wolffd@0: %
wolffd@0: % This program is free software; you can redistribute it and/or
wolffd@0: % modify it under the terms of the GNU General Public License as
wolffd@0: % published by the Free Software Foundation; either version 2 of the
wolffd@0: % License, or (at your option) any later version.
wolffd@0: % 
wolffd@0: 
wolffd@0: % This program is distributed in the hope that it will be useful, but
wolffd@0: % WITHOUT ANY WARRANTY; without even the implied warranty of
wolffd@0: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
wolffd@0: % General Public License for more details.
wolffd@0: % 
wolffd@0: % You should have received a copy of the GNU General Public License
wolffd@0: % along with this program; if not, write to the Free Software
wolffd@0: % Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
wolffd@0: % 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html.
wolffd@0: 
wolffd@0: error(nargchk(2,3,nargin));
wolffd@0: 
wolffd@0: % check inputs
wolffd@0: if size(y,2) ~= 1,
wolffd@0:   error('Input y not a column vector.');
wolffd@0: end
wolffd@0: if size(y,1) ~= size(x,1), 
wolffd@0:   error('Input x,y sizes mismatched.'); 
wolffd@0: end
wolffd@0: 
wolffd@0: % get sizes
wolffd@0: [N,k] = size(x);
wolffd@0: 
wolffd@0: % if sample weights weren't specified, set them to 1
wolffd@0: if nargin < 3, 
wolffd@0:   w = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: % normalize sample weights so max is 1
wolffd@0: w = w / max(w);
wolffd@0: 
wolffd@0: % initial guess for beta: all zeros
wolffd@0: beta = zeros(k,1);
wolffd@0: 
wolffd@0: % Newton-Raphson via IRLS,
wolffd@0: % taken from Hastie/Tibshirani/Friedman Section 4.4.
wolffd@0: iter = 0;
wolffd@0: lli = 0;
wolffd@0: while 1==1,
wolffd@0:   iter = iter + 1;
wolffd@0:   
wolffd@0:   % fitted probabilities
wolffd@0:   p = 1 ./ (1 + exp(-x*beta));	
wolffd@0:   
wolffd@0:   % log likelihood
wolffd@0:   lli_prev = lli;
wolffd@0:   lli = sum( w .* (y.*log(p+eps) + (1-y).*log(1-p+eps)) );
wolffd@0: 
wolffd@0:   % least-squares weights
wolffd@0:   wt = w .* p .* (1-p);		
wolffd@0: 
wolffd@0:   % derivatives of likelihood w.r.t. beta
wolffd@0:   deriv = x'*(w.*(y-p));
wolffd@0: 
wolffd@0:   % Hessian of likelihood w.r.t. beta
wolffd@0:   % hessian = x'Wx, where W=diag(w)
wolffd@0:   % Do it this way to be memory efficient and fast.
wolffd@0:   hess = zeros(k,k);
wolffd@0:   for i = 1:k,
wolffd@0:     wxi = wt .* x(:,i);
wolffd@0:     for j = i:k,
wolffd@0:       hij = wxi' * x(:,j);
wolffd@0:       hess(i,j) = -hij;
wolffd@0:       hess(j,i) = -hij;
wolffd@0:     end
wolffd@0:   end
wolffd@0: 
wolffd@0:   % make sure Hessian is well conditioned
wolffd@0:   if (rcond(hess) < eps), 
wolffd@0:     error(['Stopped at iteration ' num2str(iter) ...
wolffd@0:            ' because Hessian is poorly conditioned.']);
wolffd@0:     break; 
wolffd@0:   end;
wolffd@0: 
wolffd@0:   % Newton-Raphson update step
wolffd@0:   step = hess\deriv;
wolffd@0:   beta = beta - step;
wolffd@0: 
wolffd@0:   % termination criterion based on derivatives
wolffd@0:   tol = 1e-6;
wolffd@0:   if abs(deriv'*step/k) < tol, break; end;
wolffd@0: 
wolffd@0:   % termination criterion based on log likelihood
wolffd@0: %   tol = 1e-4;
wolffd@0: %   if abs((lli-lli_prev)/(lli+lli_prev)) < 0.5*tol, break; end;
wolffd@0: end;
wolffd@0: