aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--lisp/calc/calcalg3.el202
1 files changed, 140 insertions, 62 deletions
diff --git a/lisp/calc/calcalg3.el b/lisp/calc/calcalg3.el
index 9f263a2281a..611e2d7f45d 100644
--- a/lisp/calc/calcalg3.el
+++ b/lisp/calc/calcalg3.el
@@ -115,6 +115,8 @@
115 (if (calc-is-hyperbolic) 'calcFunc-efit 115 (if (calc-is-hyperbolic) 'calcFunc-efit
116 'calcFunc-fit))) 116 'calcFunc-fit)))
117 key (which 0) 117 key (which 0)
118 (nonlinear nil)
119 (plot nil)
118 n calc-curve-nvars temp data 120 n calc-curve-nvars temp data
119 (homog nil) 121 (homog nil)
120 (msgs '( "(Press ? for help)" 122 (msgs '( "(Press ? for help)"
@@ -125,7 +127,11 @@
125 "E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)" 127 "E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)"
126 "q = a + b (x-c)^2" 128 "q = a + b (x-c)^2"
127 "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)" 129 "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
130 "s = a/(1 + exp(b (x - c)))"
131 "b = a exp(b (x - c))/(1 + exp(b (x - c)))^2"
132 "o = (y/x) = a (1 - x/b)"
128 "h prefix = homogeneous model (no constant term)" 133 "h prefix = homogeneous model (no constant term)"
134 "P prefix = plot result"
129 "' = alg entry, $ = stack, u = Model1, U = Model2"))) 135 "' = alg entry, $ = stack, u = Model1, U = Model2")))
130 (while (not calc-curve-model) 136 (while (not calc-curve-model)
131 (message "Fit to model: %s:%s" 137 (message "Fit to model: %s:%s"
@@ -138,6 +144,14 @@
138 (setq which (% (1+ which) (length msgs)))) 144 (setq which (% (1+ which) (length msgs))))
139 ((memq key '(?h ?H)) 145 ((memq key '(?h ?H))
140 (setq homog (not homog))) 146 (setq homog (not homog)))
147 ((= key ?P)
148 (let ((data (calc-top 1)))
149 (if (or
150 (calc-is-hyperbolic)
151 (calc-is-inverse)
152 (not (= (length data) 3)))
153 (setq plot "Can't plot")
154 (setq plot data))))
141 ((progn 155 ((progn
142 (if (eq key ?\$) 156 (if (eq key ?\$)
143 (setq n 1) 157 (setq n 1)
@@ -164,8 +178,9 @@
164 ((= key ?1) ; linear or multilinear 178 ((= key ?1) ; linear or multilinear
165 (calc-get-fit-variables calc-curve-nvars 179 (calc-get-fit-variables calc-curve-nvars
166 (1+ calc-curve-nvars) (and homog 0)) 180 (1+ calc-curve-nvars) (and homog 0))
167 (setq calc-curve-model (math-mul calc-curve-coefnames 181 (setq calc-curve-model
168 (cons 'vec (cons 1 (cdr calc-curve-varnames)))))) 182 (math-mul calc-curve-coefnames
183 (cons 'vec (cons 1 (cdr calc-curve-varnames))))))
169 ((and (>= key ?2) (<= key ?9)) ; polynomial 184 ((and (>= key ?2) (<= key ?9)) ; polynomial
170 (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0)) 185 (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
171 (setq calc-curve-model 186 (setq calc-curve-model
@@ -180,58 +195,88 @@
180 ((= key ?p) ; power law 195 ((= key ?p) ; power law
181 (calc-get-fit-variables calc-curve-nvars 196 (calc-get-fit-variables calc-curve-nvars
182 (1+ calc-curve-nvars) (and homog 1)) 197 (1+ calc-curve-nvars) (and homog 1))
183 (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) 198 (setq calc-curve-model
184 (calcFunc-reduce 199 (math-mul
185 '(var mul var-mul) 200 (nth 1 calc-curve-coefnames)
186 (calcFunc-map 201 (calcFunc-reduce
187 '(var pow var-pow) 202 '(var mul var-mul)
188 calc-curve-varnames 203 (calcFunc-map
189 (cons 'vec (cdr (cdr calc-curve-coefnames)))))))) 204 '(var pow var-pow)
205 calc-curve-varnames
206 (cons 'vec (cdr (cdr calc-curve-coefnames))))))))
190 ((= key ?^) ; exponential law 207 ((= key ?^) ; exponential law
191 (calc-get-fit-variables calc-curve-nvars 208 (calc-get-fit-variables calc-curve-nvars
192 (1+ calc-curve-nvars) (and homog 1)) 209 (1+ calc-curve-nvars) (and homog 1))
193 (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) 210 (setq calc-curve-model
194 (calcFunc-reduce 211 (math-mul (nth 1 calc-curve-coefnames)
195 '(var mul var-mul) 212 (calcFunc-reduce
196 (calcFunc-map 213 '(var mul var-mul)
197 '(var pow var-pow) 214 (calcFunc-map
198 (cons 'vec (cdr (cdr calc-curve-coefnames))) 215 '(var pow var-pow)
199 calc-curve-varnames))))) 216 (cons 'vec (cdr (cdr calc-curve-coefnames)))
217 calc-curve-varnames)))))
218 ((= key ?s)
219 (setq nonlinear t)
220 (setq calc-curve-model t)
221 (require 'calc-nlfit)
222 (calc-fit-s-shaped-logistic-curve func))
223 ((= key ?b)
224 (setq nonlinear t)
225 (setq calc-curve-model t)
226 (require 'calc-nlfit)
227 (calc-fit-bell-shaped-logistic-curve func))
228 ((= key ?o)
229 (setq nonlinear t)
230 (setq calc-curve-model t)
231 (require 'calc-nlfit)
232 (if (and plot (not (stringp plot)))
233 (setq plot
234 (list 'vec
235 (nth 1 plot)
236 (cons
237 'vec
238 (mapcar* 'calcFunc-div
239 (cdr (nth 2 plot))
240 (cdr (nth 1 plot)))))))
241 (calc-fit-hubbert-linear-curve func))
200 ((memq key '(?e ?E)) 242 ((memq key '(?e ?E))
201 (calc-get-fit-variables calc-curve-nvars 243 (calc-get-fit-variables calc-curve-nvars
202 (1+ calc-curve-nvars) (and homog 1)) 244 (1+ calc-curve-nvars) (and homog 1))
203 (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) 245 (setq calc-curve-model
204 (calcFunc-reduce 246 (math-mul (nth 1 calc-curve-coefnames)
205 '(var mul var-mul) 247 (calcFunc-reduce
206 (calcFunc-map 248 '(var mul var-mul)
207 (if (eq key ?e) 249 (calcFunc-map
208 '(var exp var-exp) 250 (if (eq key ?e)
209 '(calcFunc-lambda 251 '(var exp var-exp)
210 (var a var-a) 252 '(calcFunc-lambda
211 (^ 10 (var a var-a)))) 253 (var a var-a)
212 (calcFunc-map 254 (^ 10 (var a var-a))))
213 '(var mul var-mul) 255 (calcFunc-map
214 (cons 'vec (cdr (cdr calc-curve-coefnames))) 256 '(var mul var-mul)
215 calc-curve-varnames)))))) 257 (cons 'vec (cdr (cdr calc-curve-coefnames)))
258 calc-curve-varnames))))))
216 ((memq key '(?x ?X)) 259 ((memq key '(?x ?X))
217 (calc-get-fit-variables calc-curve-nvars 260 (calc-get-fit-variables calc-curve-nvars
218 (1+ calc-curve-nvars) (and homog 0)) 261 (1+ calc-curve-nvars) (and homog 0))
219 (setq calc-curve-model (math-mul calc-curve-coefnames 262 (setq calc-curve-model
220 (cons 'vec (cons 1 (cdr calc-curve-varnames))))) 263 (math-mul calc-curve-coefnames
264 (cons 'vec (cons 1 (cdr calc-curve-varnames)))))
221 (setq calc-curve-model (if (eq key ?x) 265 (setq calc-curve-model (if (eq key ?x)
222 (list 'calcFunc-exp calc-curve-model) 266 (list 'calcFunc-exp calc-curve-model)
223 (list '^ 10 calc-curve-model)))) 267 (list '^ 10 calc-curve-model))))
224 ((memq key '(?l ?L)) 268 ((memq key '(?l ?L))
225 (calc-get-fit-variables calc-curve-nvars 269 (calc-get-fit-variables calc-curve-nvars
226 (1+ calc-curve-nvars) (and homog 0)) 270 (1+ calc-curve-nvars) (and homog 0))
227 (setq calc-curve-model (math-mul calc-curve-coefnames 271 (setq calc-curve-model
228 (cons 'vec 272 (math-mul calc-curve-coefnames
229 (cons 1 (cdr (calcFunc-map 273 (cons 'vec
230 (if (eq key ?l) 274 (cons 1 (cdr (calcFunc-map
231 '(var ln var-ln) 275 (if (eq key ?l)
232 '(var log10 276 '(var ln var-ln)
233 var-log10)) 277 '(var log10
234 calc-curve-varnames))))))) 278 var-log10))
279 calc-curve-varnames)))))))
235 ((= key ?q) 280 ((= key ?q)
236 (calc-get-fit-variables calc-curve-nvars 281 (calc-get-fit-variables calc-curve-nvars
237 (1+ (* 2 calc-curve-nvars)) (and homog 0)) 282 (1+ (* 2 calc-curve-nvars)) (and homog 0))
@@ -247,12 +292,14 @@
247 (list '- (car v) (nth 1 c)) 292 (list '- (car v) (nth 1 c))
248 2))))))) 293 2)))))))
249 ((= key ?g) 294 ((= key ?g)
250 (setq calc-curve-model 295 (setq
251 (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)") 296 calc-curve-model
252 calc-curve-varnames '(vec (var XFit var-XFit)) 297 (math-read-expr
253 calc-curve-coefnames '(vec (var AFit var-AFit) 298 "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
254 (var BFit var-BFit) 299 calc-curve-varnames '(vec (var XFit var-XFit))
255 (var CFit var-CFit))) 300 calc-curve-coefnames '(vec (var AFit var-AFit)
301 (var BFit var-BFit)
302 (var CFit var-CFit)))
256 (calc-get-fit-variables 1 (1- (length calc-curve-coefnames)) 303 (calc-get-fit-variables 1 (1- (length calc-curve-coefnames))
257 (and homog 1))) 304 (and homog 1)))
258 ((memq key '(?\$ ?\' ?u ?U)) 305 ((memq key '(?\$ ?\' ?u ?U))
@@ -262,8 +309,9 @@
262 (let* ((calc-dollar-values calc-arg-values) 309 (let* ((calc-dollar-values calc-arg-values)
263 (calc-dollar-used 0) 310 (calc-dollar-used 0)
264 (calc-hashes-used 0)) 311 (calc-hashes-used 0))
265 (setq calc-curve-model (calc-do-alg-entry "" "Model formula: " 312 (setq calc-curve-model
266 nil 'calc-curve-fit-history)) 313 (calc-do-alg-entry "" "Model formula: "
314 nil 'calc-curve-fit-history))
267 (if (/= (length calc-curve-model) 1) 315 (if (/= (length calc-curve-model) 1)
268 (error "Bad format")) 316 (error "Bad format"))
269 (setq calc-curve-model (car calc-curve-model) 317 (setq calc-curve-model (car calc-curve-model)
@@ -296,11 +344,13 @@
296 (or (nth 3 calc-curve-model) 344 (or (nth 3 calc-curve-model)
297 (cons 'vec 345 (cons 'vec
298 (math-all-vars-but 346 (math-all-vars-but
299 calc-curve-model calc-curve-varnames))) 347 calc-curve-model
348 calc-curve-varnames)))
300 calc-curve-model (nth 1 calc-curve-model)) 349 calc-curve-model (nth 1 calc-curve-model))
301 (error "Incorrect model specifier"))))) 350 (error "Incorrect model specifier")))))
302 (or calc-curve-varnames 351 (or calc-curve-varnames
303 (let ((with-y (eq (car-safe calc-curve-model) 'calcFunc-eq))) 352 (let ((with-y
353 (eq (car-safe calc-curve-model) 'calcFunc-eq)))
304 (if calc-curve-coefnames 354 (if calc-curve-coefnames
305 (calc-get-fit-variables 355 (calc-get-fit-variables
306 (if with-y (1+ calc-curve-nvars) calc-curve-nvars) 356 (if with-y (1+ calc-curve-nvars) calc-curve-nvars)
@@ -310,7 +360,10 @@
310 nil with-y) 360 nil with-y)
311 (let* ((coefs (math-all-vars-but calc-curve-model nil)) 361 (let* ((coefs (math-all-vars-but calc-curve-model nil))
312 (vars nil) 362 (vars nil)
313 (n (- (length coefs) calc-curve-nvars (if with-y 2 1))) 363 (n (-
364 (length coefs)
365 calc-curve-nvars
366 (if with-y 2 1)))
314 p) 367 p)
315 (if (< n 0) 368 (if (< n 0)
316 (error "Not enough variables in model")) 369 (error "Not enough variables in model"))
@@ -326,18 +379,43 @@
326 calc-curve-varnames calc-curve-coefnames) 379 calc-curve-varnames calc-curve-coefnames)
327 "modl")))) 380 "modl"))))
328 (t (beep)))) 381 (t (beep))))
329 (let ((calc-fit-to-trail t)) 382 (unless nonlinear
330 (calc-enter-result n (substring (symbol-name func) 9) 383 (let ((calc-fit-to-trail t))
331 (list func calc-curve-model 384 (calc-enter-result n (substring (symbol-name func) 9)
332 (if (= (length calc-curve-varnames) 2) 385 (list func calc-curve-model
333 (nth 1 calc-curve-varnames) 386 (if (= (length calc-curve-varnames) 2)
334 calc-curve-varnames) 387 (nth 1 calc-curve-varnames)
335 (if (= (length calc-curve-coefnames) 2) 388 calc-curve-varnames)
336 (nth 1 calc-curve-coefnames) 389 (if (= (length calc-curve-coefnames) 2)
337 calc-curve-coefnames) 390 (nth 1 calc-curve-coefnames)
338 data)) 391 calc-curve-coefnames)
339 (if (consp calc-fit-to-trail) 392 data))
340 (calc-record (calc-normalize calc-fit-to-trail) "parm")))))) 393 (if (consp calc-fit-to-trail)
394 (calc-record (calc-normalize calc-fit-to-trail) "parm"))))
395 (when plot
396 (if (stringp plot)
397 (message plot)
398 (let ((calc-graph-no-auto-view t))
399 (calc-graph-delete t)
400 (calc-graph-add-curve
401 (calc-graph-lookup (nth 1 plot))
402 (calc-graph-lookup (nth 2 plot)))
403 (unless (math-contains-sdev-p (nth 2 data))
404 (calc-graph-set-styles nil nil)
405 (calc-graph-point-style nil))
406 (setq plot (cdr (nth 1 plot)))
407 (setq plot
408 (list 'intv
409 3
410 (math-sub
411 (math-min-list (car plot) (cdr plot))
412 '(float 5 -1))
413 (math-add
414 '(float 5 -1)
415 (math-max-list (car plot) (cdr plot)))))
416 (calc-graph-add-curve (calc-graph-lookup plot)
417 (calc-graph-lookup (calc-top-n 1)))
418 (calc-graph-plot nil)))))))
341 419
342(defun calc-invent-independent-variables (n &optional but) 420(defun calc-invent-independent-variables (n &optional but)
343 (calc-invent-variables n but '(x y z t) "x")) 421 (calc-invent-variables n but '(x y z t) "x"))