annotate lisp/calc/calc-nlfit.el @ 86075:a7a49982b6d9

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