www.pudn.com > HMM1.zip > logist2.m

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