lisp/calc/calc-nlfit.el

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