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