lisp/calc/calc-nlfit.el

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