lisp/calc/calc-nlfit.el

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