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 <dmartin@eecs.berkeley.edu>
> April 16, 2002

> Copyright (C) 2002 David R. Martin <dmartin@eecs.berkeley.edu>
>
> 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;