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