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
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
16 :lisp-stat-object-system
19 :lisp-stat-compound-data
24 :lisp-stat-regression-linear
)
25 (:shadowing-import-from
:lisp-stat-object-system
26 slot-value call-method call-next-method
)
27 (:shadowing-import-from
:lisp-stat-math
28 expt
+ -
* / ** mod rem abs
1+ 1- log exp sqrt sin cos tan
29 asin acos atan sinh cosh tanh asinh acosh atanh float random
30 truncate floor ceiling round minusp zerop plusp evenp oddp
31 < <= = /= >= > ;; complex
32 conjugate realpart imagpart phase
33 min max logand logior logxor lognot ffloor fceiling
34 ftruncate fround signum cis
)
36 (:export nreg-model nreg-model-proto mean-function theta-hat epsilon
39 (in-package :lisp-stat-regression-nonlin
)
42 ;;; Nonlinear Regression Model Prototype
45 (defproto nreg-model-proto
46 '(mean-function theta-hat epsilon count-limit verbose
)
48 regression-model-proto
)
50 (defun nreg-model (mean-function y theta
61 "Args: (mean-function y theta &key (epsilon .0001) (count-limit 20)
62 (print t) parameter-names response-name case-labels
63 weights included (vetbose print))
64 Fits nonlinear regression model with MEAN-FUNCTION and response Y using initial
65 parameter guess THETA. Returns model object."
66 (let ((m (send nreg-model-proto
:new
)))
67 (send m
:mean-function mean-function
)
69 (send m
:new-initial-guess theta
)
70 (send m
:epsilon epsilon
)
71 (send m
:count-limit count-limit
)
72 (send m
:parameter-names parameter-names
)
73 (send m
:response-name response-name
)
74 (send m
:case-labels case-labels
)
75 (send m
:weights weights
)
76 (send m
:included included
)
77 (send m
:verbose verbose
)
78 (if print
(send m
:display
))
81 (defmeth nreg-model-proto
:save
()
83 Returns an expression that will reconstruct the regression model."
84 `(nreg-model ',(send self
:mean-function
)
86 ',(send self
:coef-estimates
)
87 :epsilon
',(send self
:epsilon
)
88 :count-limit
',(send self
:count-limit
)
89 :predictor-names
',(send self
:predictor-names
)
90 :response-name
',(send self
:response-name
)
91 :case-labels
',(send self
:case-labels
)
92 :weights
',(send self
:weights
)
93 :included
',(send self
:included
)
94 :verbose
',(send self
:verbose
)))
100 (defmeth nreg-model-proto
:compute
()
102 Recomputes the estimates. For internal use by other messages"
103 (let* ((y (send self
:y
))
104 (weights (send self
:weights
))
105 (inc (if-else (send self
:included
) 1 0))
106 (w (if weights
(* inc weights
) inc
)))
107 (setf (slot-value 'theta-hat
)
108 (nlreg (send self
:mean-function
)
110 (slot-value 'theta-hat
)
112 (send self
:count-limit
)
114 (send self
:verbose
)))
115 (setf (slot-value 'x
)
116 (funcall (make-jacobian (slot-value 'mean-function
)
117 (length (slot-value 'theta-hat
)))
118 (slot-value 'theta-hat
)))
119 (setf (slot-value 'intercept
) nil
)
121 (let ((r (send self
:residuals
)))
122 (setf (slot-value 'residual-sum-of-squares
)
123 (sum (* inc r r
))))))
126 ;;; Slot Accessors and Mutators
129 (defmeth nreg-model-proto
:new-initial-guess
(guess)
130 "Message args: (guess)
131 Sets a new initial uess for parmeters."
132 (setf (slot-value 'theta-hat
) guess
)
133 (send self
:needs-computing t
))
135 (defmeth nreg-model-proto
:theta-hat
()
137 Returns current parameter estimate."
138 (if (send self
:needs-computing
) (send self
:compute
))
139 (coerce (slot-value 'theta-hat
) 'list
))
141 (defmeth nreg-model-proto
:mean-function
(&optional f
)
142 "Message args: (&optional f)
143 With no argument returns the mean function as supplied to m. With an
144 argument F sets the mean function of m to F and recomputes the
146 (when (and f
(functionp f
))
147 (setf (slot-value 'mean-function
) f
)
148 (send self
:needs-computing t
))
149 (slot-value 'mean-function
))
151 (defmeth nreg-model-proto
:epsilon
(&optional eps
)
152 "Message args: (&optional eps)
153 With no argument returns the tolerance as supplied to m. With an argument
154 EPS sets the tolerance of m to EPS and recomputes the estimates."
155 (when (and eps
(numberp eps
))
156 (setf (slot-value 'epsilon
) eps
)
157 (send self
:needs-computing t
))
158 (slot-value 'epsilon
))
160 (defmeth nreg-model-proto
:count-limit
(&optional count
)
161 "Message args: (&optional new-count)
162 With no argument returns the iteration count limit as supplied to m. With
163 an argument COUNT sets the limit to COUNT and recomputes the
165 (when (and count
(numberp count
))
166 (setf (slot-value 'count-limit
) count
)
167 (send self
:needs-computing t
))
168 (slot-value 'count-limit
))
170 (defmeth nreg-model-proto
:parameter-names
(&optional
(names nil set
))
171 "Method args: (&optional names)
172 Sets or returns parameter names."
173 (if set
(setf (slot-value 'predictor-names
) names
))
174 (let ((p-names (slot-value 'predictor-names
))
175 (p (length (slot-value 'theta-hat
))))
176 (if (not (and p-names
(= p
(length p-names
))))
177 (setf (slot-value 'predictor-names
)
178 (mapcar #'(lambda (a) (format nil
"Parameter ~a" a
))
180 (slot-value 'predictor-names
))
182 (defmeth nreg-model-proto
:verbose
(&optional
(val nil set
))
183 "Method args: (&optional val)
184 Sets or retrieves verbose setting. If T iteration info is printed during
186 (if set
(setf (slot-value 'verbose
) val
))
187 (slot-value 'verbose
))
190 ;;; Overrides for Linear Regression Methods
193 (defmeth nreg-model-proto
:x
()
195 Returns the Jacobian matrix at theta-hat."
198 (defmeth nreg-model-proto
:intercept
(&rest args
)
200 Always returns nil. (For compatibility with linear regression.)"
203 (defmeth nreg-model-proto
:fit-values
()
205 Returns the fitted values for the model."
206 (coerce (funcall (send self
:mean-function
) (send self
:theta-hat
))
209 (defmeth nreg-model-proto
:coef-estimates
(&optional guess
)
210 "Message args: (&optional guess)
211 With no argument returns the current parameter estimate. With an
212 argument GUESS takes it as a new initial guess and recomputes
214 (if guess
(send self
:new-initial-guess guess
))
215 (send self
:theta-hat
))
217 (defmeth nreg-model-proto
:predictor-names
() (send self
:parameter-names
))
221 ;;;; Linear Regression Coefficients
225 (defun regression-coefficients (x y
&key
(intercept T
) weights
)
226 "Args: (x y &key (intercept T) weights)
227 Returns the coefficients of the regression of the sequence Y on the columns of
229 (let* ((m (if weights
230 (make-sweep-matrix x y weights
)
231 (make-sweep-matrix x y
)))
232 (n (array-dimension x
1)))
233 (coerce (compound-data-seq
235 (select (car (sweep-operator m
(iseq 1 n
)))
238 (select (car (sweep-operator m
(iseq 0 n
)))
245 ;;;; Nonlinear Regression Functions
248 (defun nlreg1 (f j y initial-beta epsilon count-limit weights verbose
)
249 "Args: (mean-function jacobian y initial-beta
250 epsilon count-limit weights verbose)
251 MEAN-FUNCTION returns the mean response vector for a given parameter vector.
252 JACOBIAN returns the jacobian of MEAN-FUNCTION at a given parameter vector.
253 Y is the observed response vector. Returns the estimated parameter vector
254 obtained by a Gauss-Newton algorithm with backtracking that continues until
255 the COUNT-LIMIT is reached or no component of the parameter vector changes
256 by more than EPSILON."
257 (labels ((rss (beta) ; residual sum of squares
258 (let ((res (- y
(funcall f beta
))))
259 (sum (if weights
(* res res weights
) (* res res
)))))
260 (next-beta (beta delta rss
) ; next beta by backtracking
261 (do* ((lambda 1 (/ lambda
2))
262 (new-rss (rss (+ beta delta
))
263 (rss (+ beta
(* lambda delta
)))))
264 ((or (< new-rss rss
) (< lambda
.0001))
265 (if (>= lambda
.0001)
266 (+ beta
(* lambda delta
))
268 (do* ((delbeta (regression-coefficients
269 (funcall j initial-beta
)
270 (- y
(funcall f initial-beta
))
273 (regression-coefficients
275 (- y
(funcall f beta
))
278 (beta initial-beta
(next-beta beta delbeta rss
))
279 (rss (rss beta
) (rss beta
))
280 (count 0 (1+ count
)))
281 ((or (> count count-limit
) (> epsilon
(max (abs delbeta
))))
282 (if (and verbose
(> count count-limit
))
283 (format t
"Iteration limit exceeded.~%"))
285 (if verbose
(format t
"Residual sum of squares: ~10g~%" rss
)))))
287 (defun make-jacobian (f n
)
289 F is a function of an N-vector. Returns a function that approximates the
290 jacobian function iof F by a symmetric difference."
292 (del (* h
(values-list (column-list (identity-matrix n
))))))
294 (let ((b+ (mapcar #'(lambda (x) (+ b x
)) del
))
295 (b- (mapcar #'(lambda (x) (- b x
)) del
)))
296 ;; FIXME:AJR BAD IDIOM FOLLOWS
297 (apply #'bind-columns
(/ (- (mapcar f b
+) (mapcar f b-
)) (* 2 h
)))))))
299 (defun nlreg (f y guess
&optional
300 (epsilon .0001) (count-limit 20) weights verbose
)
301 "Args: (mean-function y guess &optional
302 (epsilon .0001) (count-limit 20) weights verbose)
303 MEAN-FUNCTION returns the mean response vector for a given parameter vector.
304 Y is the observed response vector. Returns the estimated parameter vector
305 obtained by a Gauss-Newton algorithm that continues until the ITERATION-LIMIT
306 is reached or no component of the parameter vector changes by more than
307 EPSILON. The jacobian of MEAN-FUNCTION is approximated by a symmetric difference."
308 (nlreg1 f
(make-jacobian f
(length guess
)) y guess
309 epsilon count-limit weights verbose
))