Display SVG images in external <object> files
[emacs.git] / lisp / calc / calc-nlfit.el
blob67309fe8eecedaab5be78669929c2f85c7a964b1
1 ;;; calc-nlfit.el --- nonlinear curve fitting for Calc
3 ;; Copyright (C) 2007-2014 Free Software Foundation, Inc.
5 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
7 ;; This file is part of GNU Emacs.
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 of the License, or
12 ;; (at your option) any later version.
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.
19 ;; You should have received a copy of the GNU General Public License
20 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
22 ;;; Commentary:
24 ;; This code uses the Levenberg-Marquardt method, as described in
25 ;; _Numerical Analysis_ by H. R. Schwarz, to fit data to
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.
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.
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
39 ;; respective times t0, t1, ..., tn.
41 ;; The Hubbert Linearization, y/x=A*(1-x/B)
42 ;; Here, y is thought of as the rate of growth of a population
43 ;; and x represents the actual population. This is essentially
44 ;; the differential equation describing the actual population.
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
56 ;; formula, and then translating it to b and c. From this, we can
57 ;; also get approximations for the bell curve parameters.
59 ;;; Code:
61 (require 'calc-arith)
62 (require 'calcalg3)
64 ;; Declare functions which are defined elsewhere.
65 (declare-function calc-get-fit-variables "calcalg3" (nv nc &optional defv defc with-y homog))
66 (declare-function math-map-binop "calcalg3" (binop args1 args2))
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))
71 (s2data (if sdata
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)
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)))))
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
112 ;;; tdata. The following two functions will find pdata or qdata,
113 ;;; given the other..
115 ;;; First, given two lists; one of values q0, q1, ..., qn and one of
116 ;;; corresponding times t0, t1, ..., tn; return a list
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)).
123 (defun math-nlfit-get-rates-from-cumul (tdata qdata)
124 (let ((pdata (list
125 (math-div
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)))
158 ;;; Next, given two lists -- one of rates p0, p1, ..., pn and one of
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))
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
172 (math-mul
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)))
182 ;;; Given the qdata, pdata and tdata, find the parameters
183 ;;; a, b and c that fit q = a/(1+b*exp(c*t)).
184 ;;; a is found using the method described by de Sousa.
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.
190 ;;; In this case, a will be replaced by a number slightly above
191 ;;; the maximum value of q.
193 (defun math-nlfit-find-qmax (qdata pdata tdata)
194 (let* ((ratios (math-map-binop 'math-div pdata qdata))
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)))
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
210 (calcFunc-exp (nth 0 bandc))
211 (nth 1 bandc))))
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.
218 ;;; First, given pdata and tdata, find what half of qmax would be if q0=0.
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
227 (math-mul
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))
236 ;;; Next, given pdata and tdata, approximate q0.
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)))
242 (while (math-lessp (math-nlfit-find-qmax
243 (mapcar
244 (lambda (q) (math-add q0 q))
245 qdata)
246 pdata tdata)
247 (math-mul
248 '(float 5 -1)
249 (math-add
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)))
263 (if (math-lessp
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)))))
275 ;;; To improve the approximations to the parameters, we can use
276 ;;; Marquardt method as described in Schwarz's book.
278 ;;; Small numbers used in the Givens algorithm
279 (defvar math-nlfit-delta '(float 1 -8))
281 (defvar math-nlfit-epsilon '(float 1 -5))
283 ;;; Maximum number of iterations
284 (defvar math-nlfit-max-its 100)
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.
290 (defun math-nlfit-set-elt (vec i x)
291 (setcar (nthcdr (1- i) vec) x))
293 (defun math-nlfit-get-elt (vec i)
294 (nth (1- i) vec))
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)
301 (setq mat (cons (copy-sequence row) mat))
302 (setq k (1+ k)))
303 mat))
305 (defun math-nlfit-set-matx-elt (mat i j x)
306 (setcar (nthcdr (1- j) (nth (1- i) mat)) x))
308 (defun math-nlfit-get-matx-elt (mat i j)
309 (nth (1- j) (nth (1- i) mat)))
311 ;;; For solving the linearized system.
312 ;;; (The Givens method, from Schwarz.)
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))
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))
332 (if (math-lessp (calcFunc-abs cjj)
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)
340 (calcFunc-sqrt
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)))
354 (while (<= k n)
355 (let* ((cjk (math-nlfit-get-matx-elt C j k))
356 (cik (math-nlfit-get-matx-elt C i k))
357 (h (math-sub
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)))
379 (let ((i n)
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))))
389 (math-nlfit-set-elt x i
390 (math-neg
391 (math-div s
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
408 gamma (calcFunc-sqrt
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))
415 (rj (math-nlfit-get-elt r j))
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))))
430 (defun math-nlfit-jacobian (grad xlist parms &optional slist)
431 (let ((j nil))
432 (while xlist
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)))
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)))
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))
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)))))
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)))
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))
494 (setq xdata (cdr xdata))
495 (setq ydata (cdr ydata))
496 (setq sdata (cdr sdata)))
497 (reverse d)))
499 (defun math-nlfit-make-dtilda (d n)
500 (append d (make-list n 0)))
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))
519 (newparms (math-map-binop 'math-add (copy-tree parms) zeta))
520 (newchisq (math-nlfit-chi-sq xlist ylist newparms fn slist)))
521 (if (math-lessp newchisq chisq)
522 (progn
523 (if (math-lessp
524 (math-div
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)))
536 ;;; The functions that describe our models, and their gradients.
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))))))
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)))))
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)
554 (math-sqr
555 (math-add
556 1 ex)))))
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
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))))))
581 ;;; Functions to get the final covariance matrix and the sdevs
583 (defun math-nlfit-find-covar (grad xlist pparms)
584 (let ((j nil))
585 (while xlist
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)))
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)))
607 ;;; Now the Calc functions
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)))
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)))
627 ;;; Some functions to turn the parameter lists and variables
628 ;;; into the appropriate functions.
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 '+
637 (list '*
639 (calcFunc-exp
640 (list '*
642 var)))))))
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 '*
651 (calcFunc-exp
652 (list '*
654 (list '- var d))))
655 (list '^
656 (list '+
658 (calcFunc-exp
659 (list '*
661 (list '- var d))))
662 2))))
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))
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))
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))
690 (sigmacovar
691 (if sdevv
692 (math-nlfit-get-sigmas grad xdata finalparms (nth 0 fit))))
693 (sigmas
694 (if sdevv
695 (nth 0 sigmacovar)))
696 (finalparms
697 (if sigmas
698 (math-map-binop
699 (lambda (x y) (list 'sdev x y)) finalparms sigmas)
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
715 (list 'vec
716 soln
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))
724 (calcFunc-utpc (nth 0 fit)
725 (- n m))
726 '(var nan var-nan)))))
727 (math-nlfit-enter-result 1 "xfit" sln)))
729 (math-nlfit-enter-result 1 "fit" soln)))
730 (calc-record traillist "parm")))))
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))
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))
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))
759 (poverqdata (math-map-binop 'math-div pdata qdata))
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)
790 (math-mul a
791 (math-sub
793 (math-div x b))))
794 sdata)))
795 (setq sln
796 (list 'vec
797 soln
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))
810 (calcFunc-utpc
811 chisq
812 (- n 2))
813 '(var nan var-nan)))))
814 (math-nlfit-enter-result 1 "xfit" sln)))
816 (math-nlfit-enter-result 1 "fit" soln)))
817 (calc-record traillist "parm")))))
819 (provide 'calc-nlfit)