aboutsummaryrefslogtreecommitdiffstats
path: root/lisp/calc
diff options
context:
space:
mode:
authorPaul Eggert2011-11-22 23:03:56 -0800
committerPaul Eggert2011-11-22 23:03:56 -0800
commitee7683ebb70c308e596103e379ef6b91d001eebc (patch)
tree45382191619e27df7106a8d990bc03903d0fa03f /lisp/calc
parent6b21de180fba10432988d94d2b8f3e2521be5b17 (diff)
downloademacs-ee7683ebb70c308e596103e379ef6b91d001eebc.tar.gz
emacs-ee7683ebb70c308e596103e379ef6b91d001eebc.zip
Spelling fixes.
Diffstat (limited to 'lisp/calc')
-rw-r--r--lisp/calc/calc-nlfit.el91
-rw-r--r--lisp/calc/calc.el2
2 files changed, 46 insertions, 47 deletions
diff --git a/lisp/calc/calc-nlfit.el b/lisp/calc/calc-nlfit.el
index 37e6f66c1b1..bd162866c31 100644
--- a/lisp/calc/calc-nlfit.el
+++ b/lisp/calc/calc-nlfit.el
@@ -22,7 +22,7 @@
22;;; Commentary: 22;;; Commentary:
23 23
24;; This code uses the Levenberg-Marquardt method, as described in 24;; This code uses the Levenberg-Marquardt method, as described in
25;; _Numerical Analysis_ by H. R. Schwarz, to fit data to 25;; _Numerical Analysis_ by H. R. Schwarz, to fit data to
26;; nonlinear curves. Currently, the only the following curves are 26;; nonlinear curves. Currently, the only the following curves are
27;; supported: 27;; supported:
28;; The logistic S curve, y=a/(1+exp(b*(t-c))) 28;; The logistic S curve, y=a/(1+exp(b*(t-c)))
@@ -33,14 +33,14 @@
33 33
34;; The logistic bell curve, y=A*exp(B*(t-C))/(1+exp(B*(t-C)))^2 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. 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 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 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 38;; data as consisting of rates p0, p1, ..., pn and their
39;; respective times t0, t1, ..., tn. 39;; respective times t0, t1, ..., tn.
40 40
41;; The Hubbert Linearization, y/x=A*(1-x/B) 41;; The Hubbert Linearization, y/x=A*(1-x/B)
42;; Here, y is thought of as the rate of growth of a population 42;; Here, y is thought of as the rate of growth of a population
43;; and x represents the actual population. This is essentially 43;; and x represents the actual population. This is essentially
44;; the differential equation describing the actual population. 44;; the differential equation describing the actual population.
45 45
46;; The Levenberg-Marquardt method is an iterative process: it takes 46;; The Levenberg-Marquardt method is an iterative process: it takes
@@ -53,7 +53,7 @@
53;; approximations for b and c are found using least squares on the 53;; approximations for b and c are found using least squares on the
54;; linearization log((a/y)-1) = log(bb) + cc*t of 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 55;; y=a/(1+bb*exp(cc*t)), which is equivalent to the above s curve
56;; formula, and then tranlating it to b and c. From this, we can 56;; formula, and then translating it to b and c. From this, we can
57;; also get approximations for the bell curve parameters. 57;; also get approximations for the bell curve parameters.
58 58
59;;; Code: 59;;; Code:
@@ -68,7 +68,7 @@
68(defun math-nlfit-least-squares (xdata ydata &optional sdata sigmas) 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." 69 "Return the parameters A and B for the best least squares fit y=a+bx."
70 (let* ((n (length xdata)) 70 (let* ((n (length xdata))
71 (s2data (if sdata 71 (s2data (if sdata
72 (mapcar 'calcFunc-sqr sdata) 72 (mapcar 'calcFunc-sqr sdata)
73 (make-list n 1))) 73 (make-list n 1)))
74 (S (if sdata 0 n)) 74 (S (if sdata 0 n))
@@ -109,11 +109,11 @@
109;;; The methods described by de Sousa require the cumulative data qdata 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 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 111;;; qdata and the corresponding times tdata, or pdata and the corresponding
112;;; tdata. The following two functions will find pdata or qdata, 112;;; tdata. The following two functions will find pdata or qdata,
113;;; given the other.. 113;;; given the other..
114 114
115;;; First, given two lists; one of values q0, q1, ..., qn and one of 115;;; First, given two lists; one of values q0, q1, ..., qn and one of
116;;; corresponding times t0, t1, ..., tn; return a list 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. 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). 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)). 119;;; pn is the left hand derivative (qn - q(n-1))/(tn - t(n-1)).
@@ -122,7 +122,7 @@
122 122
123(defun math-nlfit-get-rates-from-cumul (tdata qdata) 123(defun math-nlfit-get-rates-from-cumul (tdata qdata)
124 (let ((pdata (list 124 (let ((pdata (list
125 (math-div 125 (math-div
126 (math-sub (nth 1 qdata) 126 (math-sub (nth 1 qdata)
127 (nth 0 qdata)) 127 (nth 0 qdata))
128 (math-sub (nth 1 tdata) 128 (math-sub (nth 1 tdata)
@@ -155,7 +155,7 @@
155 pdata)) 155 pdata))
156 (reverse pdata))) 156 (reverse pdata)))
157 157
158;;; Next, given two lists -- one of rates p0, p1, ..., pn and one of 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, 159;;; corresponding times t0, t1, ..., tn -- and an initial values q0,
160;;; return a list q0, q1, ..., qn of the cumulative values. 160;;; return a list q0, q1, ..., qn of the cumulative values.
161;;; q0 is the initial value given. 161;;; q0 is the initial value given.
@@ -169,7 +169,7 @@
169 (cons 169 (cons
170 (math-add (car qdata) 170 (math-add (car qdata)
171 (math-mul 171 (math-mul
172 (math-mul 172 (math-mul
173 '(float 5 -1) 173 '(float 5 -1)
174 (math-add (nth 1 pdata) (nth 0 pdata))) 174 (math-add (nth 1 pdata) (nth 0 pdata)))
175 (math-sub (nth 1 tdata) 175 (math-sub (nth 1 tdata)
@@ -181,13 +181,13 @@
181 181
182;;; Given the qdata, pdata and tdata, find the parameters 182;;; Given the qdata, pdata and tdata, find the parameters
183;;; a, b and c that fit q = a/(1+b*exp(c*t)). 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. 184;;; a is found using the method described by de Sousa.
185;;; b and c are found using least squares on the linearization 185;;; b and c are found using least squares on the linearization
186;;; log((a/q)-1) = log(b) + c*t 186;;; log((a/q)-1) = log(b) + c*t
187;;; In some cases (where the logistic curve may well be the wrong 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 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. 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 190;;; In this case, a will be replaced by a number slightly above
191;;; the maximum value of q. 191;;; the maximum value of q.
192 192
193(defun math-nlfit-find-qmax (qdata pdata tdata) 193(defun math-nlfit-find-qmax (qdata pdata tdata)
@@ -224,7 +224,7 @@
224 (setq qmh 224 (setq qmh
225 (math-add qmh 225 (math-add qmh
226 (math-mul 226 (math-mul
227 (math-mul 227 (math-mul
228 '(float 5 -1) 228 '(float 5 -1)
229 (math-add (nth 1 pdata) (nth 0 pdata))) 229 (math-add (nth 1 pdata) (nth 0 pdata)))
230 (math-sub (nth 1 tdata) 230 (math-sub (nth 1 tdata)
@@ -239,7 +239,7 @@
239 (let* ((qhalf (math-nlfit-find-qmaxhalf pdata tdata)) 239 (let* ((qhalf (math-nlfit-find-qmaxhalf pdata tdata))
240 (q0 (math-mul 2 qhalf)) 240 (q0 (math-mul 2 qhalf))
241 (qdata (math-nlfit-get-cumul-from-rates tdata pdata q0))) 241 (qdata (math-nlfit-get-cumul-from-rates tdata pdata q0)))
242 (while (math-lessp (math-nlfit-find-qmax 242 (while (math-lessp (math-nlfit-find-qmax
243 (mapcar 243 (mapcar
244 (lambda (q) (math-add q0 q)) 244 (lambda (q) (math-add q0 q))
245 qdata) 245 qdata)
@@ -260,7 +260,7 @@
260 (i 0)) 260 (i 0))
261 (while (< i 10) 261 (while (< i 10)
262 (setq q0 (math-mul '(float 5 -1) (math-add qmin qmax))) 262 (setq q0 (math-mul '(float 5 -1) (math-add qmin qmax)))
263 (if (math-lessp 263 (if (math-lessp
264 (math-nlfit-find-qmax 264 (math-nlfit-find-qmax
265 (mapcar 265 (mapcar
266 (lambda (q) (math-add q0 q)) 266 (lambda (q) (math-add q0 q))
@@ -272,7 +272,7 @@
272 (setq i (1+ i))) 272 (setq i (1+ i)))
273 (math-mul '(float 5 -1) (math-add qmin qmax))))) 273 (math-mul '(float 5 -1) (math-add qmin qmax)))))
274 274
275;;; To improve the approximations to the parameters, we can use 275;;; To improve the approximations to the parameters, we can use
276;;; Marquardt method as described in Schwarz's book. 276;;; Marquardt method as described in Schwarz's book.
277 277
278;;; Small numbers used in the Givens algorithm 278;;; Small numbers used in the Givens algorithm
@@ -329,7 +329,7 @@
329 (let ((cij (math-nlfit-get-matx-elt C i j)) 329 (let ((cij (math-nlfit-get-matx-elt C i j))
330 (cjj (math-nlfit-get-matx-elt C j j))) 330 (cjj (math-nlfit-get-matx-elt C j j)))
331 (when (not (math-equal 0 cij)) 331 (when (not (math-equal 0 cij))
332 (if (math-lessp (calcFunc-abs cjj) 332 (if (math-lessp (calcFunc-abs cjj)
333 (math-mul math-nlfit-delta (calcFunc-abs cij))) 333 (math-mul math-nlfit-delta (calcFunc-abs cij)))
334 (setq w (math-neg cij) 334 (setq w (math-neg cij)
335 gamma 0 335 gamma 0
@@ -337,7 +337,7 @@
337 rho 1) 337 rho 1)
338 (setq w (math-mul 338 (setq w (math-mul
339 (calcFunc-sign cjj) 339 (calcFunc-sign cjj)
340 (calcFunc-sqrt 340 (calcFunc-sqrt
341 (math-add 341 (math-add
342 (math-mul cjj cjj) 342 (math-mul cjj cjj)
343 (math-mul cij cij)))) 343 (math-mul cij cij))))
@@ -351,10 +351,10 @@
351 (math-nlfit-set-matx-elt C j j w) 351 (math-nlfit-set-matx-elt C j j w)
352 (math-nlfit-set-matx-elt C i j rho) 352 (math-nlfit-set-matx-elt C i j rho)
353 (let ((k (1+ j))) 353 (let ((k (1+ j)))
354 (while (<= k n) 354 (while (<= k n)
355 (let* ((cjk (math-nlfit-get-matx-elt C j k)) 355 (let* ((cjk (math-nlfit-get-matx-elt C j k))
356 (cik (math-nlfit-get-matx-elt C i k)) 356 (cik (math-nlfit-get-matx-elt C i k))
357 (h (math-sub 357 (h (math-sub
358 (math-mul gamma cjk) (math-mul sigma cik)))) 358 (math-mul gamma cjk) (math-mul sigma cik))))
359 (setq cik (math-add 359 (setq cik (math-add
360 (math-mul sigma cjk) 360 (math-mul sigma cjk)
@@ -386,9 +386,9 @@
386 (setq s (math-add s (math-mul (math-nlfit-get-matx-elt C i k) 386 (setq s (math-add s (math-mul (math-nlfit-get-matx-elt C i k)
387 (math-nlfit-get-elt x k)))) 387 (math-nlfit-get-elt x k))))
388 (setq k (1+ k)))) 388 (setq k (1+ k))))
389 (math-nlfit-set-elt x i 389 (math-nlfit-set-elt x i
390 (math-neg 390 (math-neg
391 (math-div s 391 (math-div s
392 (math-nlfit-get-matx-elt C i i)))) 392 (math-nlfit-get-matx-elt C i i))))
393 (setq i (1- i)))) 393 (setq i (1- i))))
394 (let ((i (1+ n))) 394 (let ((i (1+ n)))
@@ -405,7 +405,7 @@
405 sigma 1) 405 sigma 1)
406 (if (math-lessp (calcFunc-abs rho) 1) 406 (if (math-lessp (calcFunc-abs rho) 1)
407 (setq sigma rho 407 (setq sigma rho
408 gamma (calcFunc-sqrt 408 gamma (calcFunc-sqrt
409 (math-sub 1 (math-mul sigma sigma)))) 409 (math-sub 1 (math-mul sigma sigma))))
410 (setq gamma (math-div 1 (calcFunc-abs rho)) 410 (setq gamma (math-div 1 (calcFunc-abs rho))
411 sigma (math-mul (calcFunc-sign rho) 411 sigma (math-mul (calcFunc-sign rho)
@@ -429,7 +429,7 @@
429 429
430(defun math-nlfit-jacobian (grad xlist parms &optional slist) 430(defun math-nlfit-jacobian (grad xlist parms &optional slist)
431 (let ((j nil)) 431 (let ((j nil))
432 (while xlist 432 (while xlist
433 (let ((row (apply grad (car xlist) parms))) 433 (let ((row (apply grad (car xlist) parms)))
434 (setq j 434 (setq j
435 (cons 435 (cons
@@ -495,7 +495,7 @@
495 (setq ydata (cdr ydata)) 495 (setq ydata (cdr ydata))
496 (setq sdata (cdr sdata))) 496 (setq sdata (cdr sdata)))
497 (reverse d))) 497 (reverse d)))
498 498
499(defun math-nlfit-make-dtilda (d n) 499(defun math-nlfit-make-dtilda (d n)
500 (append d (make-list n 0))) 500 (append d (make-list n 0)))
501 501
@@ -520,8 +520,8 @@
520 (newchisq (math-nlfit-chi-sq xlist ylist newparms fn slist))) 520 (newchisq (math-nlfit-chi-sq xlist ylist newparms fn slist)))
521 (if (math-lessp newchisq chisq) 521 (if (math-lessp newchisq chisq)
522 (progn 522 (progn
523 (if (math-lessp 523 (if (math-lessp
524 (math-div 524 (math-div
525 (math-sub chisq newchisq) newchisq) math-nlfit-epsilon) 525 (math-sub chisq newchisq) newchisq) math-nlfit-epsilon)
526 (setq really-done t)) 526 (setq really-done t))
527 (setq lambda (math-div lambda 10)) 527 (setq lambda (math-div lambda 10))
@@ -551,7 +551,7 @@
551 (let ((ex (calcFunc-exp (math-mul c (math-sub x d))))) 551 (let ((ex (calcFunc-exp (math-mul c (math-sub x d)))))
552 (math-div 552 (math-div
553 (math-mul a ex) 553 (math-mul a ex)
554 (math-sqr 554 (math-sqr
555 (math-add 555 (math-add
556 1 ex))))) 556 1 ex)))))
557 557
@@ -582,7 +582,7 @@
582 582
583(defun math-nlfit-find-covar (grad xlist pparms) 583(defun math-nlfit-find-covar (grad xlist pparms)
584 (let ((j nil)) 584 (let ((j nil))
585 (while xlist 585 (while xlist
586 (setq j (cons (cons 'vec (apply grad (car xlist) pparms)) j)) 586 (setq j (cons (cons 'vec (apply grad (car xlist) pparms)) j))
587 (setq xlist (cdr xlist))) 587 (setq xlist (cdr xlist)))
588 (setq j (cons 'vec (reverse j))) 588 (setq j (cons 'vec (reverse j)))
@@ -603,7 +603,7 @@
603 (setq i (1+ i))) 603 (setq i (1+ i)))
604 (setq sgs (reverse sgs))) 604 (setq sgs (reverse sgs)))
605 (list sgs covar))) 605 (list sgs covar)))
606 606
607;;; Now the Calc functions 607;;; Now the Calc functions
608 608
609(defun math-nlfit-s-logistic-params (xdata ydata) 609(defun math-nlfit-s-logistic-params (xdata ydata)
@@ -687,15 +687,15 @@
687 (funcall initparms xdata ydata)) 687 (funcall initparms xdata ydata))
688 (fit (math-nlfit-fit xdata ydata parmguess fn grad sdata)) 688 (fit (math-nlfit-fit xdata ydata parmguess fn grad sdata))
689 (finalparms (nth 1 fit)) 689 (finalparms (nth 1 fit))
690 (sigmacovar 690 (sigmacovar
691 (if sdevv 691 (if sdevv
692 (math-nlfit-get-sigmas grad xdata finalparms (nth 0 fit)))) 692 (math-nlfit-get-sigmas grad xdata finalparms (nth 0 fit))))
693 (sigmas 693 (sigmas
694 (if sdevv 694 (if sdevv
695 (nth 0 sigmacovar))) 695 (nth 0 sigmacovar)))
696 (finalparms 696 (finalparms
697 (if sigmas 697 (if sigmas
698 (math-map-binop 698 (math-map-binop
699 (lambda (x y) (list 'sdev x y)) finalparms sigmas) 699 (lambda (x y) (list 'sdev x y)) finalparms sigmas)
700 finalparms)) 700 finalparms))
701 (soln (funcall solnexpr finalparms var))) 701 (soln (funcall solnexpr finalparms var)))
@@ -712,8 +712,8 @@
712 ((eq sdv 'calcFunc-xfit) 712 ((eq sdv 'calcFunc-xfit)
713 (let (sln) 713 (let (sln)
714 (setq sln 714 (setq sln
715 (list 'vec 715 (list 'vec
716 soln 716 soln
717 traillist 717 traillist
718 (nth 1 sigmacovar) 718 (nth 1 sigmacovar)
719 '(vec) 719 '(vec)
@@ -721,7 +721,7 @@
721 (let ((n (length xdata)) 721 (let ((n (length xdata))
722 (m (length finalparms))) 722 (m (length finalparms)))
723 (if (and sdata (> n m)) 723 (if (and sdata (> n m))
724 (calcFunc-utpc (nth 0 fit) 724 (calcFunc-utpc (nth 0 fit)
725 (- n m)) 725 (- n m))
726 '(var nan var-nan))))) 726 '(var nan var-nan)))))
727 (math-nlfit-enter-result 1 "xfit" sln))) 727 (math-nlfit-enter-result 1 "xfit" sln)))
@@ -787,14 +787,14 @@
787 (list (nth 1 (nth 0 finalparms)) 787 (list (nth 1 (nth 0 finalparms))
788 (nth 1 (nth 1 finalparms))) 788 (nth 1 (nth 1 finalparms)))
789 (lambda (x a b) 789 (lambda (x a b)
790 (math-mul a 790 (math-mul a
791 (math-sub 791 (math-sub
792 1 792 1
793 (math-div x b)))) 793 (math-div x b))))
794 sdata))) 794 sdata)))
795 (setq sln 795 (setq sln
796 (list 'vec 796 (list 'vec
797 soln 797 soln
798 traillist 798 traillist
799 (nth 2 parmvals) 799 (nth 2 parmvals)
800 (list 800 (list
@@ -807,7 +807,7 @@
807 chisq 807 chisq
808 (let ((n (length qdata))) 808 (let ((n (length qdata)))
809 (if (and sdata (> n 2)) 809 (if (and sdata (> n 2))
810 (calcFunc-utpc 810 (calcFunc-utpc
811 chisq 811 chisq
812 (- n 2)) 812 (- n 2))
813 '(var nan var-nan))))) 813 '(var nan var-nan)))))
@@ -817,4 +817,3 @@
817 (calc-record traillist "parm"))))) 817 (calc-record traillist "parm")))))
818 818
819(provide 'calc-nlfit) 819(provide 'calc-nlfit)
820
diff --git a/lisp/calc/calc.el b/lisp/calc/calc.el
index 626d2462b4f..23f955afe7c 100644
--- a/lisp/calc/calc.el
+++ b/lisp/calc/calc.el
@@ -1003,7 +1003,7 @@ Used by `calc-user-invocation'.")
1003(defvar calc-quick-prev-results nil 1003(defvar calc-quick-prev-results nil
1004 "Previous results from Quick Calc.") 1004 "Previous results from Quick Calc.")
1005(defvar calc-said-hello nil 1005(defvar calc-said-hello nil
1006 "Non-nil if the welcomd message has been displayed.") 1006 "Non-nil if the welcome message has been displayed.")
1007(defvar calc-executing-macro nil 1007(defvar calc-executing-macro nil
1008 "Non-nil if a keyboard macro is executing from the \"K\" key.") 1008 "Non-nil if a keyboard macro is executing from the \"K\" key.")
1009(defvar calc-any-selections nil 1009(defvar calc-any-selections nil