lisp/calc/calc-nlfit.el

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