Mercurial > emacs
changeset 82259:4e7d189611f2
(calc-curve-fit): Add support for nonlinear curves and plotting.
| author | Jay Belanger <jay.p.belanger@gmail.com> |
|---|---|
| date | Sat, 04 Aug 2007 03:52:06 +0000 |
| parents | 9f8861776d74 |
| children | c647b5ba2a46 |
| files | lisp/calc/calcalg3.el |
| diffstat | 1 files changed, 140 insertions(+), 62 deletions(-) [+] |
line wrap: on
line diff
--- a/lisp/calc/calcalg3.el Sat Aug 04 02:31:52 2007 +0000 +++ b/lisp/calc/calcalg3.el Sat Aug 04 03:52:06 2007 +0000 @@ -115,6 +115,8 @@ (if (calc-is-hyperbolic) 'calcFunc-efit 'calcFunc-fit))) key (which 0) + (nonlinear nil) + (plot nil) n calc-curve-nvars temp data (homog nil) (msgs '( "(Press ? for help)" @@ -125,7 +127,11 @@ "E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)" "q = a + b (x-c)^2" "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)" + "s = a/(1 + exp(b (x - c)))" + "b = a exp(b (x - c))/(1 + exp(b (x - c)))^2" + "o = (y/x) = a (1 - x/b)" "h prefix = homogeneous model (no constant term)" + "P prefix = plot result" "' = alg entry, $ = stack, u = Model1, U = Model2"))) (while (not calc-curve-model) (message "Fit to model: %s:%s" @@ -138,6 +144,14 @@ (setq which (% (1+ which) (length msgs)))) ((memq key '(?h ?H)) (setq homog (not homog))) + ((= key ?P) + (let ((data (calc-top 1))) + (if (or + (calc-is-hyperbolic) + (calc-is-inverse) + (not (= (length data) 3))) + (setq plot "Can't plot") + (setq plot data)))) ((progn (if (eq key ?\$) (setq n 1) @@ -164,8 +178,9 @@ ((= key ?1) ; linear or multilinear (calc-get-fit-variables calc-curve-nvars (1+ calc-curve-nvars) (and homog 0)) - (setq calc-curve-model (math-mul calc-curve-coefnames - (cons 'vec (cons 1 (cdr calc-curve-varnames)))))) + (setq calc-curve-model + (math-mul calc-curve-coefnames + (cons 'vec (cons 1 (cdr calc-curve-varnames)))))) ((and (>= key ?2) (<= key ?9)) ; polynomial (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0)) (setq calc-curve-model @@ -180,58 +195,88 @@ ((= key ?p) ; power law (calc-get-fit-variables calc-curve-nvars (1+ calc-curve-nvars) (and homog 1)) - (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) - (calcFunc-reduce - '(var mul var-mul) - (calcFunc-map - '(var pow var-pow) - calc-curve-varnames - (cons 'vec (cdr (cdr calc-curve-coefnames)))))))) + (setq calc-curve-model + (math-mul + (nth 1 calc-curve-coefnames) + (calcFunc-reduce + '(var mul var-mul) + (calcFunc-map + '(var pow var-pow) + calc-curve-varnames + (cons 'vec (cdr (cdr calc-curve-coefnames)))))))) ((= key ?^) ; exponential law (calc-get-fit-variables calc-curve-nvars (1+ calc-curve-nvars) (and homog 1)) - (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) - (calcFunc-reduce - '(var mul var-mul) - (calcFunc-map - '(var pow var-pow) - (cons 'vec (cdr (cdr calc-curve-coefnames))) - calc-curve-varnames))))) + (setq calc-curve-model + (math-mul (nth 1 calc-curve-coefnames) + (calcFunc-reduce + '(var mul var-mul) + (calcFunc-map + '(var pow var-pow) + (cons 'vec (cdr (cdr calc-curve-coefnames))) + calc-curve-varnames))))) + ((= key ?s) + (setq nonlinear t) + (setq calc-curve-model t) + (require 'calc-nlfit) + (calc-fit-s-shaped-logistic-curve func)) + ((= key ?b) + (setq nonlinear t) + (setq calc-curve-model t) + (require 'calc-nlfit) + (calc-fit-bell-shaped-logistic-curve func)) + ((= key ?o) + (setq nonlinear t) + (setq calc-curve-model t) + (require 'calc-nlfit) + (if (and plot (not (stringp plot))) + (setq plot + (list 'vec + (nth 1 plot) + (cons + 'vec + (mapcar* 'calcFunc-div + (cdr (nth 2 plot)) + (cdr (nth 1 plot))))))) + (calc-fit-hubbert-linear-curve func)) ((memq key '(?e ?E)) (calc-get-fit-variables calc-curve-nvars (1+ calc-curve-nvars) (and homog 1)) - (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) - (calcFunc-reduce - '(var mul var-mul) - (calcFunc-map - (if (eq key ?e) - '(var exp var-exp) - '(calcFunc-lambda - (var a var-a) - (^ 10 (var a var-a)))) - (calcFunc-map - '(var mul var-mul) - (cons 'vec (cdr (cdr calc-curve-coefnames))) - calc-curve-varnames)))))) + (setq calc-curve-model + (math-mul (nth 1 calc-curve-coefnames) + (calcFunc-reduce + '(var mul var-mul) + (calcFunc-map + (if (eq key ?e) + '(var exp var-exp) + '(calcFunc-lambda + (var a var-a) + (^ 10 (var a var-a)))) + (calcFunc-map + '(var mul var-mul) + (cons 'vec (cdr (cdr calc-curve-coefnames))) + calc-curve-varnames)))))) ((memq key '(?x ?X)) (calc-get-fit-variables calc-curve-nvars (1+ calc-curve-nvars) (and homog 0)) - (setq calc-curve-model (math-mul calc-curve-coefnames - (cons 'vec (cons 1 (cdr calc-curve-varnames))))) + (setq calc-curve-model + (math-mul calc-curve-coefnames + (cons 'vec (cons 1 (cdr calc-curve-varnames))))) (setq calc-curve-model (if (eq key ?x) (list 'calcFunc-exp calc-curve-model) (list '^ 10 calc-curve-model)))) ((memq key '(?l ?L)) (calc-get-fit-variables calc-curve-nvars (1+ calc-curve-nvars) (and homog 0)) - (setq calc-curve-model (math-mul calc-curve-coefnames - (cons 'vec - (cons 1 (cdr (calcFunc-map - (if (eq key ?l) - '(var ln var-ln) - '(var log10 - var-log10)) - calc-curve-varnames))))))) + (setq calc-curve-model + (math-mul calc-curve-coefnames + (cons 'vec + (cons 1 (cdr (calcFunc-map + (if (eq key ?l) + '(var ln var-ln) + '(var log10 + var-log10)) + calc-curve-varnames))))))) ((= key ?q) (calc-get-fit-variables calc-curve-nvars (1+ (* 2 calc-curve-nvars)) (and homog 0)) @@ -247,12 +292,14 @@ (list '- (car v) (nth 1 c)) 2))))))) ((= key ?g) - (setq calc-curve-model - (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)") - calc-curve-varnames '(vec (var XFit var-XFit)) - calc-curve-coefnames '(vec (var AFit var-AFit) - (var BFit var-BFit) - (var CFit var-CFit))) + (setq + calc-curve-model + (math-read-expr + "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)") + calc-curve-varnames '(vec (var XFit var-XFit)) + calc-curve-coefnames '(vec (var AFit var-AFit) + (var BFit var-BFit) + (var CFit var-CFit))) (calc-get-fit-variables 1 (1- (length calc-curve-coefnames)) (and homog 1))) ((memq key '(?\$ ?\' ?u ?U)) @@ -262,8 +309,9 @@ (let* ((calc-dollar-values calc-arg-values) (calc-dollar-used 0) (calc-hashes-used 0)) - (setq calc-curve-model (calc-do-alg-entry "" "Model formula: " - nil 'calc-curve-fit-history)) + (setq calc-curve-model + (calc-do-alg-entry "" "Model formula: " + nil 'calc-curve-fit-history)) (if (/= (length calc-curve-model) 1) (error "Bad format")) (setq calc-curve-model (car calc-curve-model) @@ -296,11 +344,13 @@ (or (nth 3 calc-curve-model) (cons 'vec (math-all-vars-but - calc-curve-model calc-curve-varnames))) + calc-curve-model + calc-curve-varnames))) calc-curve-model (nth 1 calc-curve-model)) (error "Incorrect model specifier"))))) (or calc-curve-varnames - (let ((with-y (eq (car-safe calc-curve-model) 'calcFunc-eq))) + (let ((with-y + (eq (car-safe calc-curve-model) 'calcFunc-eq))) (if calc-curve-coefnames (calc-get-fit-variables (if with-y (1+ calc-curve-nvars) calc-curve-nvars) @@ -310,7 +360,10 @@ nil with-y) (let* ((coefs (math-all-vars-but calc-curve-model nil)) (vars nil) - (n (- (length coefs) calc-curve-nvars (if with-y 2 1))) + (n (- + (length coefs) + calc-curve-nvars + (if with-y 2 1))) p) (if (< n 0) (error "Not enough variables in model")) @@ -326,18 +379,43 @@ calc-curve-varnames calc-curve-coefnames) "modl")))) (t (beep)))) - (let ((calc-fit-to-trail t)) - (calc-enter-result n (substring (symbol-name func) 9) - (list func calc-curve-model - (if (= (length calc-curve-varnames) 2) - (nth 1 calc-curve-varnames) - calc-curve-varnames) - (if (= (length calc-curve-coefnames) 2) - (nth 1 calc-curve-coefnames) - calc-curve-coefnames) - data)) - (if (consp calc-fit-to-trail) - (calc-record (calc-normalize calc-fit-to-trail) "parm")))))) + (unless nonlinear + (let ((calc-fit-to-trail t)) + (calc-enter-result n (substring (symbol-name func) 9) + (list func calc-curve-model + (if (= (length calc-curve-varnames) 2) + (nth 1 calc-curve-varnames) + calc-curve-varnames) + (if (= (length calc-curve-coefnames) 2) + (nth 1 calc-curve-coefnames) + calc-curve-coefnames) + data)) + (if (consp calc-fit-to-trail) + (calc-record (calc-normalize calc-fit-to-trail) "parm")))) + (when plot + (if (stringp plot) + (message plot) + (let ((calc-graph-no-auto-view t)) + (calc-graph-delete t) + (calc-graph-add-curve + (calc-graph-lookup (nth 1 plot)) + (calc-graph-lookup (nth 2 plot))) + (unless (math-contains-sdev-p (nth 2 data)) + (calc-graph-set-styles nil nil) + (calc-graph-point-style nil)) + (setq plot (cdr (nth 1 plot))) + (setq plot + (list 'intv + 3 + (math-sub + (math-min-list (car plot) (cdr plot)) + '(float 5 -1)) + (math-add + '(float 5 -1) + (math-max-list (car plot) (cdr plot))))) + (calc-graph-add-curve (calc-graph-lookup plot) + (calc-graph-lookup (calc-top-n 1))) + (calc-graph-plot nil))))))) (defun calc-invent-independent-variables (n &optional but) (calc-invent-variables n but '(x y z t) "x"))