crude hack to allow Carlos and my home directories. Feel free to add more.
[CommonLispStat.git] / nonlin.lsp
blobcf34ed8bd70657b81be517e0bd7f0b8c3ae8d520
1 ;;; -*- mode: lisp -*-
2 ;;;
3 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
4 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
5 ;;; Since 1991, ANSI was finally finished. Modified to match ANSI
6 ;;; Common Lisp.
8 ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
9 ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
10 ;;;; You may give out copies of this software; for conditions see the file
11 ;;;; COPYING included with this distribution.
13 (defpackage :lisp-stat-regression-nonlin
14 (:use :common-lisp
15 :lisp-stat-object-system
16 :lisp-stat-math
17 :lisp-stat-basics
18 :lisp-stat-compound-data
19 :lisp-stat-matrix
20 :lisp-stat-linalg
21 :lisp-stat-regression-linear)
22 (:shadowing-import-from :lisp-stat-object-system
23 slot-value call-method call-next-method)
24 (:shadowing-import-from :lisp-stat-math
25 expt + - * / ** mod rem abs 1+ 1- log exp sqrt sin cos tan
26 asin acos atan sinh cosh tanh asinh acosh atanh float random
27 truncate floor ceiling round minusp zerop plusp evenp oddp
28 < <= = /= >= > ;; complex
29 conjugate realpart imagpart phase
30 min max logand logior logxor lognot ffloor fceiling
31 ftruncate fround signum cis)
33 (:export nreg-model nreg-model-proto mean-function theta-hat epsilon
34 count-limit verbose))
36 (in-package :lisp-stat-regression-nonlin)
38 ;;;
39 ;;; Nonlinear Regression Model Prototype
40 ;;;
42 (defproto nreg-model-proto
43 '(mean-function theta-hat epsilon count-limit verbose)
44 '()
45 regression-model-proto)
47 (defun nreg-model (mean-function y theta
48 &key
49 (epsilon .0001)
50 (print t)
51 (count-limit 20)
52 parameter-names
53 response-name
54 case-labels
55 weights
56 included
57 (verbose print))
58 "Args: (mean-function y theta &key (epsilon .0001) (count-limit 20)
59 (print t) parameter-names response-name case-labels
60 weights included (vetbose print))
61 Fits nonlinear regression model with MEAN-FUNCTION and response Y using initial
62 parameter guess THETA. Returns model object."
63 (let ((m (send nreg-model-proto :new)))
64 (send m :mean-function mean-function)
65 (send m :y y)
66 (send m :new-initial-guess theta)
67 (send m :epsilon epsilon)
68 (send m :count-limit count-limit)
69 (send m :parameter-names parameter-names)
70 (send m :response-name response-name)
71 (send m :case-labels case-labels)
72 (send m :weights weights)
73 (send m :included included)
74 (send m :verbose verbose)
75 (if print (send m :display))
76 m))
78 (defmeth nreg-model-proto :save ()
79 "Message args: ()
80 Returns an expression that will reconstruct the regression model."
81 `(nreg-model ',(send self :mean-function)
82 ',(send self :y)
83 ',(send self :coef-estimates)
84 :epsilon ',(send self :epsilon)
85 :count-limit ',(send self :count-limit)
86 :predictor-names ',(send self :predictor-names)
87 :response-name ',(send self :response-name)
88 :case-labels ',(send self :case-labels)
89 :weights ',(send self :weights)
90 :included ',(send self :included)
91 :verbose ',(send self :verbose)))
93 ;;;
94 ;;; Computing Method
95 ;;;
97 (defmeth nreg-model-proto :compute ()
98 "Message args: ()
99 Recomputes the estimates. For internal use by other messages"
100 (let* ((y (send self :y))
101 (weights (send self :weights))
102 (inc (if-else (send self :included) 1 0))
103 (w (if weights (* inc weights) inc)))
104 (setf (slot-value 'theta-hat)
105 (nlreg (send self :mean-function)
107 (slot-value 'theta-hat)
108 (send self :epsilon)
109 (send self :count-limit)
111 (send self :verbose)))
112 (setf (slot-value 'x)
113 (funcall (make-jacobian (slot-value 'mean-function)
114 (length (slot-value 'theta-hat)))
115 (slot-value 'theta-hat)))
116 (setf (slot-value 'intercept) nil)
117 (call-next-method)
118 (let ((r (send self :residuals)))
119 (setf (slot-value 'residual-sum-of-squares)
120 (sum (* inc r r))))))
123 ;;; Slot Accessors and Mutators
126 (defmeth nreg-model-proto :new-initial-guess (guess)
127 "Message args: (guess)
128 Sets a new initial uess for parmeters."
129 (setf (slot-value 'theta-hat) guess)
130 (send self :needs-computing t))
132 (defmeth nreg-model-proto :theta-hat ()
133 "Message args: ()
134 Returns current parameter estimate."
135 (if (send self :needs-computing) (send self :compute))
136 (coerce (slot-value 'theta-hat) 'list))
138 (defmeth nreg-model-proto :mean-function (&optional f)
139 "Message args: (&optional f)
140 With no argument returns the mean function as supplied to m. With an
141 argument F sets the mean function of m to F and recomputes the
142 estimates."
143 (when (and f (functionp f))
144 (setf (slot-value 'mean-function) f)
145 (send self :needs-computing t))
146 (slot-value 'mean-function))
148 (defmeth nreg-model-proto :epsilon (&optional eps)
149 "Message args: (&optional eps)
150 With no argument returns the tolerance as supplied to m. With an argument
151 EPS sets the tolerance of m to EPS and recomputes the estimates."
152 (when (and eps (numberp eps))
153 (setf (slot-value 'epsilon) eps)
154 (send self :needs-computing t))
155 (slot-value 'epsilon))
157 (defmeth nreg-model-proto :count-limit (&optional count)
158 "Message args: (&optional new-count)
159 With no argument returns the iteration count limit as supplied to m. With
160 an argument COUNT sets the limit to COUNT and recomputes the
161 estimates."
162 (when (and count (numberp count))
163 (setf (slot-value 'count-limit) count)
164 (send self :needs-computing t))
165 (slot-value 'count-limit))
167 (defmeth nreg-model-proto :parameter-names (&optional (names nil set))
168 "Method args: (&optional names)
169 Sets or returns parameter names."
170 (if set (setf (slot-value 'predictor-names) names))
171 (let ((p-names (slot-value 'predictor-names))
172 (p (length (slot-value 'theta-hat))))
173 (if (not (and p-names (= p (length p-names))))
174 (setf (slot-value 'predictor-names)
175 (mapcar #'(lambda (a) (format nil "Parameter ~a" a))
176 (iseq 0 (- p 1))))))
177 (slot-value 'predictor-names))
179 (defmeth nreg-model-proto :verbose (&optional (val nil set))
180 "Method args: (&optional val)
181 Sets or retrieves verbose setting. If T iteration info is printed during
182 optimization."
183 (if set (setf (slot-value 'verbose) val))
184 (slot-value 'verbose))
187 ;;; Overrides for Linear Regression Methods
190 (defmeth nreg-model-proto :x ()
191 "Message args: ()
192 Returns the Jacobian matrix at theta-hat."
193 (call-next-method))
195 (defmeth nreg-model-proto :intercept (&rest args)
196 "Message args: ()
197 Always returns nil. (For compatibility with linear regression.)"
198 nil)
200 (defmeth nreg-model-proto :fit-values ()
201 "Message args: ()
202 Returns the fitted values for the model."
203 (coerce (funcall (send self :mean-function) (send self :theta-hat))
204 'list))
206 (defmeth nreg-model-proto :coef-estimates (&optional guess)
207 "Message args: (&optional guess)
208 With no argument returns the current parameter estimate. With an
209 argument GUESS takes it as a new initial guess and recomputes
210 the estimates."
211 (if guess (send self :new-initial-guess guess))
212 (send self :theta-hat))
214 (defmeth nreg-model-proto :predictor-names () (send self :parameter-names))
216 ;;;;
217 ;;;;
218 ;;;; Linear Regression Coefficients
219 ;;;;
220 ;;;;
222 (defun regression-coefficients (x y &key (intercept T) weights)
223 "Args: (x y &key (intercept T) weights)
224 Returns the coefficients of the regression of the sequence Y on the columns of
225 the matrix X."
226 (let* ((m (if weights
227 (make-sweep-matrix x y weights)
228 (make-sweep-matrix x y)))
229 (n (array-dimension x 1)))
230 (coerce (compound-data-seq
231 (if intercept
232 (select (car (sweep-operator m (iseq 1 n)))
233 (1+ n)
234 (iseq 0 n))
235 (select (car (sweep-operator m (iseq 0 n)))
236 (1+ n)
237 (iseq 1 n))))
238 'vector)))
240 ;;;;
241 ;;;;
242 ;;;; Nonlinear Regression Functions
243 ;;;;
244 ;;;;
245 (defun nlreg1 (f j y initial-beta epsilon count-limit weights verbose)
246 "Args: (mean-function jacobian y initial-beta
247 epsilon count-limit weights verbose)
248 MEAN-FUNCTION returns the mean response vector for a given parameter vector.
249 JACOBIAN returns the jacobian of MEAN-FUNCTION at a given parameter vector.
250 Y is the observed response vector. Returns the estimated parameter vector
251 obtained by a Gauss-Newton algorithm with backtracking that continues until
252 the COUNT-LIMIT is reached or no component of the parameter vector changes
253 by more than EPSILON."
254 (labels ((rss (beta) ; residual sum of squares
255 (let ((res (- y (funcall f beta))))
256 (sum (if weights (* res res weights) (* res res)))))
257 (next-beta (beta delta rss) ; next beta by backtracking
258 (do* ((lambda 1 (/ lambda 2))
259 (new-rss (rss (+ beta delta))
260 (rss (+ beta (* lambda delta)))))
261 ((or (< new-rss rss) (< lambda .0001))
262 (if (>= lambda .0001)
263 (+ beta (* lambda delta))
264 beta)))))
265 (do* ((delbeta (regression-coefficients
266 (funcall j initial-beta)
267 (- y (funcall f initial-beta))
268 :intercept nil
269 :weights weights)
270 (regression-coefficients
271 (funcall j beta)
272 (- y (funcall f beta))
273 :intercept nil
274 :weights weights))
275 (beta initial-beta (next-beta beta delbeta rss))
276 (rss (rss beta) (rss beta))
277 (count 0 (1+ count)))
278 ((or (> count count-limit) (> epsilon (max (abs delbeta))))
279 (if (and verbose (> count count-limit))
280 (format t "Iteration limit exceeded.~%"))
281 beta)
282 (if verbose (format t "Residual sum of squares: ~10g~%" rss)))))
284 (defun make-jacobian (f n)
285 "Args: (f n)
286 F is a function of an N-vector. Returns a function that approximates the
287 jacobian function iof F by a symmetric difference."
288 (let* ((h .0001)
289 (del (* h (values-list (column-list (identity-matrix n))))))
290 #'(lambda (b)
291 (let ((b+ (mapcar #'(lambda (x) (+ b x)) del))
292 (b- (mapcar #'(lambda (x) (- b x)) del)))
293 ;; FIXME:AJR BAD IDIOM FOLLOWS
294 (apply #'bind-columns (/ (- (mapcar f b+) (mapcar f b-)) (* 2 h)))))))
296 (defun nlreg (f y guess &optional
297 (epsilon .0001) (count-limit 20) weights verbose)
298 "Args: (mean-function y guess &optional
299 (epsilon .0001) (count-limit 20) weights verbose)
300 MEAN-FUNCTION returns the mean response vector for a given parameter vector.
301 Y is the observed response vector. Returns the estimated parameter vector
302 obtained by a Gauss-Newton algorithm that continues until the ITERATION-LIMIT
303 is reached or no component of the parameter vector changes by more than
304 EPSILON. The jacobian of MEAN-FUNCTION is approximated by a symmetric difference."
305 (nlreg1 f (make-jacobian f (length guess)) y guess
306 epsilon count-limit weights verbose))