[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Emacs-diffs] Changes to calc-nlfit.el
|
From: |
Jay Belanger |
|
Subject: |
[Emacs-diffs] Changes to calc-nlfit.el |
|
Date: |
Sat, 04 Aug 2007 03:54:14 +0000 |
CVSROOT: /cvsroot/emacs
Module name: emacs
Changes by: Jay Belanger <jpb> 07/08/04 03:54:14
Index: calc-nlfit.el
===================================================================
RCS file: calc-nlfit.el
diff -N calc-nlfit.el
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ calc-nlfit.el 4 Aug 2007 03:54:14 -0000 1.1
@@ -0,0 +1,815 @@
+;;; calc-nlfit.el --- nonlinear curve fitting for Calc
+
+;; Copyright (C) 2007 Free Software Foundation, Inc.
+
+;; Maintainer: Jay Belanger <address@hidden>
+
+;; This file is part of GNU Emacs.
+
+;; GNU Emacs is free software; you can redistribute it and/or modify
+;; it under the terms of the GNU General Public License as published by
+;; the Free Software Foundation; either version 3, or (at your option)
+;; any later version.
+
+;; GNU Emacs is distributed in the hope that it will be useful,
+;; but WITHOUT ANY WARRANTY; without even the implied warranty of
+;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+;; GNU General Public License for more details.
+
+;; You should have received a copy of the GNU General Public License
+;; along with GNU Emacs; see the file COPYING. If not, write to the
+;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+;; Boston, MA 02110-1301, USA.
+
+;;; Commentary:
+
+;; This code uses the Levenberg-Marquardt method, as described in
+;; _Numerical Analysis_ by H. R. Schwarz, to fit data to
+;; nonlinear curves. Currently, the only the following curves are
+;; supported:
+;; The logistic S curve, y=a/(1+exp(b*(t-c)))
+;; Here, y is usually interpreted as the population of some
+;; quantity at time t. So we will think of the data as consisting
+;; of quantities q0, q1, ..., qn and their respective times
+;; t0, t1, ..., tn.
+
+;; The logistic bell curve, y=A*exp(B*(t-C))/(1+exp(B*(t-C)))^2
+;; Note that this is the derivative of the formula for the S curve.
+;; We get A=-a*b, B=b and C=c. Here, y is interpreted as the rate
+;; of growth of a population at time t. So we will think of the
+;; data as consisting of rates p0, p1, ..., pn and their
+;; respective times t0, t1, ..., tn.
+
+;; The Hubbert Linearization, y/x=A*(1-x/B)
+;; Here, y is thought of as the rate of growth of a population
+;; and x represents the actual population. This is essentially
+;; the differential equation describing the actual population.
+
+;; The Levenberg-Marquardt method is an iterative process: it takes
+;; an initial guess for the parameters and refines them. To get an
+;; initial guess for the parameters, we'll use a method described by
+;; Luis de Sousa in "Hubbert's Peak Mathematics". The idea is that
+;; given quantities Q and the corresponding rates P, they should
+;; satisfy P/Q= mQ+a. We can use the parameter a for an
+;; approximation for the parameter a in the S curve, and
+;; approximations for b and c are found using least squares on the
+;; linearization log((a/y)-1) = log(bb) + cc*t of
+;; y=a/(1+bb*exp(cc*t)), which is equivalent to the above s curve
+;; formula, and then tranlating it to b and c. From this, we can
+;; also get approximations for the bell curve parameters.
+
+;;; Code:
+
+(require 'calc-arith)
+
+(defun math-nlfit-least-squares (xdata ydata &optional sdata sigmas)
+ "Return the parameters A and B for the best least squares fit y=a+bx."
+ (let* ((n (length xdata))
+ (s2data (if sdata
+ (mapcar 'calcFunc-sqr sdata)
+ (make-list n 1)))
+ (S (if sdata 0 n))
+ (Sx 0)
+ (Sy 0)
+ (Sxx 0)
+ (Sxy 0)
+ D)
+ (while xdata
+ (let ((x (car xdata))
+ (y (car ydata))
+ (s (car s2data)))
+ (setq Sx (math-add Sx (if s (math-div x s) x)))
+ (setq Sy (math-add Sy (if s (math-div y s) y)))
+ (setq Sxx (math-add Sxx (if s (math-div (math-mul x x) s)
+ (math-mul x x))))
+ (setq Sxy (math-add Sxy (if s (math-div (math-mul x y) s)
+ (math-mul x y))))
+ (if sdata
+ (setq S (math-add S (math-div 1 s)))))
+ (setq xdata (cdr xdata))
+ (setq ydata (cdr ydata))
+ (setq s2data (cdr s2data)))
+ (setq D (math-sub (math-mul S Sxx) (math-mul Sx Sx)))
+ (let ((A (math-div (math-sub (math-mul Sxx Sy) (math-mul Sx Sxy)) D))
+ (B (math-div (math-sub (math-mul S Sxy) (math-mul Sx Sy)) D)))
+ (if sigmas
+ (let ((C11 (math-div Sxx D))
+ (C12 (math-neg (math-div Sx D)))
+ (C22 (math-div S D)))
+ (list (list 'sdev A (calcFunc-sqrt C11))
+ (list 'sdev B (calcFunc-sqrt C22))
+ (list 'vec
+ (list 'vec C11 C12)
+ (list 'vec C12 C22))))
+ (list A B)))))
+
+;;; The methods described by de Sousa require the cumulative data qdata
+;;; and the rates pdata. We will assume that we are given either
+;;; qdata and the corresponding times tdata, or pdata and the corresponding
+;;; tdata. The following two functions will find pdata or qdata,
+;;; given the other..
+
+;;; First, given two lists; one of values q0, q1, ..., qn and one of
+;;; corresponding times t0, t1, ..., tn; return a list
+;;; p0, p1, ..., pn of the rates of change of the qi with respect to t.
+;;; p0 is the right hand derivative (q1 - q0)/(t1 - t0).
+;;; pn is the left hand derivative (qn - q(n-1))/(tn - t(n-1)).
+;;; The other pis are the averages of the two:
+;;; (1/2)((qi - q(i-1))/(ti - t(i-1)) + (q(i+1) - qi)/(t(i+1) - ti)).
+
+(defun math-nlfit-get-rates-from-cumul (tdata qdata)
+ (let ((pdata (list
+ (math-div
+ (math-sub (nth 1 qdata)
+ (nth 0 qdata))
+ (math-sub (nth 1 tdata)
+ (nth 0 tdata))))))
+ (while (> (length qdata) 2)
+ (setq pdata
+ (cons
+ (math-mul
+ '(float 5 -1)
+ (math-add
+ (math-div
+ (math-sub (nth 2 qdata)
+ (nth 1 qdata))
+ (math-sub (nth 2 tdata)
+ (nth 1 tdata)))
+ (math-div
+ (math-sub (nth 1 qdata)
+ (nth 0 qdata))
+ (math-sub (nth 1 tdata)
+ (nth 0 tdata)))))
+ pdata))
+ (setq qdata (cdr qdata)))
+ (setq pdata
+ (cons
+ (math-div
+ (math-sub (nth 1 qdata)
+ (nth 0 qdata))
+ (math-sub (nth 1 tdata)
+ (nth 0 tdata)))
+ pdata))
+ (reverse pdata)))
+
+;;; Next, given two lists -- one of rates p0, p1, ..., pn and one of
+;;; corresponding times t0, t1, ..., tn -- and an initial values q0,
+;;; return a list q0, q1, ..., qn of the cumulative values.
+;;; q0 is the initial value given.
+;;; For i>0, qi is computed using the trapezoid rule:
+;;; qi = q(i-1) + (1/2)(pi + p(i-1))(ti - t(i-1))
+
+(defun math-nlfit-get-cumul-from-rates (tdata pdata q0)
+ (let* ((qdata (list q0)))
+ (while (cdr pdata)
+ (setq qdata
+ (cons
+ (math-add (car qdata)
+ (math-mul
+ (math-mul
+ '(float 5 -1)
+ (math-add (nth 1 pdata) (nth 0 pdata)))
+ (math-sub (nth 1 tdata)
+ (nth 0 tdata))))
+ qdata))
+ (setq pdata (cdr pdata))
+ (setq tdata (cdr tdata)))
+ (reverse qdata)))
+
+;;; Given the qdata, pdata and tdata, find the parameters
+;;; a, b and c that fit q = a/(1+b*exp(c*t)).
+;;; a is found using the method described by de Sousa.
+;;; b and c are found using least squares on the linearization
+;;; log((a/q)-1) = log(b) + c*t
+;;; In some cases (where the logistic curve may well be the wrong
+;;; model), the computed a will be less than or equal to the maximum
+;;; value of q in qdata; in which case the above linearization won't work.
+;;; In this case, a will be replaced by a number slightly above
+;;; the maximum value of q.
+
+(defun math-nlfit-find-qmax (qdata pdata tdata)
+ (let* ((ratios (mapcar* 'math-div pdata qdata))
+ (lsdata (math-nlfit-least-squares ratios tdata))
+ (qmax (math-max-list (car qdata) (cdr qdata)))
+ (a (math-neg (math-div (nth 1 lsdata) (nth 0 lsdata)))))
+ (if (math-lessp a qmax)
+ (math-add '(float 5 -1) qmax)
+ a)))
+
+(defun math-nlfit-find-logistic-parameters (qdata pdata tdata)
+ (let* ((a (math-nlfit-find-qmax qdata pdata tdata))
+ (newqdata
+ (mapcar (lambda (q) (calcFunc-ln (math-sub (math-div a q) 1)))
+ qdata))
+ (bandc (math-nlfit-least-squares tdata newqdata)))
+ (list
+ a
+ (calcFunc-exp (nth 0 bandc))
+ (nth 1 bandc))))
+
+;;; Next, given the pdata and tdata, we can find the qdata if we know q0.
+;;; We first try to find q0, using the fact that when p takes on its largest
+;;; value, q is half of its maximum value. So we'll find the maximum value
+;;; of q given various q0, and use bisection to approximate the correct q0.
+
+;;; First, given pdata and tdata, find what half of qmax would be if q0=0.
+
+(defun math-nlfit-find-qmaxhalf (pdata tdata)
+ (let ((pmax (math-max-list (car pdata) (cdr pdata)))
+ (qmh 0))
+ (while (math-lessp (car pdata) pmax)
+ (setq qmh
+ (math-add qmh
+ (math-mul
+ (math-mul
+ '(float 5 -1)
+ (math-add (nth 1 pdata) (nth 0 pdata)))
+ (math-sub (nth 1 tdata)
+ (nth 0 tdata)))))
+ (setq pdata (cdr pdata))
+ (setq tdata (cdr tdata)))
+ qmh))
+
+;;; Next, given pdata and tdata, approximate q0.
+
+(defun math-nlfit-find-q0 (pdata tdata)
+ (let* ((qhalf (math-nlfit-find-qmaxhalf pdata tdata))
+ (q0 (math-mul 2 qhalf))
+ (qdata (math-nlfit-get-cumul-from-rates tdata pdata q0)))
+ (while (math-lessp (math-nlfit-find-qmax
+ (mapcar
+ (lambda (q) (math-add q0 q))
+ qdata)
+ pdata tdata)
+ (math-mul
+ '(float 5 -1)
+ (math-add
+ q0
+ qhalf)))
+ (setq q0 (math-add q0 qhalf)))
+ (let* ((qmin (math-sub q0 qhalf))
+ (qmax q0)
+ (qt (math-nlfit-find-qmax
+ (mapcar
+ (lambda (q) (math-add q0 q))
+ qdata)
+ pdata tdata))
+ (i 0))
+ (while (< i 10)
+ (setq q0 (math-mul '(float 5 -1) (math-add qmin qmax)))
+ (if (math-lessp
+ (math-nlfit-find-qmax
+ (mapcar
+ (lambda (q) (math-add q0 q))
+ qdata)
+ pdata tdata)
+ (math-mul '(float 5 -1) (math-add qhalf q0)))
+ (setq qmin q0)
+ (setq qmax q0))
+ (setq i (1+ i)))
+ (math-mul '(float 5 -1) (math-add qmin qmax)))))
+
+;;; To improve the approximations to the parameters, we can use
+;;; Marquardt method as described in Schwarz's book.
+
+;;; Small numbers used in the Givens algorithm
+(defvar math-nlfit-delta '(float 1 -8))
+
+(defvar math-nlfit-epsilon '(float 1 -5))
+
+;;; Maximum number of iterations
+(defvar math-nlfit-max-its 100)
+
+;;; Next, we need some functions for dealing with vectors and
+;;; matrices. For convenience, we'll work with Emacs lists
+;;; as vectors, rather than Calc's vectors.
+
+(defun math-nlfit-set-elt (vec i x)
+ (setcar (nthcdr (1- i) vec) x))
+
+(defun math-nlfit-get-elt (vec i)
+ (nth (1- i) vec))
+
+(defun math-nlfit-make-matrix (i j)
+ (let ((row (make-list j 0))
+ (mat nil)
+ (k 0))
+ (while (< k i)
+ (setq mat (cons (copy-list row) mat))
+ (setq k (1+ k)))
+ mat))
+
+(defun math-nlfit-set-matx-elt (mat i j x)
+ (setcar (nthcdr (1- j) (nth (1- i) mat)) x))
+
+(defun math-nlfit-get-matx-elt (mat i j)
+ (nth (1- j) (nth (1- i) mat)))
+
+;;; For solving the linearized system.
+;;; (The Givens method, from Schwarz.)
+
+(defun math-nlfit-givens (C d)
+ (let* ((C (copy-tree C))
+ (d (copy-tree d))
+ (n (length (car C)))
+ (N (length C))
+ (j 1)
+ (r (make-list N 0))
+ (x (make-list N 0))
+ w
+ gamma
+ sigma
+ rho)
+ (while (<= j n)
+ (let ((i (1+ j)))
+ (while (<= i N)
+ (let ((cij (math-nlfit-get-matx-elt C i j))
+ (cjj (math-nlfit-get-matx-elt C j j)))
+ (when (not (math-equal 0 cij))
+ (if (math-lessp (calcFunc-abs cjj)
+ (math-mul math-nlfit-delta (calcFunc-abs cij)))
+ (setq w (math-neg cij)
+ gamma 0
+ sigma 1
+ rho 1)
+ (setq w (math-mul
+ (calcFunc-sign cjj)
+ (calcFunc-sqrt
+ (math-add
+ (math-mul cjj cjj)
+ (math-mul cij cij))))
+ gamma (math-div cjj w)
+ sigma (math-neg (math-div cij w)))
+ (if (math-lessp (calcFunc-abs sigma) gamma)
+ (setq rho sigma)
+ (setq rho (math-div (calcFunc-sign sigma) gamma))))
+ (setq cjj w
+ cij rho)
+ (math-nlfit-set-matx-elt C j j w)
+ (math-nlfit-set-matx-elt C i j rho)
+ (let ((k (1+ j)))
+ (while (<= k n)
+ (let* ((cjk (math-nlfit-get-matx-elt C j k))
+ (cik (math-nlfit-get-matx-elt C i k))
+ (h (math-sub
+ (math-mul gamma cjk) (math-mul sigma cik))))
+ (setq cik (math-add
+ (math-mul sigma cjk)
+ (math-mul gamma cik)))
+ (setq cjk h)
+ (math-nlfit-set-matx-elt C i k cik)
+ (math-nlfit-set-matx-elt C j k cjk)
+ (setq k (1+ k)))))
+ (let* ((di (math-nlfit-get-elt d i))
+ (dj (math-nlfit-get-elt d j))
+ (h (math-sub
+ (math-mul gamma dj)
+ (math-mul sigma di))))
+ (setq di (math-add
+ (math-mul sigma dj)
+ (math-mul gamma di)))
+ (setq dj h)
+ (math-nlfit-set-elt d i di)
+ (math-nlfit-set-elt d j dj))))
+ (setq i (1+ i))))
+ (setq j (1+ j)))
+ (let ((i n))
+ (while (>= i 1)
+ (math-nlfit-set-elt r i 0)
+ (setq s (math-nlfit-get-elt d i))
+ (let ((k (1+ i)))
+ (while (<= k n)
+ (setq s (math-add s (math-mul (math-nlfit-get-matx-elt C i k)
+ (math-nlfit-get-elt x k))))
+ (setq k (1+ k))))
+ (math-nlfit-set-elt x i
+ (math-neg
+ (math-div s
+ (math-nlfit-get-matx-elt C i i))))
+ (setq i (1- i))))
+ (let ((i (1+ n)))
+ (while (<= i N)
+ (math-nlfit-set-elt r i (math-nlfit-get-elt d i))
+ (setq i (1+ i))))
+ (let ((j n))
+ (while (>= j 1)
+ (let ((i N))
+ (while (>= i (1+ j))
+ (setq rho (math-nlfit-get-matx-elt C i j))
+ (if (math-equal rho 1)
+ (setq gamma 0
+ sigma 1)
+ (if (math-lessp (calcFunc-abs rho) 1)
+ (setq sigma rho
+ gamma (calcFunc-sqrt
+ (math-sub 1 (math-mul sigma sigma))))
+ (setq gamma (math-div 1 (calcFunc-abs rho))
+ sigma (math-mul (calcFunc-sign rho)
+ (calcFunc-sqrt
+ (math-sub 1 (math-mul gamma
gamma)))))))
+ (let ((ri (math-nlfit-get-elt r i))
+ (rj (math-nlfit-get-elt r j)))
+ (setq h (math-add (math-mul gamma rj)
+ (math-mul sigma ri)))
+ (setq ri (math-sub
+ (math-mul gamma ri)
+ (math-mul sigma rj)))
+ (setq rj h)
+ (math-nlfit-set-elt r i ri)
+ (math-nlfit-set-elt r j rj))
+ (setq i (1- i))))
+ (setq j (1- j))))
+
+ x))
+
+(defun math-nlfit-jacobian (grad xlist parms &optional slist)
+ (let ((j nil))
+ (while xlist
+ (let ((row (apply grad (car xlist) parms)))
+ (setq j
+ (cons
+ (if slist
+ (mapcar (lambda (x) (math-div x (car slist))) row)
+ row)
+ j)))
+ (setq slist (cdr slist))
+ (setq xlist (cdr xlist)))
+ (reverse j)))
+
+(defun math-nlfit-make-ident (l n)
+ (let ((m (math-nlfit-make-matrix n n))
+ (i 1))
+ (while (<= i n)
+ (math-nlfit-set-matx-elt m i i l)
+ (setq i (1+ i)))
+ m))
+
+(defun math-nlfit-chi-sq (xlist ylist parms fn &optional slist)
+ (let ((cs 0))
+ (while xlist
+ (let ((c
+ (math-sub
+ (apply fn (car xlist) parms)
+ (car ylist))))
+ (if slist
+ (setq c (math-div c (car slist))))
+ (setq cs
+ (math-add cs
+ (math-mul c c))))
+ (setq xlist (cdr xlist))
+ (setq ylist (cdr ylist))
+ (setq slist (cdr slist)))
+ cs))
+
+(defun math-nlfit-init-lambda (C)
+ (let ((l 0)
+ (n (length (car C)))
+ (N (length C)))
+ (while C
+ (let ((row (car C)))
+ (while row
+ (setq l (math-add l (math-mul (car row) (car row))))
+ (setq row (cdr row))))
+ (setq C (cdr C)))
+ (calcFunc-sqrt (math-div l (math-mul n N)))))
+
+(defun math-nlfit-make-Ctilda (C l)
+ (let* ((n (length (car C)))
+ (bot (math-nlfit-make-ident l n)))
+ (append C bot)))
+
+(defun math-nlfit-make-d (fn xdata ydata parms &optional sdata)
+ (let ((d nil))
+ (while xdata
+ (setq d (cons
+ (let ((dd (math-sub (apply fn (car xdata) parms)
+ (car ydata))))
+ (if sdata (math-div dd (car sdata)) dd))
+ d))
+ (setq xdata (cdr xdata))
+ (setq ydata (cdr ydata))
+ (setq sdata (cdr sdata)))
+ (reverse d)))
+
+(defun math-nlfit-make-dtilda (d n)
+ (append d (make-list n 0)))
+
+(defun math-nlfit-fit (xlist ylist parms fn grad &optional slist)
+ (let*
+ ((C (math-nlfit-jacobian grad xlist parms slist))
+ (d (math-nlfit-make-d fn xlist ylist parms slist))
+ (chisq (math-nlfit-chi-sq xlist ylist parms fn slist))
+ (lambda (math-nlfit-init-lambda C))
+ (really-done nil)
+ (iters 0))
+ (while (and
+ (not really-done)
+ (< iters math-nlfit-max-its))
+ (setq iters (1+ iters))
+ (let ((done nil))
+ (while (not done)
+ (let* ((Ctilda (math-nlfit-make-Ctilda C lambda))
+ (dtilda (math-nlfit-make-dtilda d (length (car C))))
+ (zeta (math-nlfit-givens Ctilda dtilda))
+ (newparms (mapcar* 'math-add (copy-tree parms) zeta))
+ (newchisq (math-nlfit-chi-sq xlist ylist newparms fn slist)))
+ (if (math-lessp newchisq chisq)
+ (progn
+ (if (math-lessp
+ (math-div
+ (math-sub chisq newchisq) newchisq) math-nlfit-epsilon)
+ (setq really-done t))
+ (setq lambda (math-div lambda 10))
+ (setq chisq newchisq)
+ (setq parms newparms)
+ (setq done t))
+ (setq lambda (math-mul lambda 10)))))
+ (setq C (math-nlfit-jacobian grad xlist parms slist))
+ (setq d (math-nlfit-make-d fn xlist ylist parms slist))))
+ (list chisq parms)))
+
+;;; The functions that describe our models, and their gradients.
+
+(defun math-nlfit-s-logistic-fn (x a b c)
+ (math-div a (math-add 1 (math-mul b (calcFunc-exp (math-mul c x))))))
+
+(defun math-nlfit-s-logistic-grad (x a b c)
+ (let* ((ep (calcFunc-exp (math-mul c x)))
+ (d (math-add 1 (math-mul b ep)))
+ (d2 (math-mul d d)))
+ (list
+ (math-div 1 d)
+ (math-neg (math-div (math-mul a ep) d2))
+ (math-neg (math-div (math-mul a (math-mul b (math-mul x ep))) d2)))))
+
+(defun math-nlfit-b-logistic-fn (x a c d)
+ (let ((ex (calcFunc-exp (math-mul c (math-sub x d)))))
+ (math-div
+ (math-mul a ex)
+ (math-sqr
+ (math-add
+ 1 ex)))))
+
+(defun math-nlfit-b-logistic-grad (x a c d)
+ (let* ((ex (calcFunc-exp (math-mul c (math-sub x d))))
+ (ex1 (math-add 1 ex))
+ (xd (math-sub x d)))
+ (list
+ (math-div
+ ex
+ (math-sqr ex1))
+ (math-sub
+ (math-div
+ (math-mul a (math-mul xd ex))
+ (math-sqr ex1))
+ (math-div
+ (math-mul 2 (math-mul a (math-mul xd (math-sqr ex))))
+ (math-pow ex1 3)))
+ (math-sub
+ (math-div
+ (math-mul 2 (math-mul a (math-mul c (math-sqr ex))))
+ (math-pow ex1 3))
+ (math-div
+ (math-mul a (math-mul c ex))
+ (math-sqr ex1))))))
+
+;;; Functions to get the final covariance matrix and the sdevs
+
+(defun math-nlfit-find-covar (grad xlist pparms)
+ (let ((j nil))
+ (while xlist
+ (setq j (cons (cons 'vec (apply grad (car xlist) pparms)) j))
+ (setq xlist (cdr xlist)))
+ (setq j (cons 'vec (reverse j)))
+ (setq j
+ (math-mul
+ (calcFunc-trn j) j))
+ (calcFunc-inv j)))
+
+(defun math-nlfit-get-sigmas (grad xlist pparms chisq)
+ (let* ((sgs nil)
+ (covar (math-nlfit-find-covar grad xlist pparms))
+ (n (1- (length covar)))
+ (N (length xlist))
+ (i 1))
+ (when (> N n)
+ (while (<= i n)
+ (setq sgs (cons (calcFunc-sqrt (nth i (nth i covar))) sgs))
+ (setq i (1+ i)))
+ (setq sgs (reverse sgs)))
+ (list sgs covar)))
+
+;;; Now the Calc functions
+
+(defun math-nlfit-s-logistic-params (xdata ydata)
+ (let ((pdata (math-nlfit-get-rates-from-cumul xdata ydata)))
+ (math-nlfit-find-logistic-parameters ydata pdata xdata)))
+
+(defun math-nlfit-b-logistic-params (xdata ydata)
+ (let* ((q0 (math-nlfit-find-q0 ydata xdata))
+ (qdata (math-nlfit-get-cumul-from-rates xdata ydata q0))
+ (abc (math-nlfit-find-logistic-parameters qdata ydata xdata))
+ (B (nth 1 abc))
+ (C (nth 2 abc))
+ (A (math-neg
+ (math-mul
+ (nth 0 abc)
+ (math-mul B C))))
+ (D (math-neg (math-div (calcFunc-ln B) C)))
+ (A (math-div A B)))
+ (list A C D)))
+
+;;; Some functions to turn the parameter lists and variables
+;;; into the appropriate functions.
+
+(defun math-nlfit-s-logistic-solnexpr (pms var)
+ (let ((a (nth 0 pms))
+ (b (nth 1 pms))
+ (c (nth 2 pms)))
+ (list '/ a
+ (list '+
+ 1
+ (list '*
+ b
+ (calcFunc-exp
+ (list '*
+ c
+ var)))))))
+
+(defun math-nlfit-b-logistic-solnexpr (pms var)
+ (let ((a (nth 0 pms))
+ (c (nth 1 pms))
+ (d (nth 2 pms)))
+ (list '/
+ (list '*
+ a
+ (calcFunc-exp
+ (list '*
+ c
+ (list '- var d))))
+ (list '^
+ (list '+
+ 1
+ (calcFunc-exp
+ (list '*
+ c
+ (list '- var d))))
+ 2))))
+
+(defun math-nlfit-enter-result (n prefix vals)
+ (setq calc-aborted-prefix prefix)
+ (calc-pop-push-record-list n prefix vals)
+ (calc-handle-whys))
+
+(defun math-nlfit-fit-curve (fn grad solnexpr initparms &optional sdv)
+ (calc-slow-wrapper
+ (let* ((sdevv (or (eq sdv 'calcFunc-efit) (eq sdv 'calcFunc-xfit)))
+ (calc-display-working-message nil)
+ (data (calc-top 1))
+ (xdata (cdr (car (cdr data))))
+ (ydata (cdr (car (cdr (cdr data)))))
+ (sdata (if (math-contains-sdev-p ydata)
+ (mapcar (lambda (x) (math-get-sdev x t)) ydata)
+ nil))
+ (ydata (mapcar (lambda (x) (math-get-value x)) ydata))
+ (zzz (progn (setq j1 xdata j2 ydata j3 sdata) 1))
+
+ (calc-curve-varnames nil)
+ (calc-curve-coefnames nil)
+ (calc-curve-nvars 1)
+ (fitvars (calc-get-fit-variables 1 3))
+ (var (nth 1 calc-curve-varnames))
+ (parms (cdr calc-curve-coefnames))
+ (parmguess
+ (funcall initparms xdata ydata))
+ (fit (math-nlfit-fit xdata ydata parmguess fn grad sdata))
+ (finalparms (nth 1 fit))
+ (sigmacovar
+ (if sdevv
+ (math-nlfit-get-sigmas grad xdata finalparms (nth 0 fit))))
+ (sigmas
+ (if sdevv
+ (nth 0 sigmacovar)))
+ (finalparms
+ (if sigmas
+ (mapcar* (lambda (x y) (list 'sdev x y)) finalparms sigmas)
+ finalparms))
+ (soln (funcall solnexpr finalparms var)))
+ (let ((calc-fit-to-trail t)
+ (traillist nil))
+ (while parms
+ (setq traillist (cons (list 'calcFunc-eq (car parms) (car finalparms))
+ traillist))
+ (setq finalparms (cdr finalparms))
+ (setq parms (cdr parms)))
+ (setq traillist (calc-normalize (cons 'vec (nreverse traillist))))
+ (cond ((eq sdv 'calcFunc-efit)
+ (math-nlfit-enter-result 1 "efit" soln))
+ ((eq sdv 'calcFunc-xfit)
+ (let (sln)
+ (setq sln
+ (list 'vec
+ soln
+ traillist
+ (nth 1 sigmacovar)
+ '(vec)
+ (nth 0 fit)
+ (let ((n (length xdata))
+ (m (length finalparms)))
+ (if (and sdata (> n m))
+ (calcFunc-utpc (nth 0 fit)
+ (- n m))
+ '(var nan var-nan)))))
+ (math-nlfit-enter-result 1 "xfit" sln)))
+ (t
+ (math-nlfit-enter-result 1 "fit" soln)))
+ (calc-record traillist "parm")))))
+
+(defun calc-fit-s-shaped-logistic-curve (arg)
+ (interactive "P")
+ (math-nlfit-fit-curve 'math-nlfit-s-logistic-fn
+ 'math-nlfit-s-logistic-grad
+ 'math-nlfit-s-logistic-solnexpr
+ 'math-nlfit-s-logistic-params
+ arg))
+
+(defun calc-fit-bell-shaped-logistic-curve (arg)
+ (interactive "P")
+ (math-nlfit-fit-curve 'math-nlfit-b-logistic-fn
+ 'math-nlfit-b-logistic-grad
+ 'math-nlfit-b-logistic-solnexpr
+ 'math-nlfit-b-logistic-params
+ arg))
+
+(defun calc-fit-hubbert-linear-curve (&optional sdv)
+ (calc-slow-wrapper
+ (let* ((sdevv (or (eq sdv 'calcFunc-efit) (eq sdv 'calcFunc-xfit)))
+ (calc-display-working-message nil)
+ (data (calc-top 1))
+ (qdata (cdr (car (cdr data))))
+ (pdata (cdr (car (cdr (cdr data)))))
+ (sdata (if (math-contains-sdev-p pdata)
+ (mapcar (lambda (x) (math-get-sdev x t)) pdata)
+ nil))
+ (pdata (mapcar (lambda (x) (math-get-value x)) pdata))
+ (poverqdata (mapcar* 'math-div pdata qdata))
+ (parmvals (math-nlfit-least-squares qdata poverqdata sdata sdevv))
+ (finalparms (list (nth 0 parmvals)
+ (math-neg
+ (math-div (nth 0 parmvals)
+ (nth 1 parmvals)))))
+ (calc-curve-varnames nil)
+ (calc-curve-coefnames nil)
+ (calc-curve-nvars 1)
+ (fitvars (calc-get-fit-variables 1 2))
+ (var (nth 1 calc-curve-varnames))
+ (parms (cdr calc-curve-coefnames))
+ (soln (list '* (nth 0 finalparms)
+ (list '- 1
+ (list '/ var (nth 1 finalparms))))))
+ (let ((calc-fit-to-trail t)
+ (traillist nil))
+ (setq traillist
+ (list 'vec
+ (list 'calcFunc-eq (nth 0 parms) (nth 0 finalparms))
+ (list 'calcFunc-eq (nth 1 parms) (nth 1 finalparms))))
+ (cond ((eq sdv 'calcFunc-efit)
+ (math-nlfit-enter-result 1 "efit" soln))
+ ((eq sdv 'calcFunc-xfit)
+ (let (sln
+ (chisq
+ (math-nlfit-chi-sq
+ qdata poverqdata
+ (list (nth 1 (nth 0 finalparms))
+ (nth 1 (nth 1 finalparms)))
+ (lambda (x a b)
+ (math-mul a
+ (math-sub
+ 1
+ (math-div x b))))
+ sdata)))
+ (setq sln
+ (list 'vec
+ soln
+ traillist
+ (nth 2 parmvals)
+ (list
+ 'vec
+ '(calcFunc-fitdummy 1)
+ (list 'calcFunc-neg
+ (list '/
+ '(calcFunc-fitdummy 1)
+ '(calcFunc-fitdummy 2))))
+ chisq
+ (let ((n (length qdata)))
+ (if (and sdata (> n 2))
+ (calcFunc-utpc
+ chisq
+ (- n 2))
+ '(var nan var-nan)))))
+ (math-nlfit-enter-result 1 "xfit" sln)))
+ (t
+ (math-nlfit-enter-result 1 "fit" soln)))
+ (calc-record traillist "parm")))))
+
+(provide 'calc-nlfit)
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Emacs-diffs] Changes to calc-nlfit.el,
Jay Belanger <=