Mercurial > hg > octave-nkf
diff scripts/statistics/models/logistic_regression.m @ 3426:f8dde1807dee
[project @ 2000-01-13 08:40:00 by jwe]
author | jwe |
---|---|
date | Thu, 13 Jan 2000 08:40:53 +0000 |
parents | 041ea33fbbf4 |
children | d8b731d3f7a3 |
line wrap: on
line diff
--- a/scripts/statistics/models/logistic_regression.m +++ b/scripts/statistics/models/logistic_regression.m @@ -1,15 +1,15 @@ ## Copyright (C) 1995, 1996, 1997 Kurt Hornik -## +## ## 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, 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. -## +## General Public License for more details. +## ## You should have received a copy of the GNU General Public License ## along with this file. If not, write to the Free Software Foundation, ## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. @@ -50,11 +50,11 @@ ## ## `p' holds estimates for the conditional distribution of Y given x. -## Original for MATLAB written by Gordon K Smyth <gks@maths.uq.oz.au>, +## Original for MATLAB written by Gordon K Smyth <gks@maths.uq.oz.au>, ## U of Queensland, Australia, on Nov 19, 1990. Last revision Aug 3, ## 1992. -## Author: Gordon K Smyth <gks@maths.uq.oz.au>, +## Author: Gordon K Smyth <gks@maths.uq.oz.au>, ## Adapted-By: KH <Kurt.Hornik@ci.tuwien.ac.at> ## Description: Ordinal logistic regression @@ -63,38 +63,38 @@ function [theta, beta, dev, dl, d2l, p] ... = logistic_regression (y, x, print, theta, beta) - + ## check input - y = round (vec (y)); - [my, ny] = size (y); + y = round (vec (y)); + [my, ny] = size (y); if (nargin < 2) - x = zeros (my, 0); + x = zeros (my, 0); endif; [mx, nx] = size (x); if (mx != my) error ("x and y must have the same number of observations"); endif - + ## initial calculations x = -x; tol = 1e-6; incr = 10; decr = 2; ymin = min (y); ymax = max (y); yrange = ymax - ymin; z = (y * ones (1, yrange)) == ((y * 0 + 1) * (ymin : (ymax - 1))); z1 = (y * ones (1, yrange)) == ((y * 0 + 1) * ((ymin + 1) : ymax)); - z = z(:, any (z)); - z1 = z1 (:, any(z1)); + z = z(:, any (z)); + z1 = z1 (:, any(z1)); [mz, nz] = size (z); - + ## starting values if (nargin < 3) - print = 0; + print = 0; endif; - if (nargin < 4) - beta = zeros (nx, 1); + if (nargin < 4) + beta = zeros (nx, 1); endif; - if (nargin < 5) - g = cumsum (sum (z))' ./ my; - theta = log (g ./ (1 - g)); + if (nargin < 5) + g = cumsum (sum (z))' ./ my; + theta = log (g ./ (1 - g)); endif; tb = [theta; beta]; @@ -102,7 +102,7 @@ [g, g1, p, dev] = logistic_regression_likelihood (y, x, tb, z, z1); [dl, d2l] = logistic_regression_derivatives (x, z, z1, g, g1, p); epsilon = std (vec (d2l)) / 1000; - + ## maximize likelihood using Levenberg modified Newton's method iter = 0; while (abs (dl' * (d2l \ dl) / length (dl)) > tol) @@ -115,12 +115,12 @@ epsilon = epsilon / decr; else while ((dev - devold) / (dl' * (tb - tbold)) > 0) - epsilon = epsilon * incr; + epsilon = epsilon * incr; if (epsilon > 1e+15) - error ("epsilon too large"); + error ("epsilon too large"); endif - tb = tbold - (d2l - epsilon * eye (size (d2l))) \ dl; - [g, g1, p, dev] = logistic_regression_likelihood (y, x, tb, z, z1); + tb = tbold - (d2l - epsilon * eye (size (d2l))) \ dl; + [g, g1, p, dev] = logistic_regression_likelihood (y, x, tb, z, z1); disp ("epsilon"); disp (epsilon); endwhile endif @@ -141,19 +141,19 @@ if (print >= 1) printf ("\n"); printf ("Logistic Regression Results:\n"); - printf ("\n"); + printf ("\n"); printf ("Number of Iterations: %d\n", iter); printf ("Deviance: %f\n", dev); printf ("Parameter Estimates:\n"); printf (" Theta S.E.\n"); - se = sqrt (diag (inv (-d2l))); + se = sqrt (diag (inv (-d2l))); for i = 1 : nz printf (" %8.4f %8.4f\n", tb (i), se (i)); endfor if (nx > 0) printf (" Beta S.E.\n"); for i = (nz + 1) : (nz + nx) - printf (" %8.4f %8.4f\n", tb (i), se (i)); + printf (" %8.4f %8.4f\n", tb (i), se (i)); endfor endif endif @@ -166,5 +166,5 @@ endif gamma = diff ([(y * 0), (exp (e) ./ (1 + exp (e))), (y * 0 + 1)]')'; endif - + endfunction