1;;; calc-nlfit.el --- nonlinear curve fitting for Calc
2
3;; Copyright (C) 2007-2021 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 <https://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)
818