Function: math-general-fit
math-general-fit is a byte-compiled function defined in
calcalg3.el.gz.
Signature
(math-general-fit EXPR VARS COEFS DATA MODE)
Source Code
;; Defined in /usr/src/emacs/lisp/calc/calcalg3.el.gz
(defun math-general-fit (expr vars coefs data mode)
(defvar var-YVAL) (defvar var-YVALX)
(let ((calc-simplify-mode nil)
(math-dummy-counter math-dummy-counter)
(math-in-fit 1)
(extended (eq mode 'full))
(math-fit-first-coef math-dummy-counter)
math-fit-first-var
(plain-expr expr)
orig-expr
have-sdevs need-chisq chisq
(x-funcs nil)
(y-filter nil)
y-dummy
(coef-filters nil)
math-fit-new-coefs
(xy-values nil)
(weights nil)
(var-YVAL nil) (var-YVALX nil)
covar beta
n m mm v dummy p) ;; nn
;; Validate and parse arguments.
(or data
(if coefs
(setq data coefs
coefs nil)
(if (math-vectorp expr)
(if (memq (length expr) '(3 4))
(setq data vars
vars (nth 2 expr)
coefs (nth 3 expr)
expr (nth 1 expr))
(math-dimension-error))
(setq data vars
vars nil
coefs nil))))
(or (math-matrixp data) (math-reject-arg data 'matrixp))
(setq v (1- (length data))
n (1- (length (nth 1 data))))
(or (math-vectorp vars) (null vars)
(setq vars (list 'vec vars)))
(or (math-vectorp coefs) (null coefs)
(setq coefs (list 'vec coefs)))
(or coefs
(setq coefs (cons 'vec (math-all-vars-but expr vars))))
(or vars
(if (<= (1- (length coefs)) v)
(math-reject-arg coefs "*Not enough variables in model")
(setq coefs (copy-sequence coefs))
(let ((p (nthcdr (- (length coefs) v
(if (eq (car-safe expr) 'calcFunc-eq) 1 0))
coefs)))
(setq vars (cons 'vec (cdr p)))
(setcdr p nil))))
(or (= (1- (length vars)) v)
(= (length vars) v)
(math-reject-arg vars "*Number of variables does not match data"))
(setq m (1- (length coefs)))
(if (< m 1)
(math-reject-arg coefs "*Need at least one parameter"))
;; Rewrite expr in terms of fitparam and fitvar, make into an equation.
(setq p coefs)
(while (setq p (cdr p))
(or (eq (car-safe (car p)) 'var)
(math-reject-arg (car p) "*Expected a variable"))
(setq dummy (math-dummy-variable)
expr (math-expr-subst expr (car p)
(list 'calcFunc-fitparam
(- math-dummy-counter math-fit-first-coef)))))
(setq math-fit-first-var math-dummy-counter
p vars)
(while (setq p (cdr p))
(or (eq (car-safe (car p)) 'var)
(math-reject-arg (car p) "*Expected a variable"))
(setq dummy (math-dummy-variable)
expr (math-expr-subst expr (car p)
(list 'calcFunc-fitvar
(- math-dummy-counter math-fit-first-var)))))
(if (< math-dummy-counter (+ math-fit-first-var v))
(setq dummy (math-dummy-variable))) ; dependent variable may be unnamed
(setq y-dummy dummy
orig-expr expr)
(or (eq (car-safe expr) 'calcFunc-eq)
(setq expr (list 'calcFunc-eq (list 'calcFunc-fitvar v) expr)))
(let ((calc-symbolic-mode nil))
;; Apply rewrites to put expr into a linear-like form.
(setq expr (math-evaluate-expr expr)
expr (math-rewrite (list 'calcFunc-fitmodel expr)
'(var FitRules var-FitRules))
math-in-fit 2
expr (math-evaluate-expr expr))
(or (and (eq (car-safe expr) 'calcFunc-fitsystem)
(= (length expr) 4)
(math-vectorp (nth 2 expr))
(math-vectorp (nth 3 expr))
(> (length (nth 2 expr)) 1)
(= (length (nth 3 expr)) (1+ m)))
(math-reject-arg plain-expr "*Model expression is too complex"))
(setq y-filter (nth 1 expr)
x-funcs (vconcat (cdr (nth 2 expr)))
coef-filters (nth 3 expr)
mm (length x-funcs))
(if (equal y-filter y-dummy)
(setq y-filter nil))
;; Build the (square) system of linear equations to be solved.
(setq beta (cons 'vec (make-list mm 0))
covar (cons 'vec (mapcar 'copy-sequence (make-list mm beta))))
(let* ((ptrs (vconcat (cdr data)))
(isigsq 1)
(xvals (make-vector mm 0))
(i 0)
j k xval yval sigmasqr wt covj covjk covk betaj) ;; lud
(while (<= (setq i (1+ i)) n)
;; Assign various independent variables for this data point.
(setq j 0
sigmasqr nil)
(while (< j v)
(aset ptrs j (cdr (aref ptrs j)))
(setq xval (car (aref ptrs j)))
(if (= j (1- v))
(if sigmasqr
(progn
(if (eq (car-safe xval) 'sdev)
(setq sigmasqr (math-add (math-sqr (nth 2 xval))
sigmasqr)
xval (nth 1 xval)))
(if y-filter
(setq xval (math-make-sdev xval
(math-sqrt sigmasqr))))))
(if (eq (car-safe xval) 'sdev)
(setq sigmasqr (math-add (math-sqr (nth 2 xval))
(or sigmasqr 0))
xval (nth 1 xval))))
(set (nth 2 (aref math-dummy-vars (+ math-fit-first-var j))) xval)
(setq j (1+ j)))
;; Compute Y value for this data point.
(if y-filter
(setq yval (math-evaluate-expr y-filter))
(setq yval (symbol-value (nth 2 y-dummy))))
(if (eq (car-safe yval) 'sdev)
(setq sigmasqr (math-sqr (nth 2 yval))
yval (nth 1 yval)))
(if (= i 1)
(setq have-sdevs sigmasqr
need-chisq (or extended
(and (eq mode 'sdev) (not have-sdevs)))))
(if have-sdevs
(if sigmasqr
(progn
(setq isigsq (math-div 1 sigmasqr))
(if need-chisq
(setq weights (cons isigsq weights))))
(math-reject-arg yval "*Mixed error forms and plain numbers"))
(if sigmasqr
(math-reject-arg yval "*Mixed error forms and plain numbers")))
;; Compute X values for this data point and update covar and beta.
(if (eq (car-safe xval) 'sdev)
(set (nth 2 y-dummy) (nth 1 xval)))
(setq j 0
covj covar
betaj beta)
(while (< j mm)
(setq wt (math-evaluate-expr (aref x-funcs j)))
(aset xvals j wt)
(setq wt (math-mul wt isigsq)
betaj (cdr betaj)
covjk (car (setq covj (cdr covj)))
k 0)
(while (<= k j)
(setq covjk (cdr covjk))
(setcar covjk (math-add (car covjk)
(math-mul wt (aref xvals k))))
(setq k (1+ k)))
(setcar betaj (math-add (car betaj) (math-mul wt yval)))
(setq j (1+ j)))
(if need-chisq
(setq xy-values (cons (append xvals (list yval)) xy-values))))
;; Fill in symmetric half of covar matrix.
(setq j 0
covj covar)
(while (< j (1- mm))
(setq k j
j (1+ j)
covjk (nthcdr j (car (setq covj (cdr covj))))
covk (nthcdr j covar))
(while (< (setq k (1+ k)) mm)
(setq covjk (cdr covjk)
covk (cdr covk))
(setcar covjk (nth j (car covk))))))
;; Solve the linear system.
(if mode
(progn
(setq covar (math-matrix-inv-raw covar))
(if covar
(setq beta (math-mul covar beta))
(if (math-zerop (math-abs beta))
(setq covar (calcFunc-diag 0 (1- (length beta))))
(math-reject-arg orig-expr "*Singular matrix")))
(or (math-vectorp covar)
(setq covar (list 'vec (list 'vec covar)))))
(setq beta (math-div beta covar)))
;; Compute chi-square statistic if necessary.
(if need-chisq
(let (bp xp sum)
(setq chisq 0)
(while xy-values
(setq bp beta
xp (car xy-values)
sum 0)
(while (setq bp (cdr bp))
(setq sum (math-add sum (math-mul (car bp) (car xp)))
xp (cdr xp)))
(setq sum (math-sqr (math-sub (car xp) sum)))
(if weights (setq sum (math-mul sum (car weights))))
(setq chisq (math-add chisq sum)
weights (cdr weights)
xy-values (cdr xy-values)))))
;; Convert coefficients back into original terms.
(setq math-fit-new-coefs (copy-sequence beta))
(let* ((bp math-fit-new-coefs)
(cp covar)
(sigdat 1)
(math-in-fit 3)
(j 0))
(and mode (not have-sdevs)
(setq sigdat (if (<= n mm)
0
(math-div chisq (- n mm)))))
(if mode
(while (setq bp (cdr bp))
(setcar bp (math-make-sdev
(car bp)
(math-sqrt (math-mul (nth (setq j (1+ j))
(car (setq cp (cdr cp))))
sigdat))))))
(setq math-fit-new-coefs (math-evaluate-expr coef-filters))
(if calc-fit-to-trail
(let ((bp math-fit-new-coefs)
(cp coefs)
(vec nil))
(while (setq bp (cdr bp) cp (cdr cp))
(setq vec (cons (list 'calcFunc-eq (car cp) (car bp)) vec)))
(setq calc-fit-to-trail (cons 'vec (nreverse vec)))))))
;; Substitute best-fit coefficients back into original formula.
(setq expr (math-multi-subst
orig-expr
(let ((n v)
(vec nil))
(while (>= n 1)
(setq vec (cons (list 'calcFunc-fitvar n) vec)
n (1- n)))
(setq n m)
(while (>= n 1)
(setq vec (cons (list 'calcFunc-fitparam n) vec)
n (1- n)))
vec)
(append (cdr math-fit-new-coefs) (cdr vars))))
;; Package the result.
(math-normalize
(if extended
(list 'vec expr beta covar
(let ((p coef-filters)
(n 0))
(while (and (setq n (1+ n) p (cdr p))
(eq (car-safe (car p)) 'calcFunc-fitdummy)
(eq (nth 1 (car p)) n)))
(if p
coef-filters
(list 'vec)))
chisq
(if (and have-sdevs (> n mm))
(list 'calcFunc-utpc chisq (- n mm))
'(var nan var-nan)))
expr))))