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