annotate lisp/calc/calc-nlfit.el @ 87023:5405672da978

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