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