src/lsp/nonlin.lsp

   1 ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
   2 ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
   3 ;;;; You may give out copies of this software; for conditions see the file
   4 ;;;; COPYING included with this distribution.
   5
   6 (provide "nonlin")
   7
   8 (require "regress")
   9
  10 ;;;;
  11 ;;;;
  12 ;;;; Nonlinear Regression Model Prototype
  13 ;;;;
  14 ;;;;
  15
  16 (defproto nreg-model-proto
  17           '(mean-function theta-hat epsilon count-limit verbose)
  18           '()
  19           regression-model-proto)
  20
  21 (defun nreg-model (mean-function y theta
  22                                  &key
  23                                  (epsilon .0001)
  24                                  (print t)
  25                                  (count-limit 20)
  26                                  parameter-names
  27                                  response-name
  28                                  case-labels
  29                                  weights
  30                                  included
  31                                  (verbose print))
  32 "Args: (mean-function y theta &key (epsilon .0001) (count-limit 20)
  33                       (print t) parameter-names response-name case-labels
  34                       weights included (verbose print))
  35 Fits nonlinear regression model with MEAN-FUNCTION and response Y using initial
  36 parameter guess THETA. Returns model object."
  37   (let ((m (send nreg-model-proto :new)))
  38     (send m :mean-function mean-function)
  39     (send m :y y)
  40     (send m :new-initial-guess theta)
  41     (send m :epsilon epsilon)
  42     (send m :count-limit count-limit)
  43     (send m :parameter-names parameter-names)
  44     (send m :response-name response-name)
  45     (send m :case-labels case-labels)
  46     (send m :weights weights)
  47     (send m :included included)
  48     (send m :verbose verbose)
  49     (if print (send m :display))
  50     m))
  51
  52 (defmeth nreg-model-proto :save ()
  53 "Message args: ()
  54 Returns an expression that will reconstruct the regression model."
  55   `(nreg-model ',(send self :mean-function)
  56                ',(send self :y)
  57                ',(send self :coef-estimates)
  58                :epsilon ',(send self :epsilon)
  59                :count-limit ',(send self :count-limit)
  60                :predictor-names ',(send self :predictor-names)
  61                :response-name ',(send self :response-name)
  62                :case-labels ',(send self :case-labels)
  63                :weights ',(send self :weights)
  64                :included ',(send self :included)
  65                :verbose ',(send self :verbose)))
  66
  67 ;;;
  68 ;;; Computing Method
  69 ;;;
  70
  71 (defmeth nreg-model-proto :compute ()
  72 "Message args: ()
  73 Recomputes the estimates. For internal use by other messages"
  74   (let* ((y (send self :y))
  75          (weights (send self :weights))
  76          (inc (if-else (send self :included) 1 0))
  77          (w (if weights (* inc weights) inc)))
  78     (setf (slot-value 'theta-hat)
  79           (nlreg (send self :mean-function)
  80                  y
  81                  (slot-value 'theta-hat)
  82                  (send self :epsilon)
  83                  (send self :count-limit)
  84                  w
  85                  (send self :verbose)))
  86     (setf (slot-value 'x)
  87           (funcall (make-jacobian (slot-value 'mean-function)
  88                                   (length (slot-value 'theta-hat)))
  89                    (slot-value 'theta-hat)))
  90     (setf (slot-value 'intercept) nil)
  91     (call-next-method)
  92     (let ((r (send self :residuals)))
  93       (setf (slot-value 'residual-sum-of-squares)
  94             (sum (* inc r r))))))
  95
  96 ;;;
  97 ;;; Slot Accessors and Mutators
  98 ;;;
  99
 100 (defmeth nreg-model-proto :new-initial-guess (guess)
 101 "Message args: (guess)
 102 Sets a new initial uess for parmeters."
 103   (setf (slot-value 'theta-hat) guess)
 104   (send self :needs-computing t))
 105
 106 (defmeth nreg-model-proto :theta-hat ()
 107 "Message args: ()
 108 Returns current parameter estimate."
 109   (if (send self :needs-computing) (send self :compute))
 110   (coerce (slot-value 'theta-hat) 'list))
 111
 112 (defmeth nreg-model-proto :mean-function (&optional f)
 113 "Message args: (&optional f)
 114 With no argument returns the mean function as supplied to m. With an
 115 argument F sets the mean function of m to F and recomputes the
 116 estimates."
 117   (when (and f (functionp f))
 118         (setf (slot-value 'mean-function) f)
 119         (send self :needs-computing t))
 120   (slot-value 'mean-function))
 121
 122 (defmeth nreg-model-proto :epsilon (&optional eps)
 123 "Message args: (&optional eps)
 124 With no argument returns the tolerance as supplied to m. With an argument
 125 EPS sets the tolerance of m to EPS and recomputes the estimates."
 126   (when (and eps (numberp eps))
 127         (setf (slot-value 'epsilon) eps)
 128         (send self :needs-computing t))
 129    (slot-value 'epsilon))
 130
 131 (defmeth nreg-model-proto :count-limit (&optional count)
 132 "Message args: (&optional new-count)
 133 With no argument returns the iteration count limit as supplied to m. With
 134 an argument COUNT sets the limit to COUNT and recomputes the
 135 estimates."
 136   (when (and count (numberp count))
 137         (setf (slot-value 'count-limit) count)
 138         (send self :needs-computing t))
 139   (slot-value 'count-limit))
 140
 141 (defmeth nreg-model-proto :parameter-names (&optional (names nil set))
 142 "Method args: (&optional names)
 143 Sets or returns parameter names."
 144   (if set (setf (slot-value 'predictor-names) names))
 145   (let ((p-names (slot-value 'predictor-names))
 146         (p (length (slot-value 'theta-hat))))
 147     (if (not (and p-names (= p (length p-names))))
 148         (setf (slot-value 'predictor-names)
 149               (mapcar #'(lambda (a) (format nil "Parameter ~a" a))
 150                       (iseq 0 (- p 1))))))
 151   (slot-value 'predictor-names))
 152
 153 (defmeth nreg-model-proto :verbose (&optional (val nil set))
 154 "Method args: (&optional val)
 155 Sets or retrieves verbose setting. If T iteration info is printed during
 156 optimization."
 157   (if set (setf (slot-value 'verbose) val))
 158   (slot-value 'verbose))
 159
 160 ;;;
 161 ;;; Overrides for Linear Regression Methods
 162 ;;;
 163
 164 (defmeth nreg-model-proto :x ()
 165 "Message args: ()
 166 Returns the Jacobian matrix at theta-hat."
 167    (call-next-method))
 168
 169 (defmeth nreg-model-proto :intercept (&rest args)
 170 "Message args: ()
 171 Always returns nil. (For compatibility with linear regression.)"
 172   (declare (ignore args))
 173   nil)
 174
 175 (defmeth nreg-model-proto :fit-values ()
 176 "Message args: ()
 177 Returns the fitted values for the model."
 178   (coerce (funcall (send self :mean-function) (send self :theta-hat))
 179           'list))
 180
 181 (defmeth nreg-model-proto :coef-estimates (&optional guess)
 182 "Message args: (&optional guess)
 183 With no argument returns the current parameter estimate. With an
 184 argument GUESS takes it as a new initial guess and recomputes
 185 the estimates."
 186   (if guess (send self :new-initial-guess guess))
 187   (send self :theta-hat))
 188
 189 (defmeth nreg-model-proto :predictor-names () (send self :parameter-names))
 190
 191 ;;;;
 192 ;;;;
 193 ;;;; Linear Regression Coefficients
 194 ;;;;
 195 ;;;;
 196
 197 (defun regression-coefficients (x y &key (intercept T) weights)
 198 "Args: (x y &key (intercept T) weights)
 199 Returns the coefficients of the regression of the sequence Y on the columns of
 200 the matrix X."
 201   (let* ((m (if weights
 202                 (make-sweep-matrix x y weights)
 203                 (make-sweep-matrix x y)))
 204          (n (array-dimension x 1)))
 205     (coerce (compound-data-seq
 206              (if intercept
 207                  (select (car (sweep-operator m (iseq 1 n)))
 208                          (1+ n)
 209                          (iseq 0 n))
 210                  (select (car (sweep-operator m (iseq 0 n)))
 211                          (1+ n)
 212                          (iseq 1 n))))
 213             'vector)))
 214
 215 ;;;;
 216 ;;;;
 217 ;;;; Nonlinear Regression Functions
 218 ;;;;
 219 ;;;;
 220 (defun nlreg1 (f j y initial-beta epsilon count-limit weights verbose)
 221 "Args: (mean-function jacobian y initial-beta
 222         epsilon count-limit weights verbose)
 223 MEAN-FUNCTION returns the mean response vector for a given parameter vector.
 224 JACOBIAN returns the jacobian of MEAN-FUNCTION at a given parameter vector.
 225 Y is the observed response vector. Returns the estimated parameter vector
 226 obtained by a Gauss-Newton algorithm with backtracking that continues until
 227 the COUNT-LIMIT is reached or no component of the parameter vector changes
 228 by more than EPSILON."
 229   (labels ((rss (beta)                   ; residual sum of squares
 230                 (let ((res (- y (funcall f beta))))
 231                   (sum (if weights (* res res weights) (* res res)))))
 232            (next-beta (beta delta rss)   ; next beta by backtracking
 233                       (do* ((lambda 1 (/ lambda 2))
 234                             (new-rss (rss (+ beta delta))
 235                                      (rss (+ beta (* lambda delta)))))
 236                            ((or (< new-rss rss) (< lambda .0001))
 237                             (if (>= lambda .0001)
 238                                 (+ beta (* lambda delta))
 239                                 beta)))))
 240     (do* ((delbeta (regression-coefficients
 241                     (funcall j initial-beta)
 242                     (- y (funcall f initial-beta))
 243                     :intercept nil
 244                     :weights weights)
 245                    (regression-coefficients
 246                     (funcall j beta)
 247                     (- y (funcall f beta))
 248                     :intercept nil
 249                     :weights weights))
 250           (beta initial-beta (next-beta beta delbeta rss))
 251           (rss (rss beta) (rss beta))
 252           (count 0 (1+ count)))
 253          ((or (> count count-limit) (> epsilon (max (abs delbeta))))
 254           (if (and verbose (> count count-limit))
 255               (format t "Iteration limit exceeded.~%"))
 256           beta)
 257          (if verbose (format t "Residual sum of squares: ~,6g~%" rss)))))
 258
 259 (defun make-jacobian (f n)
 260 "Args: (f n)
 261 F is a function of an N-vector. Returns a function that approximates the
 262 jacobian function iof F by a symmetric difference."
 263   (let* ((h .0001)
 264          (del (* h (column-list (identity-matrix n)))))
 265     #'(lambda (b)
 266       (let ((b+ (mapcar #'(lambda (x) (+ b x)) del))
 267             (b- (mapcar #'(lambda (x) (- b x)) del)))
 268         (apply #'bind-columns (/ (- (mapcar f b+) (mapcar f b-)) (* 2 h)))))))
 269
 270 (defun nlreg (f y guess &optional
 271                 (epsilon .0001) (count-limit 20) weights verbose)
 272 "Args: (mean-function y guess &optional
 273         (epsilon .0001) (count-limit 20) weights verbose)
 274 MEAN-FUNCTION returns the mean response vector for a given parameter vector.
 275 Y is the observed response vector. Returns the estimated parameter vector
 276 obtained by a Gauss-Newton algorithm that continues until the ITERATION-LIMIT
 277 is reached or no component of the parameter vector changes by more than
 278 EPSILON. The jacobian of MEAN-FUNCTION is approximated by a symmetric difference."
 279   (nlreg1 f (make-jacobian f (length guess)) y guess
 280           epsilon count-limit weights verbose))