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