Mercurial > emacs
annotate lisp/calc/calc-nlfit.el @ 87498:566965810b9e
*** empty log message ***
author | Magnus Henoch <mange@freemail.hu> |
---|---|
date | Mon, 31 Dec 2007 00:57:51 +0000 |
parents | 73907d4ce6a6 |
children | b9e8ab94c460 |
rev | line source |
---|---|
82260 | 1 ;;; calc-nlfit.el --- nonlinear curve fitting for Calc |
2 | |
3 ;; Copyright (C) 2007 Free Software Foundation, Inc. | |
4 | |
5 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com> | |
6 | |
7 ;; This file is part of GNU Emacs. | |
8 | |
9 ;; GNU Emacs is free software; you can redistribute it and/or modify | |
10 ;; it under the terms of the GNU General Public License as published by | |
11 ;; the Free Software Foundation; either version 3, or (at your option) | |
12 ;; any later version. | |
13 | |
14 ;; GNU Emacs is distributed in the hope that it will be useful, | |
15 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of | |
16 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
17 ;; GNU General Public License for more details. | |
18 | |
19 ;; You should have received a copy of the GNU General Public License | |
20 ;; along with GNU Emacs; see the file COPYING. If not, write to the | |
21 ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, | |
22 ;; Boston, MA 02110-1301, USA. | |
23 | |
24 ;;; Commentary: | |
25 | |
26 ;; This code uses the Levenberg-Marquardt method, as described in | |
27 ;; _Numerical Analysis_ by H. R. Schwarz, to fit data to | |
28 ;; nonlinear curves. Currently, the only the following curves are | |
29 ;; supported: | |
30 ;; The logistic S curve, y=a/(1+exp(b*(t-c))) | |
31 ;; Here, y is usually interpreted as the population of some | |
32 ;; quantity at time t. So we will think of the data as consisting | |
33 ;; of quantities q0, q1, ..., qn and their respective times | |
34 ;; t0, t1, ..., tn. | |
35 | |
36 ;; The logistic bell curve, y=A*exp(B*(t-C))/(1+exp(B*(t-C)))^2 | |
37 ;; Note that this is the derivative of the formula for the S curve. | |
38 ;; We get A=-a*b, B=b and C=c. Here, y is interpreted as the rate | |
39 ;; of growth of a population at time t. So we will think of the | |
40 ;; data as consisting of rates p0, p1, ..., pn and their | |
41 ;; respective times t0, t1, ..., tn. | |
42 | |
43 ;; The Hubbert Linearization, y/x=A*(1-x/B) | |
44 ;; Here, y is thought of as the rate of growth of a population | |
45 ;; and x represents the actual population. This is essentially | |
46 ;; the differential equation describing the actual population. | |
47 | |
48 ;; The Levenberg-Marquardt method is an iterative process: it takes | |
49 ;; an initial guess for the parameters and refines them. To get an | |
50 ;; initial guess for the parameters, we'll use a method described by | |
51 ;; Luis de Sousa in "Hubbert's Peak Mathematics". The idea is that | |
52 ;; given quantities Q and the corresponding rates P, they should | |
53 ;; satisfy P/Q= mQ+a. We can use the parameter a for an | |
54 ;; approximation for the parameter a in the S curve, and | |
55 ;; approximations for b and c are found using least squares on the | |
56 ;; linearization log((a/y)-1) = log(bb) + cc*t of | |
57 ;; y=a/(1+bb*exp(cc*t)), which is equivalent to the above s curve | |
58 ;; formula, and then tranlating it to b and c. From this, we can | |
59 ;; also get approximations for the bell curve parameters. | |
60 | |
61 ;;; Code: | |
62 | |
63 (require 'calc-arith) | |
86501
9c534c544c8c
(math-map-binop): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
86480
diff
changeset
|
64 (require 'calcalg3) |
82260 | 65 |
86480
9a67bab483db
(calc-get-fit-variables): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
82279
diff
changeset
|
66 ;; Declare functions which are defined elsewhere. |
9a67bab483db
(calc-get-fit-variables): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
82279
diff
changeset
|
67 (declare-function calc-get-fit-variables "calcalg3" (nv nc &optional defv defc with-y homog)) |
86506
73907d4ce6a6
(math-map-binop): Fix declaration.
Glenn Morris <rgm@gnu.org>
parents:
86501
diff
changeset
|
68 (declare-function math-map-binop "calcalg3" (binop args1 args2)) |
86480
9a67bab483db
(calc-get-fit-variables): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
82279
diff
changeset
|
69 |
82260 | 70 (defun math-nlfit-least-squares (xdata ydata &optional sdata sigmas) |
71 "Return the parameters A and B for the best least squares fit y=a+bx." | |
72 (let* ((n (length xdata)) | |
73 (s2data (if sdata | |
74 (mapcar 'calcFunc-sqr sdata) | |
75 (make-list n 1))) | |
76 (S (if sdata 0 n)) | |
77 (Sx 0) | |
78 (Sy 0) | |
79 (Sxx 0) | |
80 (Sxy 0) | |
81 D) | |
82 (while xdata | |
83 (let ((x (car xdata)) | |
84 (y (car ydata)) | |
85 (s (car s2data))) | |
86 (setq Sx (math-add Sx (if s (math-div x s) x))) | |
87 (setq Sy (math-add Sy (if s (math-div y s) y))) | |
88 (setq Sxx (math-add Sxx (if s (math-div (math-mul x x) s) | |
89 (math-mul x x)))) | |
90 (setq Sxy (math-add Sxy (if s (math-div (math-mul x y) s) | |
91 (math-mul x y)))) | |
92 (if sdata | |
93 (setq S (math-add S (math-div 1 s))))) | |
94 (setq xdata (cdr xdata)) | |
95 (setq ydata (cdr ydata)) | |
96 (setq s2data (cdr s2data))) | |
97 (setq D (math-sub (math-mul S Sxx) (math-mul Sx Sx))) | |
98 (let ((A (math-div (math-sub (math-mul Sxx Sy) (math-mul Sx Sxy)) D)) | |
99 (B (math-div (math-sub (math-mul S Sxy) (math-mul Sx Sy)) D))) | |
100 (if sigmas | |
101 (let ((C11 (math-div Sxx D)) | |
102 (C12 (math-neg (math-div Sx D))) | |
103 (C22 (math-div S D))) | |
104 (list (list 'sdev A (calcFunc-sqrt C11)) | |
105 (list 'sdev B (calcFunc-sqrt C22)) | |
106 (list 'vec | |
107 (list 'vec C11 C12) | |
108 (list 'vec C12 C22)))) | |
109 (list A B))))) | |
110 | |
111 ;;; The methods described by de Sousa require the cumulative data qdata | |
112 ;;; and the rates pdata. We will assume that we are given either | |
113 ;;; qdata and the corresponding times tdata, or pdata and the corresponding | |
114 ;;; tdata. The following two functions will find pdata or qdata, | |
115 ;;; given the other.. | |
116 | |
117 ;;; First, given two lists; one of values q0, q1, ..., qn and one of | |
118 ;;; corresponding times t0, t1, ..., tn; return a list | |
119 ;;; p0, p1, ..., pn of the rates of change of the qi with respect to t. | |
120 ;;; p0 is the right hand derivative (q1 - q0)/(t1 - t0). | |
121 ;;; pn is the left hand derivative (qn - q(n-1))/(tn - t(n-1)). | |
122 ;;; The other pis are the averages of the two: | |
123 ;;; (1/2)((qi - q(i-1))/(ti - t(i-1)) + (q(i+1) - qi)/(t(i+1) - ti)). | |
124 | |
125 (defun math-nlfit-get-rates-from-cumul (tdata qdata) | |
126 (let ((pdata (list | |
127 (math-div | |
128 (math-sub (nth 1 qdata) | |
129 (nth 0 qdata)) | |
130 (math-sub (nth 1 tdata) | |
131 (nth 0 tdata)))))) | |
132 (while (> (length qdata) 2) | |
133 (setq pdata | |
134 (cons | |
135 (math-mul | |
136 '(float 5 -1) | |
137 (math-add | |
138 (math-div | |
139 (math-sub (nth 2 qdata) | |
140 (nth 1 qdata)) | |
141 (math-sub (nth 2 tdata) | |
142 (nth 1 tdata))) | |
143 (math-div | |
144 (math-sub (nth 1 qdata) | |
145 (nth 0 qdata)) | |
146 (math-sub (nth 1 tdata) | |
147 (nth 0 tdata))))) | |
148 pdata)) | |
149 (setq qdata (cdr qdata))) | |
150 (setq pdata | |
151 (cons | |
152 (math-div | |
153 (math-sub (nth 1 qdata) | |
154 (nth 0 qdata)) | |
155 (math-sub (nth 1 tdata) | |
156 (nth 0 tdata))) | |
157 pdata)) | |
158 (reverse pdata))) | |
159 | |
160 ;;; Next, given two lists -- one of rates p0, p1, ..., pn and one of | |
161 ;;; corresponding times t0, t1, ..., tn -- and an initial values q0, | |
162 ;;; return a list q0, q1, ..., qn of the cumulative values. | |
163 ;;; q0 is the initial value given. | |
164 ;;; For i>0, qi is computed using the trapezoid rule: | |
165 ;;; qi = q(i-1) + (1/2)(pi + p(i-1))(ti - t(i-1)) | |
166 | |
167 (defun math-nlfit-get-cumul-from-rates (tdata pdata q0) | |
168 (let* ((qdata (list q0))) | |
169 (while (cdr pdata) | |
170 (setq qdata | |
171 (cons | |
172 (math-add (car qdata) | |
173 (math-mul | |
174 (math-mul | |
175 '(float 5 -1) | |
176 (math-add (nth 1 pdata) (nth 0 pdata))) | |
177 (math-sub (nth 1 tdata) | |
178 (nth 0 tdata)))) | |
179 qdata)) | |
180 (setq pdata (cdr pdata)) | |
181 (setq tdata (cdr tdata))) | |
182 (reverse qdata))) | |
183 | |
184 ;;; Given the qdata, pdata and tdata, find the parameters | |
185 ;;; a, b and c that fit q = a/(1+b*exp(c*t)). | |
186 ;;; a is found using the method described by de Sousa. | |
187 ;;; b and c are found using least squares on the linearization | |
188 ;;; log((a/q)-1) = log(b) + c*t | |
189 ;;; In some cases (where the logistic curve may well be the wrong | |
190 ;;; model), the computed a will be less than or equal to the maximum | |
191 ;;; value of q in qdata; in which case the above linearization won't work. | |
192 ;;; In this case, a will be replaced by a number slightly above | |
193 ;;; the maximum value of q. | |
194 | |
195 (defun math-nlfit-find-qmax (qdata pdata tdata) | |
86501
9c534c544c8c
(math-map-binop): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
86480
diff
changeset
|
196 (let* ((ratios (math-map-binop 'math-div pdata qdata)) |
82260 | 197 (lsdata (math-nlfit-least-squares ratios tdata)) |
198 (qmax (math-max-list (car qdata) (cdr qdata))) | |
199 (a (math-neg (math-div (nth 1 lsdata) (nth 0 lsdata))))) | |
200 (if (math-lessp a qmax) | |
201 (math-add '(float 5 -1) qmax) | |
202 a))) | |
203 | |
204 (defun math-nlfit-find-logistic-parameters (qdata pdata tdata) | |
205 (let* ((a (math-nlfit-find-qmax qdata pdata tdata)) | |
206 (newqdata | |
207 (mapcar (lambda (q) (calcFunc-ln (math-sub (math-div a q) 1))) | |
208 qdata)) | |
209 (bandc (math-nlfit-least-squares tdata newqdata))) | |
210 (list | |
211 a | |
212 (calcFunc-exp (nth 0 bandc)) | |
213 (nth 1 bandc)))) | |
214 | |
215 ;;; Next, given the pdata and tdata, we can find the qdata if we know q0. | |
216 ;;; We first try to find q0, using the fact that when p takes on its largest | |
217 ;;; value, q is half of its maximum value. So we'll find the maximum value | |
218 ;;; of q given various q0, and use bisection to approximate the correct q0. | |
219 | |
220 ;;; First, given pdata and tdata, find what half of qmax would be if q0=0. | |
221 | |
222 (defun math-nlfit-find-qmaxhalf (pdata tdata) | |
223 (let ((pmax (math-max-list (car pdata) (cdr pdata))) | |
224 (qmh 0)) | |
225 (while (math-lessp (car pdata) pmax) | |
226 (setq qmh | |
227 (math-add qmh | |
228 (math-mul | |
229 (math-mul | |
230 '(float 5 -1) | |
231 (math-add (nth 1 pdata) (nth 0 pdata))) | |
232 (math-sub (nth 1 tdata) | |
233 (nth 0 tdata))))) | |
234 (setq pdata (cdr pdata)) | |
235 (setq tdata (cdr tdata))) | |
236 qmh)) | |
237 | |
238 ;;; Next, given pdata and tdata, approximate q0. | |
239 | |
240 (defun math-nlfit-find-q0 (pdata tdata) | |
241 (let* ((qhalf (math-nlfit-find-qmaxhalf pdata tdata)) | |
242 (q0 (math-mul 2 qhalf)) | |
243 (qdata (math-nlfit-get-cumul-from-rates tdata pdata q0))) | |
244 (while (math-lessp (math-nlfit-find-qmax | |
245 (mapcar | |
246 (lambda (q) (math-add q0 q)) | |
247 qdata) | |
248 pdata tdata) | |
249 (math-mul | |
250 '(float 5 -1) | |
251 (math-add | |
252 q0 | |
253 qhalf))) | |
254 (setq q0 (math-add q0 qhalf))) | |
255 (let* ((qmin (math-sub q0 qhalf)) | |
256 (qmax q0) | |
257 (qt (math-nlfit-find-qmax | |
258 (mapcar | |
259 (lambda (q) (math-add q0 q)) | |
260 qdata) | |
261 pdata tdata)) | |
262 (i 0)) | |
263 (while (< i 10) | |
264 (setq q0 (math-mul '(float 5 -1) (math-add qmin qmax))) | |
265 (if (math-lessp | |
266 (math-nlfit-find-qmax | |
267 (mapcar | |
268 (lambda (q) (math-add q0 q)) | |
269 qdata) | |
270 pdata tdata) | |
271 (math-mul '(float 5 -1) (math-add qhalf q0))) | |
272 (setq qmin q0) | |
273 (setq qmax q0)) | |
274 (setq i (1+ i))) | |
275 (math-mul '(float 5 -1) (math-add qmin qmax))))) | |
276 | |
277 ;;; To improve the approximations to the parameters, we can use | |
278 ;;; Marquardt method as described in Schwarz's book. | |
279 | |
280 ;;; Small numbers used in the Givens algorithm | |
281 (defvar math-nlfit-delta '(float 1 -8)) | |
282 | |
283 (defvar math-nlfit-epsilon '(float 1 -5)) | |
284 | |
285 ;;; Maximum number of iterations | |
286 (defvar math-nlfit-max-its 100) | |
287 | |
288 ;;; Next, we need some functions for dealing with vectors and | |
289 ;;; matrices. For convenience, we'll work with Emacs lists | |
290 ;;; as vectors, rather than Calc's vectors. | |
291 | |
292 (defun math-nlfit-set-elt (vec i x) | |
293 (setcar (nthcdr (1- i) vec) x)) | |
294 | |
295 (defun math-nlfit-get-elt (vec i) | |
296 (nth (1- i) vec)) | |
297 | |
298 (defun math-nlfit-make-matrix (i j) | |
299 (let ((row (make-list j 0)) | |
300 (mat nil) | |
301 (k 0)) | |
302 (while (< k i) | |
86501
9c534c544c8c
(math-map-binop): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
86480
diff
changeset
|
303 (setq mat (cons (copy-sequence row) mat)) |
82260 | 304 (setq k (1+ k))) |
305 mat)) | |
306 | |
307 (defun math-nlfit-set-matx-elt (mat i j x) | |
308 (setcar (nthcdr (1- j) (nth (1- i) mat)) x)) | |
309 | |
310 (defun math-nlfit-get-matx-elt (mat i j) | |
311 (nth (1- j) (nth (1- i) mat))) | |
312 | |
313 ;;; For solving the linearized system. | |
314 ;;; (The Givens method, from Schwarz.) | |
315 | |
316 (defun math-nlfit-givens (C d) | |
317 (let* ((C (copy-tree C)) | |
318 (d (copy-tree d)) | |
319 (n (length (car C))) | |
320 (N (length C)) | |
321 (j 1) | |
322 (r (make-list N 0)) | |
323 (x (make-list N 0)) | |
324 w | |
325 gamma | |
326 sigma | |
327 rho) | |
328 (while (<= j n) | |
329 (let ((i (1+ j))) | |
330 (while (<= i N) | |
331 (let ((cij (math-nlfit-get-matx-elt C i j)) | |
332 (cjj (math-nlfit-get-matx-elt C j j))) | |
333 (when (not (math-equal 0 cij)) | |
334 (if (math-lessp (calcFunc-abs cjj) | |
335 (math-mul math-nlfit-delta (calcFunc-abs cij))) | |
336 (setq w (math-neg cij) | |
337 gamma 0 | |
338 sigma 1 | |
339 rho 1) | |
340 (setq w (math-mul | |
341 (calcFunc-sign cjj) | |
342 (calcFunc-sqrt | |
343 (math-add | |
344 (math-mul cjj cjj) | |
345 (math-mul cij cij)))) | |
346 gamma (math-div cjj w) | |
347 sigma (math-neg (math-div cij w))) | |
348 (if (math-lessp (calcFunc-abs sigma) gamma) | |
349 (setq rho sigma) | |
350 (setq rho (math-div (calcFunc-sign sigma) gamma)))) | |
351 (setq cjj w | |
352 cij rho) | |
353 (math-nlfit-set-matx-elt C j j w) | |
354 (math-nlfit-set-matx-elt C i j rho) | |
355 (let ((k (1+ j))) | |
356 (while (<= k n) | |
357 (let* ((cjk (math-nlfit-get-matx-elt C j k)) | |
358 (cik (math-nlfit-get-matx-elt C i k)) | |
359 (h (math-sub | |
360 (math-mul gamma cjk) (math-mul sigma cik)))) | |
361 (setq cik (math-add | |
362 (math-mul sigma cjk) | |
363 (math-mul gamma cik))) | |
364 (setq cjk h) | |
365 (math-nlfit-set-matx-elt C i k cik) | |
366 (math-nlfit-set-matx-elt C j k cjk) | |
367 (setq k (1+ k))))) | |
368 (let* ((di (math-nlfit-get-elt d i)) | |
369 (dj (math-nlfit-get-elt d j)) | |
370 (h (math-sub | |
371 (math-mul gamma dj) | |
372 (math-mul sigma di)))) | |
373 (setq di (math-add | |
374 (math-mul sigma dj) | |
375 (math-mul gamma di))) | |
376 (setq dj h) | |
377 (math-nlfit-set-elt d i di) | |
378 (math-nlfit-set-elt d j dj)))) | |
379 (setq i (1+ i)))) | |
380 (setq j (1+ j))) | |
82279
3f6660bf6595
(math-nlfit-curve): Remove unnecessary variables.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
82274
diff
changeset
|
381 (let ((i n) |
3f6660bf6595
(math-nlfit-curve): Remove unnecessary variables.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
82274
diff
changeset
|
382 s) |
82260 | 383 (while (>= i 1) |
384 (math-nlfit-set-elt r i 0) | |
385 (setq s (math-nlfit-get-elt d i)) | |
386 (let ((k (1+ i))) | |
387 (while (<= k n) | |
388 (setq s (math-add s (math-mul (math-nlfit-get-matx-elt C i k) | |
389 (math-nlfit-get-elt x k)))) | |
390 (setq k (1+ k)))) | |
391 (math-nlfit-set-elt x i | |
392 (math-neg | |
393 (math-div s | |
394 (math-nlfit-get-matx-elt C i i)))) | |
395 (setq i (1- i)))) | |
396 (let ((i (1+ n))) | |
397 (while (<= i N) | |
398 (math-nlfit-set-elt r i (math-nlfit-get-elt d i)) | |
399 (setq i (1+ i)))) | |
400 (let ((j n)) | |
401 (while (>= j 1) | |
402 (let ((i N)) | |
403 (while (>= i (1+ j)) | |
404 (setq rho (math-nlfit-get-matx-elt C i j)) | |
405 (if (math-equal rho 1) | |
406 (setq gamma 0 | |
407 sigma 1) | |
408 (if (math-lessp (calcFunc-abs rho) 1) | |
409 (setq sigma rho | |
410 gamma (calcFunc-sqrt | |
411 (math-sub 1 (math-mul sigma sigma)))) | |
412 (setq gamma (math-div 1 (calcFunc-abs rho)) | |
413 sigma (math-mul (calcFunc-sign rho) | |
414 (calcFunc-sqrt | |
415 (math-sub 1 (math-mul gamma gamma))))))) | |
416 (let ((ri (math-nlfit-get-elt r i)) | |
82279
3f6660bf6595
(math-nlfit-curve): Remove unnecessary variables.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
82274
diff
changeset
|
417 (rj (math-nlfit-get-elt r j)) |
3f6660bf6595
(math-nlfit-curve): Remove unnecessary variables.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
82274
diff
changeset
|
418 h) |
82260 | 419 (setq h (math-add (math-mul gamma rj) |
420 (math-mul sigma ri))) | |
421 (setq ri (math-sub | |
422 (math-mul gamma ri) | |
423 (math-mul sigma rj))) | |
424 (setq rj h) | |
425 (math-nlfit-set-elt r i ri) | |
426 (math-nlfit-set-elt r j rj)) | |
427 (setq i (1- i)))) | |
428 (setq j (1- j)))) | |
429 | |
430 x)) | |
431 | |
432 (defun math-nlfit-jacobian (grad xlist parms &optional slist) | |
433 (let ((j nil)) | |
434 (while xlist | |
435 (let ((row (apply grad (car xlist) parms))) | |
436 (setq j | |
437 (cons | |
438 (if slist | |
439 (mapcar (lambda (x) (math-div x (car slist))) row) | |
440 row) | |
441 j))) | |
442 (setq slist (cdr slist)) | |
443 (setq xlist (cdr xlist))) | |
444 (reverse j))) | |
445 | |
446 (defun math-nlfit-make-ident (l n) | |
447 (let ((m (math-nlfit-make-matrix n n)) | |
448 (i 1)) | |
449 (while (<= i n) | |
450 (math-nlfit-set-matx-elt m i i l) | |
451 (setq i (1+ i))) | |
452 m)) | |
453 | |
454 (defun math-nlfit-chi-sq (xlist ylist parms fn &optional slist) | |
455 (let ((cs 0)) | |
456 (while xlist | |
457 (let ((c | |
458 (math-sub | |
459 (apply fn (car xlist) parms) | |
460 (car ylist)))) | |
461 (if slist | |
462 (setq c (math-div c (car slist)))) | |
463 (setq cs | |
464 (math-add cs | |
465 (math-mul c c)))) | |
466 (setq xlist (cdr xlist)) | |
467 (setq ylist (cdr ylist)) | |
468 (setq slist (cdr slist))) | |
469 cs)) | |
470 | |
471 (defun math-nlfit-init-lambda (C) | |
472 (let ((l 0) | |
473 (n (length (car C))) | |
474 (N (length C))) | |
475 (while C | |
476 (let ((row (car C))) | |
477 (while row | |
478 (setq l (math-add l (math-mul (car row) (car row)))) | |
479 (setq row (cdr row)))) | |
480 (setq C (cdr C))) | |
481 (calcFunc-sqrt (math-div l (math-mul n N))))) | |
482 | |
483 (defun math-nlfit-make-Ctilda (C l) | |
484 (let* ((n (length (car C))) | |
485 (bot (math-nlfit-make-ident l n))) | |
486 (append C bot))) | |
487 | |
488 (defun math-nlfit-make-d (fn xdata ydata parms &optional sdata) | |
489 (let ((d nil)) | |
490 (while xdata | |
491 (setq d (cons | |
492 (let ((dd (math-sub (apply fn (car xdata) parms) | |
493 (car ydata)))) | |
494 (if sdata (math-div dd (car sdata)) dd)) | |
495 d)) | |
496 (setq xdata (cdr xdata)) | |
497 (setq ydata (cdr ydata)) | |
498 (setq sdata (cdr sdata))) | |
499 (reverse d))) | |
500 | |
501 (defun math-nlfit-make-dtilda (d n) | |
502 (append d (make-list n 0))) | |
503 | |
504 (defun math-nlfit-fit (xlist ylist parms fn grad &optional slist) | |
505 (let* | |
506 ((C (math-nlfit-jacobian grad xlist parms slist)) | |
507 (d (math-nlfit-make-d fn xlist ylist parms slist)) | |
508 (chisq (math-nlfit-chi-sq xlist ylist parms fn slist)) | |
509 (lambda (math-nlfit-init-lambda C)) | |
510 (really-done nil) | |
511 (iters 0)) | |
512 (while (and | |
513 (not really-done) | |
514 (< iters math-nlfit-max-its)) | |
515 (setq iters (1+ iters)) | |
516 (let ((done nil)) | |
517 (while (not done) | |
518 (let* ((Ctilda (math-nlfit-make-Ctilda C lambda)) | |
519 (dtilda (math-nlfit-make-dtilda d (length (car C)))) | |
520 (zeta (math-nlfit-givens Ctilda dtilda)) | |
86501
9c534c544c8c
(math-map-binop): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
86480
diff
changeset
|
521 (newparms (math-map-binop 'math-add (copy-tree parms) zeta)) |
82260 | 522 (newchisq (math-nlfit-chi-sq xlist ylist newparms fn slist))) |
523 (if (math-lessp newchisq chisq) | |
524 (progn | |
525 (if (math-lessp | |
526 (math-div | |
527 (math-sub chisq newchisq) newchisq) math-nlfit-epsilon) | |
528 (setq really-done t)) | |
529 (setq lambda (math-div lambda 10)) | |
530 (setq chisq newchisq) | |
531 (setq parms newparms) | |
532 (setq done t)) | |
533 (setq lambda (math-mul lambda 10))))) | |
534 (setq C (math-nlfit-jacobian grad xlist parms slist)) | |
535 (setq d (math-nlfit-make-d fn xlist ylist parms slist)))) | |
536 (list chisq parms))) | |
537 | |
538 ;;; The functions that describe our models, and their gradients. | |
539 | |
540 (defun math-nlfit-s-logistic-fn (x a b c) | |
541 (math-div a (math-add 1 (math-mul b (calcFunc-exp (math-mul c x)))))) | |
542 | |
543 (defun math-nlfit-s-logistic-grad (x a b c) | |
544 (let* ((ep (calcFunc-exp (math-mul c x))) | |
545 (d (math-add 1 (math-mul b ep))) | |
546 (d2 (math-mul d d))) | |
547 (list | |
548 (math-div 1 d) | |
549 (math-neg (math-div (math-mul a ep) d2)) | |
550 (math-neg (math-div (math-mul a (math-mul b (math-mul x ep))) d2))))) | |
551 | |
552 (defun math-nlfit-b-logistic-fn (x a c d) | |
553 (let ((ex (calcFunc-exp (math-mul c (math-sub x d))))) | |
554 (math-div | |
555 (math-mul a ex) | |
556 (math-sqr | |
557 (math-add | |
558 1 ex))))) | |
559 | |
560 (defun math-nlfit-b-logistic-grad (x a c d) | |
561 (let* ((ex (calcFunc-exp (math-mul c (math-sub x d)))) | |
562 (ex1 (math-add 1 ex)) | |
563 (xd (math-sub x d))) | |
564 (list | |
565 (math-div | |
566 ex | |
567 (math-sqr ex1)) | |
568 (math-sub | |
569 (math-div | |
570 (math-mul a (math-mul xd ex)) | |
571 (math-sqr ex1)) | |
572 (math-div | |
573 (math-mul 2 (math-mul a (math-mul xd (math-sqr ex)))) | |
574 (math-pow ex1 3))) | |
575 (math-sub | |
576 (math-div | |
577 (math-mul 2 (math-mul a (math-mul c (math-sqr ex)))) | |
578 (math-pow ex1 3)) | |
579 (math-div | |
580 (math-mul a (math-mul c ex)) | |
581 (math-sqr ex1)))))) | |
582 | |
583 ;;; Functions to get the final covariance matrix and the sdevs | |
584 | |
585 (defun math-nlfit-find-covar (grad xlist pparms) | |
586 (let ((j nil)) | |
587 (while xlist | |
588 (setq j (cons (cons 'vec (apply grad (car xlist) pparms)) j)) | |
589 (setq xlist (cdr xlist))) | |
590 (setq j (cons 'vec (reverse j))) | |
591 (setq j | |
592 (math-mul | |
593 (calcFunc-trn j) j)) | |
594 (calcFunc-inv j))) | |
595 | |
596 (defun math-nlfit-get-sigmas (grad xlist pparms chisq) | |
597 (let* ((sgs nil) | |
598 (covar (math-nlfit-find-covar grad xlist pparms)) | |
599 (n (1- (length covar))) | |
600 (N (length xlist)) | |
601 (i 1)) | |
602 (when (> N n) | |
603 (while (<= i n) | |
604 (setq sgs (cons (calcFunc-sqrt (nth i (nth i covar))) sgs)) | |
605 (setq i (1+ i))) | |
606 (setq sgs (reverse sgs))) | |
607 (list sgs covar))) | |
608 | |
609 ;;; Now the Calc functions | |
610 | |
611 (defun math-nlfit-s-logistic-params (xdata ydata) | |
612 (let ((pdata (math-nlfit-get-rates-from-cumul xdata ydata))) | |
613 (math-nlfit-find-logistic-parameters ydata pdata xdata))) | |
614 | |
615 (defun math-nlfit-b-logistic-params (xdata ydata) | |
616 (let* ((q0 (math-nlfit-find-q0 ydata xdata)) | |
617 (qdata (math-nlfit-get-cumul-from-rates xdata ydata q0)) | |
618 (abc (math-nlfit-find-logistic-parameters qdata ydata xdata)) | |
619 (B (nth 1 abc)) | |
620 (C (nth 2 abc)) | |
621 (A (math-neg | |
622 (math-mul | |
623 (nth 0 abc) | |
624 (math-mul B C)))) | |
625 (D (math-neg (math-div (calcFunc-ln B) C))) | |
626 (A (math-div A B))) | |
627 (list A C D))) | |
628 | |
629 ;;; Some functions to turn the parameter lists and variables | |
630 ;;; into the appropriate functions. | |
631 | |
632 (defun math-nlfit-s-logistic-solnexpr (pms var) | |
633 (let ((a (nth 0 pms)) | |
634 (b (nth 1 pms)) | |
635 (c (nth 2 pms))) | |
636 (list '/ a | |
637 (list '+ | |
638 1 | |
639 (list '* | |
640 b | |
641 (calcFunc-exp | |
642 (list '* | |
643 c | |
644 var))))))) | |
645 | |
646 (defun math-nlfit-b-logistic-solnexpr (pms var) | |
647 (let ((a (nth 0 pms)) | |
648 (c (nth 1 pms)) | |
649 (d (nth 2 pms))) | |
650 (list '/ | |
651 (list '* | |
652 a | |
653 (calcFunc-exp | |
654 (list '* | |
655 c | |
656 (list '- var d)))) | |
657 (list '^ | |
658 (list '+ | |
659 1 | |
660 (calcFunc-exp | |
661 (list '* | |
662 c | |
663 (list '- var d)))) | |
664 2)))) | |
665 | |
666 (defun math-nlfit-enter-result (n prefix vals) | |
667 (setq calc-aborted-prefix prefix) | |
668 (calc-pop-push-record-list n prefix vals) | |
669 (calc-handle-whys)) | |
670 | |
671 (defun math-nlfit-fit-curve (fn grad solnexpr initparms &optional sdv) | |
672 (calc-slow-wrapper | |
673 (let* ((sdevv (or (eq sdv 'calcFunc-efit) (eq sdv 'calcFunc-xfit))) | |
674 (calc-display-working-message nil) | |
675 (data (calc-top 1)) | |
676 (xdata (cdr (car (cdr data)))) | |
677 (ydata (cdr (car (cdr (cdr data))))) | |
678 (sdata (if (math-contains-sdev-p ydata) | |
679 (mapcar (lambda (x) (math-get-sdev x t)) ydata) | |
680 nil)) | |
681 (ydata (mapcar (lambda (x) (math-get-value x)) ydata)) | |
682 (calc-curve-varnames nil) | |
683 (calc-curve-coefnames nil) | |
684 (calc-curve-nvars 1) | |
685 (fitvars (calc-get-fit-variables 1 3)) | |
686 (var (nth 1 calc-curve-varnames)) | |
687 (parms (cdr calc-curve-coefnames)) | |
688 (parmguess | |
689 (funcall initparms xdata ydata)) | |
690 (fit (math-nlfit-fit xdata ydata parmguess fn grad sdata)) | |
691 (finalparms (nth 1 fit)) | |
692 (sigmacovar | |
693 (if sdevv | |
694 (math-nlfit-get-sigmas grad xdata finalparms (nth 0 fit)))) | |
695 (sigmas | |
696 (if sdevv | |
697 (nth 0 sigmacovar))) | |
698 (finalparms | |
699 (if sigmas | |
86501
9c534c544c8c
(math-map-binop): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
86480
diff
changeset
|
700 (math-map-binop |
9c534c544c8c
(math-map-binop): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
86480
diff
changeset
|
701 (lambda (x y) (list 'sdev x y)) finalparms sigmas) |
82260 | 702 finalparms)) |
703 (soln (funcall solnexpr finalparms var))) | |
704 (let ((calc-fit-to-trail t) | |
705 (traillist nil)) | |
706 (while parms | |
707 (setq traillist (cons (list 'calcFunc-eq (car parms) (car finalparms)) | |
708 traillist)) | |
709 (setq finalparms (cdr finalparms)) | |
710 (setq parms (cdr parms))) | |
711 (setq traillist (calc-normalize (cons 'vec (nreverse traillist)))) | |
712 (cond ((eq sdv 'calcFunc-efit) | |
713 (math-nlfit-enter-result 1 "efit" soln)) | |
714 ((eq sdv 'calcFunc-xfit) | |
715 (let (sln) | |
716 (setq sln | |
717 (list 'vec | |
718 soln | |
719 traillist | |
720 (nth 1 sigmacovar) | |
721 '(vec) | |
722 (nth 0 fit) | |
723 (let ((n (length xdata)) | |
724 (m (length finalparms))) | |
725 (if (and sdata (> n m)) | |
726 (calcFunc-utpc (nth 0 fit) | |
727 (- n m)) | |
728 '(var nan var-nan))))) | |
729 (math-nlfit-enter-result 1 "xfit" sln))) | |
730 (t | |
731 (math-nlfit-enter-result 1 "fit" soln))) | |
732 (calc-record traillist "parm"))))) | |
733 | |
734 (defun calc-fit-s-shaped-logistic-curve (arg) | |
735 (interactive "P") | |
736 (math-nlfit-fit-curve 'math-nlfit-s-logistic-fn | |
737 'math-nlfit-s-logistic-grad | |
738 'math-nlfit-s-logistic-solnexpr | |
739 'math-nlfit-s-logistic-params | |
740 arg)) | |
741 | |
742 (defun calc-fit-bell-shaped-logistic-curve (arg) | |
743 (interactive "P") | |
744 (math-nlfit-fit-curve 'math-nlfit-b-logistic-fn | |
745 'math-nlfit-b-logistic-grad | |
746 'math-nlfit-b-logistic-solnexpr | |
747 'math-nlfit-b-logistic-params | |
748 arg)) | |
749 | |
750 (defun calc-fit-hubbert-linear-curve (&optional sdv) | |
751 (calc-slow-wrapper | |
752 (let* ((sdevv (or (eq sdv 'calcFunc-efit) (eq sdv 'calcFunc-xfit))) | |
753 (calc-display-working-message nil) | |
754 (data (calc-top 1)) | |
755 (qdata (cdr (car (cdr data)))) | |
756 (pdata (cdr (car (cdr (cdr data))))) | |
757 (sdata (if (math-contains-sdev-p pdata) | |
758 (mapcar (lambda (x) (math-get-sdev x t)) pdata) | |
759 nil)) | |
760 (pdata (mapcar (lambda (x) (math-get-value x)) pdata)) | |
86501
9c534c544c8c
(math-map-binop): Declare as a function.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
86480
diff
changeset
|
761 (poverqdata (math-map-binop 'math-div pdata qdata)) |
82260 | 762 (parmvals (math-nlfit-least-squares qdata poverqdata sdata sdevv)) |
763 (finalparms (list (nth 0 parmvals) | |
764 (math-neg | |
765 (math-div (nth 0 parmvals) | |
766 (nth 1 parmvals))))) | |
767 (calc-curve-varnames nil) | |
768 (calc-curve-coefnames nil) | |
769 (calc-curve-nvars 1) | |
770 (fitvars (calc-get-fit-variables 1 2)) | |
771 (var (nth 1 calc-curve-varnames)) | |
772 (parms (cdr calc-curve-coefnames)) | |
773 (soln (list '* (nth 0 finalparms) | |
774 (list '- 1 | |
775 (list '/ var (nth 1 finalparms)))))) | |
776 (let ((calc-fit-to-trail t) | |
777 (traillist nil)) | |
778 (setq traillist | |
779 (list 'vec | |
780 (list 'calcFunc-eq (nth 0 parms) (nth 0 finalparms)) | |
781 (list 'calcFunc-eq (nth 1 parms) (nth 1 finalparms)))) | |
782 (cond ((eq sdv 'calcFunc-efit) | |
783 (math-nlfit-enter-result 1 "efit" soln)) | |
784 ((eq sdv 'calcFunc-xfit) | |
785 (let (sln | |
786 (chisq | |
787 (math-nlfit-chi-sq | |
788 qdata poverqdata | |
789 (list (nth 1 (nth 0 finalparms)) | |
790 (nth 1 (nth 1 finalparms))) | |
791 (lambda (x a b) | |
792 (math-mul a | |
793 (math-sub | |
794 1 | |
795 (math-div x b)))) | |
796 sdata))) | |
797 (setq sln | |
798 (list 'vec | |
799 soln | |
800 traillist | |
801 (nth 2 parmvals) | |
802 (list | |
803 'vec | |
804 '(calcFunc-fitdummy 1) | |
805 (list 'calcFunc-neg | |
806 (list '/ | |
807 '(calcFunc-fitdummy 1) | |
808 '(calcFunc-fitdummy 2)))) | |
809 chisq | |
810 (let ((n (length qdata))) | |
811 (if (and sdata (> n 2)) | |
812 (calcFunc-utpc | |
813 chisq | |
814 (- n 2)) | |
815 '(var nan var-nan))))) | |
816 (math-nlfit-enter-result 1 "xfit" sln))) | |
817 (t | |
818 (math-nlfit-enter-result 1 "fit" soln))) | |
819 (calc-record traillist "parm"))))) | |
820 | |
821 (provide 'calc-nlfit) | |
82274 | 822 |
823 ;; arch-tag: 6eba3cd6-f48b-4a84-8174-10c15a024928 |