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