1 ;;; calc-nlfit.el --- nonlinear curve fitting for Calc
3 ;; Copyright (C) 2007, 2008 Free Software Foundation, Inc.
5 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
7 ;; This file is part of GNU Emacs.
9 ;; GNU Emacs is free software; you can redistribute it and/or modify
10 ;; it under the terms of the GNU General Public License as published by
11 ;; the Free Software Foundation; either version 3, or (at your option)
14 ;; GNU Emacs is distributed in the hope that it will be useful,
15 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 ;; GNU General Public License for more details.
19 ;; You should have received a copy of the GNU General Public License
20 ;; along with GNU Emacs; see the file COPYING. If not, write to the
21 ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 ;; Boston, MA 02110-1301, USA.
26 ;; This code uses the Levenberg-Marquardt method, as described in
27 ;; _Numerical Analysis_ by H. R. Schwarz, to fit data to
28 ;; nonlinear curves. Currently, the only the following curves are
30 ;; The logistic S curve, y=a/(1+exp(b*(t-c)))
31 ;; Here, y is usually interpreted as the population of some
32 ;; quantity at time t. So we will think of the data as consisting
33 ;; of quantities q0, q1, ..., qn and their respective times
36 ;; The logistic bell curve, y=A*exp(B*(t-C))/(1+exp(B*(t-C)))^2
37 ;; Note that this is the derivative of the formula for the S curve.
38 ;; We get A=-a*b, B=b and C=c. Here, y is interpreted as the rate
39 ;; of growth of a population at time t. So we will think of the
40 ;; data as consisting of rates p0, p1, ..., pn and their
41 ;; respective times t0, t1, ..., tn.
43 ;; The Hubbert Linearization, y/x=A*(1-x/B)
44 ;; Here, y is thought of as the rate of growth of a population
45 ;; and x represents the actual population. This is essentially
46 ;; the differential equation describing the actual population.
48 ;; The Levenberg-Marquardt method is an iterative process: it takes
49 ;; an initial guess for the parameters and refines them. To get an
50 ;; initial guess for the parameters, we'll use a method described by
51 ;; Luis de Sousa in "Hubbert's Peak Mathematics". The idea is that
52 ;; given quantities Q and the corresponding rates P, they should
53 ;; satisfy P/Q= mQ+a. We can use the parameter a for an
54 ;; approximation for the parameter a in the S curve, and
55 ;; approximations for b and c are found using least squares on the
56 ;; linearization log((a/y)-1) = log(bb) + cc*t of
57 ;; y=a/(1+bb*exp(cc*t)), which is equivalent to the above s curve
58 ;; formula, and then tranlating it to b and c. From this, we can
59 ;; also get approximations for the bell curve parameters.
66 ;; Declare functions which are defined elsewhere.
67 (declare-function calc-get-fit-variables
"calcalg3" (nv nc
&optional defv defc with-y homog
))
68 (declare-function math-map-binop
"calcalg3" (binop args1 args2
))
70 (defun math-nlfit-least-squares (xdata ydata
&optional sdata sigmas
)
71 "Return the parameters A and B for the best least squares fit y=a+bx."
72 (let* ((n (length xdata
))
74 (mapcar 'calcFunc-sqr sdata
)
86 (setq Sx
(math-add Sx
(if s
(math-div x s
) x
)))
87 (setq Sy
(math-add Sy
(if s
(math-div y s
) y
)))
88 (setq Sxx
(math-add Sxx
(if s
(math-div (math-mul x x
) s
)
90 (setq Sxy
(math-add Sxy
(if s
(math-div (math-mul x y
) s
)
93 (setq S
(math-add S
(math-div 1 s
)))))
94 (setq xdata
(cdr xdata
))
95 (setq ydata
(cdr ydata
))
96 (setq s2data
(cdr s2data
)))
97 (setq D
(math-sub (math-mul S Sxx
) (math-mul Sx Sx
)))
98 (let ((A (math-div (math-sub (math-mul Sxx Sy
) (math-mul Sx Sxy
)) D
))
99 (B (math-div (math-sub (math-mul S Sxy
) (math-mul Sx Sy
)) D
)))
101 (let ((C11 (math-div Sxx D
))
102 (C12 (math-neg (math-div Sx D
)))
103 (C22 (math-div S D
)))
104 (list (list 'sdev A
(calcFunc-sqrt C11
))
105 (list 'sdev B
(calcFunc-sqrt C22
))
108 (list 'vec C12 C22
))))
111 ;;; The methods described by de Sousa require the cumulative data qdata
112 ;;; and the rates pdata. We will assume that we are given either
113 ;;; qdata and the corresponding times tdata, or pdata and the corresponding
114 ;;; tdata. The following two functions will find pdata or qdata,
115 ;;; given the other..
117 ;;; First, given two lists; one of values q0, q1, ..., qn and one of
118 ;;; corresponding times t0, t1, ..., tn; return a list
119 ;;; p0, p1, ..., pn of the rates of change of the qi with respect to t.
120 ;;; p0 is the right hand derivative (q1 - q0)/(t1 - t0).
121 ;;; pn is the left hand derivative (qn - q(n-1))/(tn - t(n-1)).
122 ;;; The other pis are the averages of the two:
123 ;;; (1/2)((qi - q(i-1))/(ti - t(i-1)) + (q(i+1) - qi)/(t(i+1) - ti)).
125 (defun math-nlfit-get-rates-from-cumul (tdata qdata
)
128 (math-sub (nth 1 qdata
)
130 (math-sub (nth 1 tdata
)
132 (while (> (length qdata
) 2)
139 (math-sub (nth 2 qdata
)
141 (math-sub (nth 2 tdata
)
144 (math-sub (nth 1 qdata
)
146 (math-sub (nth 1 tdata
)
149 (setq qdata
(cdr qdata
)))
153 (math-sub (nth 1 qdata
)
155 (math-sub (nth 1 tdata
)
160 ;;; Next, given two lists -- one of rates p0, p1, ..., pn and one of
161 ;;; corresponding times t0, t1, ..., tn -- and an initial values q0,
162 ;;; return a list q0, q1, ..., qn of the cumulative values.
163 ;;; q0 is the initial value given.
164 ;;; For i>0, qi is computed using the trapezoid rule:
165 ;;; qi = q(i-1) + (1/2)(pi + p(i-1))(ti - t(i-1))
167 (defun math-nlfit-get-cumul-from-rates (tdata pdata q0
)
168 (let* ((qdata (list q0
)))
172 (math-add (car qdata
)
176 (math-add (nth 1 pdata
) (nth 0 pdata
)))
177 (math-sub (nth 1 tdata
)
180 (setq pdata
(cdr pdata
))
181 (setq tdata
(cdr tdata
)))
184 ;;; Given the qdata, pdata and tdata, find the parameters
185 ;;; a, b and c that fit q = a/(1+b*exp(c*t)).
186 ;;; a is found using the method described by de Sousa.
187 ;;; b and c are found using least squares on the linearization
188 ;;; log((a/q)-1) = log(b) + c*t
189 ;;; In some cases (where the logistic curve may well be the wrong
190 ;;; model), the computed a will be less than or equal to the maximum
191 ;;; value of q in qdata; in which case the above linearization won't work.
192 ;;; In this case, a will be replaced by a number slightly above
193 ;;; the maximum value of q.
195 (defun math-nlfit-find-qmax (qdata pdata tdata
)
196 (let* ((ratios (math-map-binop 'math-div pdata qdata
))
197 (lsdata (math-nlfit-least-squares ratios tdata
))
198 (qmax (math-max-list (car qdata
) (cdr qdata
)))
199 (a (math-neg (math-div (nth 1 lsdata
) (nth 0 lsdata
)))))
200 (if (math-lessp a qmax
)
201 (math-add '(float 5 -
1) qmax
)
204 (defun math-nlfit-find-logistic-parameters (qdata pdata tdata
)
205 (let* ((a (math-nlfit-find-qmax qdata pdata tdata
))
207 (mapcar (lambda (q) (calcFunc-ln (math-sub (math-div a q
) 1)))
209 (bandc (math-nlfit-least-squares tdata newqdata
)))
212 (calcFunc-exp (nth 0 bandc
))
215 ;;; Next, given the pdata and tdata, we can find the qdata if we know q0.
216 ;;; We first try to find q0, using the fact that when p takes on its largest
217 ;;; value, q is half of its maximum value. So we'll find the maximum value
218 ;;; of q given various q0, and use bisection to approximate the correct q0.
220 ;;; First, given pdata and tdata, find what half of qmax would be if q0=0.
222 (defun math-nlfit-find-qmaxhalf (pdata tdata
)
223 (let ((pmax (math-max-list (car pdata
) (cdr pdata
)))
225 (while (math-lessp (car pdata
) pmax
)
231 (math-add (nth 1 pdata
) (nth 0 pdata
)))
232 (math-sub (nth 1 tdata
)
234 (setq pdata
(cdr pdata
))
235 (setq tdata
(cdr tdata
)))
238 ;;; Next, given pdata and tdata, approximate q0.
240 (defun math-nlfit-find-q0 (pdata tdata
)
241 (let* ((qhalf (math-nlfit-find-qmaxhalf pdata tdata
))
242 (q0 (math-mul 2 qhalf
))
243 (qdata (math-nlfit-get-cumul-from-rates tdata pdata q0
)))
244 (while (math-lessp (math-nlfit-find-qmax
246 (lambda (q) (math-add q0 q
))
254 (setq q0
(math-add q0 qhalf
)))
255 (let* ((qmin (math-sub q0 qhalf
))
257 (qt (math-nlfit-find-qmax
259 (lambda (q) (math-add q0 q
))
264 (setq q0
(math-mul '(float 5 -
1) (math-add qmin qmax
)))
266 (math-nlfit-find-qmax
268 (lambda (q) (math-add q0 q
))
271 (math-mul '(float 5 -
1) (math-add qhalf q0
)))
275 (math-mul '(float 5 -
1) (math-add qmin qmax
)))))
277 ;;; To improve the approximations to the parameters, we can use
278 ;;; Marquardt method as described in Schwarz's book.
280 ;;; Small numbers used in the Givens algorithm
281 (defvar math-nlfit-delta
'(float 1 -
8))
283 (defvar math-nlfit-epsilon
'(float 1 -
5))
285 ;;; Maximum number of iterations
286 (defvar math-nlfit-max-its
100)
288 ;;; Next, we need some functions for dealing with vectors and
289 ;;; matrices. For convenience, we'll work with Emacs lists
290 ;;; as vectors, rather than Calc's vectors.
292 (defun math-nlfit-set-elt (vec i x
)
293 (setcar (nthcdr (1- i
) vec
) x
))
295 (defun math-nlfit-get-elt (vec i
)
298 (defun math-nlfit-make-matrix (i j
)
299 (let ((row (make-list j
0))
303 (setq mat
(cons (copy-sequence row
) mat
))
307 (defun math-nlfit-set-matx-elt (mat i j x
)
308 (setcar (nthcdr (1- j
) (nth (1- i
) mat
)) x
))
310 (defun math-nlfit-get-matx-elt (mat i j
)
311 (nth (1- j
) (nth (1- i
) mat
)))
313 ;;; For solving the linearized system.
314 ;;; (The Givens method, from Schwarz.)
316 (defun math-nlfit-givens (C d
)
317 (let* ((C (copy-tree C
))
331 (let ((cij (math-nlfit-get-matx-elt C i j
))
332 (cjj (math-nlfit-get-matx-elt C j j
)))
333 (when (not (math-equal 0 cij
))
334 (if (math-lessp (calcFunc-abs cjj
)
335 (math-mul math-nlfit-delta
(calcFunc-abs cij
)))
336 (setq w
(math-neg cij
)
345 (math-mul cij cij
))))
346 gamma
(math-div cjj w
)
347 sigma
(math-neg (math-div cij w
)))
348 (if (math-lessp (calcFunc-abs sigma
) gamma
)
350 (setq rho
(math-div (calcFunc-sign sigma
) gamma
))))
353 (math-nlfit-set-matx-elt C j j w
)
354 (math-nlfit-set-matx-elt C i j rho
)
357 (let* ((cjk (math-nlfit-get-matx-elt C j k
))
358 (cik (math-nlfit-get-matx-elt C i k
))
360 (math-mul gamma cjk
) (math-mul sigma cik
))))
363 (math-mul gamma cik
)))
365 (math-nlfit-set-matx-elt C i k cik
)
366 (math-nlfit-set-matx-elt C j k cjk
)
368 (let* ((di (math-nlfit-get-elt d i
))
369 (dj (math-nlfit-get-elt d j
))
372 (math-mul sigma di
))))
375 (math-mul gamma di
)))
377 (math-nlfit-set-elt d i di
)
378 (math-nlfit-set-elt d j dj
))))
384 (math-nlfit-set-elt r i
0)
385 (setq s
(math-nlfit-get-elt d i
))
388 (setq s
(math-add s
(math-mul (math-nlfit-get-matx-elt C i k
)
389 (math-nlfit-get-elt x k
))))
391 (math-nlfit-set-elt x i
394 (math-nlfit-get-matx-elt C i i
))))
398 (math-nlfit-set-elt r i
(math-nlfit-get-elt d i
))
404 (setq rho
(math-nlfit-get-matx-elt C i j
))
405 (if (math-equal rho
1)
408 (if (math-lessp (calcFunc-abs rho
) 1)
411 (math-sub 1 (math-mul sigma sigma
))))
412 (setq gamma
(math-div 1 (calcFunc-abs rho
))
413 sigma
(math-mul (calcFunc-sign rho
)
415 (math-sub 1 (math-mul gamma gamma
)))))))
416 (let ((ri (math-nlfit-get-elt r i
))
417 (rj (math-nlfit-get-elt r j
))
419 (setq h
(math-add (math-mul gamma rj
)
420 (math-mul sigma ri
)))
423 (math-mul sigma rj
)))
425 (math-nlfit-set-elt r i ri
)
426 (math-nlfit-set-elt r j rj
))
432 (defun math-nlfit-jacobian (grad xlist parms
&optional slist
)
435 (let ((row (apply grad
(car xlist
) parms
)))
439 (mapcar (lambda (x) (math-div x
(car slist
))) row
)
442 (setq slist
(cdr slist
))
443 (setq xlist
(cdr xlist
)))
446 (defun math-nlfit-make-ident (l n
)
447 (let ((m (math-nlfit-make-matrix n n
))
450 (math-nlfit-set-matx-elt m i i l
)
454 (defun math-nlfit-chi-sq (xlist ylist parms fn
&optional slist
)
459 (apply fn
(car xlist
) parms
)
462 (setq c
(math-div c
(car slist
))))
466 (setq xlist
(cdr xlist
))
467 (setq ylist
(cdr ylist
))
468 (setq slist
(cdr slist
)))
471 (defun math-nlfit-init-lambda (C)
478 (setq l
(math-add l
(math-mul (car row
) (car row
))))
479 (setq row
(cdr row
))))
481 (calcFunc-sqrt (math-div l
(math-mul n N
)))))
483 (defun math-nlfit-make-Ctilda (C l
)
484 (let* ((n (length (car C
)))
485 (bot (math-nlfit-make-ident l n
)))
488 (defun math-nlfit-make-d (fn xdata ydata parms
&optional sdata
)
492 (let ((dd (math-sub (apply fn
(car xdata
) parms
)
494 (if sdata
(math-div dd
(car sdata
)) dd
))
496 (setq xdata
(cdr xdata
))
497 (setq ydata
(cdr ydata
))
498 (setq sdata
(cdr sdata
)))
501 (defun math-nlfit-make-dtilda (d n
)
502 (append d
(make-list n
0)))
504 (defun math-nlfit-fit (xlist ylist parms fn grad
&optional slist
)
506 ((C (math-nlfit-jacobian grad xlist parms slist
))
507 (d (math-nlfit-make-d fn xlist ylist parms slist
))
508 (chisq (math-nlfit-chi-sq xlist ylist parms fn slist
))
509 (lambda (math-nlfit-init-lambda C
))
514 (< iters math-nlfit-max-its
))
515 (setq iters
(1+ iters
))
518 (let* ((Ctilda (math-nlfit-make-Ctilda C lambda
))
519 (dtilda (math-nlfit-make-dtilda d
(length (car C
))))
520 (zeta (math-nlfit-givens Ctilda dtilda
))
521 (newparms (math-map-binop 'math-add
(copy-tree parms
) zeta
))
522 (newchisq (math-nlfit-chi-sq xlist ylist newparms fn slist
)))
523 (if (math-lessp newchisq chisq
)
527 (math-sub chisq newchisq
) newchisq
) math-nlfit-epsilon
)
528 (setq really-done t
))
529 (setq lambda
(math-div lambda
10))
530 (setq chisq newchisq
)
531 (setq parms newparms
)
533 (setq lambda
(math-mul lambda
10)))))
534 (setq C
(math-nlfit-jacobian grad xlist parms slist
))
535 (setq d
(math-nlfit-make-d fn xlist ylist parms slist
))))
538 ;;; The functions that describe our models, and their gradients.
540 (defun math-nlfit-s-logistic-fn (x a b c
)
541 (math-div a
(math-add 1 (math-mul b
(calcFunc-exp (math-mul c x
))))))
543 (defun math-nlfit-s-logistic-grad (x a b c
)
544 (let* ((ep (calcFunc-exp (math-mul c x
)))
545 (d (math-add 1 (math-mul b ep
)))
549 (math-neg (math-div (math-mul a ep
) d2
))
550 (math-neg (math-div (math-mul a
(math-mul b
(math-mul x ep
))) d2
)))))
552 (defun math-nlfit-b-logistic-fn (x a c d
)
553 (let ((ex (calcFunc-exp (math-mul c
(math-sub x d
)))))
560 (defun math-nlfit-b-logistic-grad (x a c d
)
561 (let* ((ex (calcFunc-exp (math-mul c
(math-sub x d
))))
562 (ex1 (math-add 1 ex
))
570 (math-mul a
(math-mul xd ex
))
573 (math-mul 2 (math-mul a
(math-mul xd
(math-sqr ex
))))
577 (math-mul 2 (math-mul a
(math-mul c
(math-sqr ex
))))
580 (math-mul a
(math-mul c ex
))
583 ;;; Functions to get the final covariance matrix and the sdevs
585 (defun math-nlfit-find-covar (grad xlist pparms
)
588 (setq j
(cons (cons 'vec
(apply grad
(car xlist
) pparms
)) j
))
589 (setq xlist
(cdr xlist
)))
590 (setq j
(cons 'vec
(reverse j
)))
596 (defun math-nlfit-get-sigmas (grad xlist pparms chisq
)
598 (covar (math-nlfit-find-covar grad xlist pparms
))
599 (n (1- (length covar
)))
604 (setq sgs
(cons (calcFunc-sqrt (nth i
(nth i covar
))) sgs
))
606 (setq sgs
(reverse sgs
)))
609 ;;; Now the Calc functions
611 (defun math-nlfit-s-logistic-params (xdata ydata
)
612 (let ((pdata (math-nlfit-get-rates-from-cumul xdata ydata
)))
613 (math-nlfit-find-logistic-parameters ydata pdata xdata
)))
615 (defun math-nlfit-b-logistic-params (xdata ydata
)
616 (let* ((q0 (math-nlfit-find-q0 ydata xdata
))
617 (qdata (math-nlfit-get-cumul-from-rates xdata ydata q0
))
618 (abc (math-nlfit-find-logistic-parameters qdata ydata xdata
))
625 (D (math-neg (math-div (calcFunc-ln B
) C
)))
629 ;;; Some functions to turn the parameter lists and variables
630 ;;; into the appropriate functions.
632 (defun math-nlfit-s-logistic-solnexpr (pms var
)
633 (let ((a (nth 0 pms
))
646 (defun math-nlfit-b-logistic-solnexpr (pms var
)
647 (let ((a (nth 0 pms
))
666 (defun math-nlfit-enter-result (n prefix vals
)
667 (setq calc-aborted-prefix prefix
)
668 (calc-pop-push-record-list n prefix vals
)
671 (defun math-nlfit-fit-curve (fn grad solnexpr initparms
&optional sdv
)
673 (let* ((sdevv (or (eq sdv
'calcFunc-efit
) (eq sdv
'calcFunc-xfit
)))
674 (calc-display-working-message nil
)
676 (xdata (cdr (car (cdr data
))))
677 (ydata (cdr (car (cdr (cdr data
)))))
678 (sdata (if (math-contains-sdev-p ydata
)
679 (mapcar (lambda (x) (math-get-sdev x t
)) ydata
)
681 (ydata (mapcar (lambda (x) (math-get-value x
)) ydata
))
682 (calc-curve-varnames nil
)
683 (calc-curve-coefnames nil
)
685 (fitvars (calc-get-fit-variables 1 3))
686 (var (nth 1 calc-curve-varnames
))
687 (parms (cdr calc-curve-coefnames
))
689 (funcall initparms xdata ydata
))
690 (fit (math-nlfit-fit xdata ydata parmguess fn grad sdata
))
691 (finalparms (nth 1 fit
))
694 (math-nlfit-get-sigmas grad xdata finalparms
(nth 0 fit
))))
701 (lambda (x y
) (list 'sdev x y
)) finalparms sigmas
)
703 (soln (funcall solnexpr finalparms var
)))
704 (let ((calc-fit-to-trail t
)
707 (setq traillist
(cons (list 'calcFunc-eq
(car parms
) (car finalparms
))
709 (setq finalparms
(cdr finalparms
))
710 (setq parms
(cdr parms
)))
711 (setq traillist
(calc-normalize (cons 'vec
(nreverse traillist
))))
712 (cond ((eq sdv
'calcFunc-efit
)
713 (math-nlfit-enter-result 1 "efit" soln
))
714 ((eq sdv
'calcFunc-xfit
)
723 (let ((n (length xdata
))
724 (m (length finalparms
)))
725 (if (and sdata
(> n m
))
726 (calcFunc-utpc (nth 0 fit
)
728 '(var nan var-nan
)))))
729 (math-nlfit-enter-result 1 "xfit" sln
)))
731 (math-nlfit-enter-result 1 "fit" soln
)))
732 (calc-record traillist
"parm")))))
734 (defun calc-fit-s-shaped-logistic-curve (arg)
736 (math-nlfit-fit-curve 'math-nlfit-s-logistic-fn
737 'math-nlfit-s-logistic-grad
738 'math-nlfit-s-logistic-solnexpr
739 'math-nlfit-s-logistic-params
742 (defun calc-fit-bell-shaped-logistic-curve (arg)
744 (math-nlfit-fit-curve 'math-nlfit-b-logistic-fn
745 'math-nlfit-b-logistic-grad
746 'math-nlfit-b-logistic-solnexpr
747 'math-nlfit-b-logistic-params
750 (defun calc-fit-hubbert-linear-curve (&optional sdv
)
752 (let* ((sdevv (or (eq sdv
'calcFunc-efit
) (eq sdv
'calcFunc-xfit
)))
753 (calc-display-working-message nil
)
755 (qdata (cdr (car (cdr data
))))
756 (pdata (cdr (car (cdr (cdr data
)))))
757 (sdata (if (math-contains-sdev-p pdata
)
758 (mapcar (lambda (x) (math-get-sdev x t
)) pdata
)
760 (pdata (mapcar (lambda (x) (math-get-value x
)) pdata
))
761 (poverqdata (math-map-binop 'math-div pdata qdata
))
762 (parmvals (math-nlfit-least-squares qdata poverqdata sdata sdevv
))
763 (finalparms (list (nth 0 parmvals
)
765 (math-div (nth 0 parmvals
)
767 (calc-curve-varnames nil
)
768 (calc-curve-coefnames nil
)
770 (fitvars (calc-get-fit-variables 1 2))
771 (var (nth 1 calc-curve-varnames
))
772 (parms (cdr calc-curve-coefnames
))
773 (soln (list '* (nth 0 finalparms
)
775 (list '/ var
(nth 1 finalparms
))))))
776 (let ((calc-fit-to-trail t
)
780 (list 'calcFunc-eq
(nth 0 parms
) (nth 0 finalparms
))
781 (list 'calcFunc-eq
(nth 1 parms
) (nth 1 finalparms
))))
782 (cond ((eq sdv
'calcFunc-efit
)
783 (math-nlfit-enter-result 1 "efit" soln
))
784 ((eq sdv
'calcFunc-xfit
)
789 (list (nth 1 (nth 0 finalparms
))
790 (nth 1 (nth 1 finalparms
)))
804 '(calcFunc-fitdummy 1)
807 '(calcFunc-fitdummy 1)
808 '(calcFunc-fitdummy 2))))
810 (let ((n (length qdata
)))
811 (if (and sdata
(> n
2))
815 '(var nan var-nan
)))))
816 (math-nlfit-enter-result 1 "xfit" sln
)))
818 (math-nlfit-enter-result 1 "fit" soln
)))
819 (calc-record traillist
"parm")))))
821 (provide 'calc-nlfit
)
823 ;; arch-tag: 6eba3cd6-f48b-4a84-8174-10c15a024928