40785
|
1 ;; Calculator for GNU Emacs, part II [calc-stat.el]
|
|
2 ;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
|
|
3 ;; Written by Dave Gillespie, daveg@synaptics.com.
|
|
4
|
|
5 ;; This file is part of GNU Emacs.
|
|
6
|
|
7 ;; GNU Emacs is distributed in the hope that it will be useful,
|
|
8 ;; but WITHOUT ANY WARRANTY. No author or distributor
|
|
9 ;; accepts responsibility to anyone for the consequences of using it
|
|
10 ;; or for whether it serves any particular purpose or works at all,
|
|
11 ;; unless he says so in writing. Refer to the GNU Emacs General Public
|
|
12 ;; License for full details.
|
|
13
|
|
14 ;; Everyone is granted permission to copy, modify and redistribute
|
|
15 ;; GNU Emacs, but only under the conditions described in the
|
|
16 ;; GNU Emacs General Public License. A copy of this license is
|
|
17 ;; supposed to have been given to you along with GNU Emacs so you
|
|
18 ;; can know your rights and responsibilities. It should be in a
|
|
19 ;; file named COPYING. Among other things, the copyright notice
|
|
20 ;; and this notice must be preserved on all copies.
|
|
21
|
|
22
|
|
23
|
|
24 ;; This file is autoloaded from calc-ext.el.
|
|
25 (require 'calc-ext)
|
|
26
|
|
27 (require 'calc-macs)
|
|
28
|
|
29 (defun calc-Need-calc-stat () nil)
|
|
30
|
|
31
|
|
32 ;;; Statistical operations on vectors.
|
|
33
|
|
34 (defun calc-vector-count (arg)
|
|
35 (interactive "P")
|
|
36 (calc-slow-wrapper
|
|
37 (calc-vector-op "coun" 'calcFunc-vcount arg))
|
|
38 )
|
|
39
|
|
40 (defun calc-vector-sum (arg)
|
|
41 (interactive "P")
|
|
42 (calc-slow-wrapper
|
|
43 (if (calc-is-hyperbolic)
|
|
44 (calc-vector-op "vprd" 'calcFunc-vprod arg)
|
|
45 (calc-vector-op "vsum" 'calcFunc-vsum arg)))
|
|
46 )
|
|
47
|
|
48 (defun calc-vector-product (arg)
|
|
49 (interactive "P")
|
|
50 (calc-hyperbolic-func)
|
|
51 (calc-vector-sum arg)
|
|
52 )
|
|
53
|
|
54 (defun calc-vector-max (arg)
|
|
55 (interactive "P")
|
|
56 (calc-slow-wrapper
|
|
57 (if (calc-is-inverse)
|
|
58 (calc-vector-op "vmin" 'calcFunc-vmin arg)
|
|
59 (calc-vector-op "vmax" 'calcFunc-vmax arg)))
|
|
60 )
|
|
61
|
|
62 (defun calc-vector-min (arg)
|
|
63 (interactive "P")
|
|
64 (calc-invert-func)
|
|
65 (calc-vector-max arg)
|
|
66 )
|
|
67
|
|
68 (defun calc-vector-mean (arg)
|
|
69 (interactive "P")
|
|
70 (calc-slow-wrapper
|
|
71 (if (calc-is-hyperbolic)
|
|
72 (if (calc-is-inverse)
|
|
73 (calc-vector-op "harm" 'calcFunc-vhmean arg)
|
|
74 (calc-vector-op "medn" 'calcFunc-vmedian arg))
|
|
75 (if (calc-is-inverse)
|
|
76 (calc-vector-op "meae" 'calcFunc-vmeane arg)
|
|
77 (calc-vector-op "mean" 'calcFunc-vmean arg))))
|
|
78 )
|
|
79
|
|
80 (defun calc-vector-mean-error (arg)
|
|
81 (interactive "P")
|
|
82 (calc-invert-func)
|
|
83 (calc-vector-mean arg)
|
|
84 )
|
|
85
|
|
86 (defun calc-vector-median (arg)
|
|
87 (interactive "P")
|
|
88 (calc-hyperbolic-func)
|
|
89 (calc-vector-mean arg)
|
|
90 )
|
|
91
|
|
92 (defun calc-vector-harmonic-mean (arg)
|
|
93 (interactive "P")
|
|
94 (calc-invert-func)
|
|
95 (calc-hyperbolic-func)
|
|
96 (calc-vector-mean arg)
|
|
97 )
|
|
98
|
|
99 (defun calc-vector-geometric-mean (arg)
|
|
100 (interactive "P")
|
|
101 (calc-slow-wrapper
|
|
102 (if (calc-is-hyperbolic)
|
|
103 (calc-binary-op "geom" 'calcFunc-agmean arg)
|
|
104 (calc-vector-op "geom" 'calcFunc-vgmean arg)))
|
|
105 )
|
|
106
|
|
107 (defun calc-vector-sdev (arg)
|
|
108 (interactive "P")
|
|
109 (calc-slow-wrapper
|
|
110 (if (calc-is-hyperbolic)
|
|
111 (if (calc-is-inverse)
|
|
112 (calc-vector-op "pvar" 'calcFunc-vpvar arg)
|
|
113 (calc-vector-op "var" 'calcFunc-vvar arg))
|
|
114 (if (calc-is-inverse)
|
|
115 (calc-vector-op "psdv" 'calcFunc-vpsdev arg)
|
|
116 (calc-vector-op "sdev" 'calcFunc-vsdev arg))))
|
|
117 )
|
|
118
|
|
119 (defun calc-vector-pop-sdev (arg)
|
|
120 (interactive "P")
|
|
121 (calc-invert-func)
|
|
122 (calc-vector-sdev arg)
|
|
123 )
|
|
124
|
|
125 (defun calc-vector-variance (arg)
|
|
126 (interactive "P")
|
|
127 (calc-hyperbolic-func)
|
|
128 (calc-vector-sdev arg)
|
|
129 )
|
|
130
|
|
131 (defun calc-vector-pop-variance (arg)
|
|
132 (interactive "P")
|
|
133 (calc-invert-func)
|
|
134 (calc-hyperbolic-func)
|
|
135 (calc-vector-sdev arg)
|
|
136 )
|
|
137
|
|
138 (defun calc-vector-covariance (arg)
|
|
139 (interactive "P")
|
|
140 (calc-slow-wrapper
|
|
141 (let ((n (if (eq arg 1) 1 2)))
|
|
142 (if (calc-is-hyperbolic)
|
|
143 (calc-enter-result n "corr" (cons 'calcFunc-vcorr
|
|
144 (calc-top-list-n n)))
|
|
145 (if (calc-is-inverse)
|
|
146 (calc-enter-result n "pcov" (cons 'calcFunc-vpcov
|
|
147 (calc-top-list-n n)))
|
|
148 (calc-enter-result n "cov" (cons 'calcFunc-vcov
|
|
149 (calc-top-list-n n)))))))
|
|
150 )
|
|
151
|
|
152 (defun calc-vector-pop-covariance (arg)
|
|
153 (interactive "P")
|
|
154 (calc-invert-func)
|
|
155 (calc-vector-covariance arg)
|
|
156 )
|
|
157
|
|
158 (defun calc-vector-correlation (arg)
|
|
159 (interactive "P")
|
|
160 (calc-hyperbolic-func)
|
|
161 (calc-vector-covariance arg)
|
|
162 )
|
|
163
|
|
164 (defun calc-vector-op (name func arg)
|
|
165 (setq calc-aborted-prefix name
|
|
166 arg (prefix-numeric-value arg))
|
|
167 (if (< arg 0)
|
|
168 (error "Negative arguments not allowed"))
|
|
169 (calc-enter-result arg name (cons func (calc-top-list-n arg)))
|
|
170 )
|
|
171
|
|
172
|
|
173
|
|
174
|
|
175 ;;; Useful statistical functions
|
|
176
|
|
177 ;;; Sum, product, etc., of one or more values or vectors.
|
|
178 ;;; Each argument must be either a number or a vector. Vectors
|
|
179 ;;; are flattened, but variables inside are assumed to represent
|
|
180 ;;; non-vectors.
|
|
181
|
|
182 (defun calcFunc-vsum (&rest vecs)
|
|
183 (math-reduce-many-vecs 'calcFunc-add 'calcFunc-vsum vecs 0)
|
|
184 )
|
|
185
|
|
186 (defun calcFunc-vprod (&rest vecs)
|
|
187 (math-reduce-many-vecs 'calcFunc-mul 'calcFunc-vprod vecs 1)
|
|
188 )
|
|
189
|
|
190 (defun calcFunc-vmax (&rest vecs)
|
|
191 (if (eq (car-safe (car vecs)) 'sdev)
|
|
192 '(var inf var-inf)
|
|
193 (if (eq (car-safe (car vecs)) 'intv)
|
|
194 (nth 3 (math-fix-int-intv (car vecs)))
|
|
195 (math-reduce-many-vecs 'calcFunc-max 'calcFunc-vmax vecs
|
|
196 '(neg (var inf var-inf)))))
|
|
197 )
|
|
198
|
|
199 (defun calcFunc-vmin (&rest vecs)
|
|
200 (if (eq (car-safe (car vecs)) 'sdev)
|
|
201 '(neg (var inf var-inf))
|
|
202 (if (eq (car-safe (car vecs)) 'intv)
|
|
203 (nth 2 (math-fix-int-intv (car vecs)))
|
|
204 (math-reduce-many-vecs 'calcFunc-min 'calcFunc-vmin vecs
|
|
205 '(var inf var-inf))))
|
|
206 )
|
|
207
|
|
208 (defun math-reduce-many-vecs (func whole-func vecs ident)
|
|
209 (let ((const-part nil)
|
|
210 (symb-part nil)
|
|
211 val vec)
|
|
212 (let ((calc-internal-prec (+ calc-internal-prec 2)))
|
|
213 (while vecs
|
|
214 (setq val (car vecs))
|
|
215 (and (eq (car-safe val) 'var)
|
|
216 (eq (car-safe (calc-var-value (nth 2 val))) 'vec)
|
|
217 (setq val (symbol-value (nth 2 val))))
|
|
218 (cond ((Math-vectorp val)
|
|
219 (setq vec (append (and const-part (list const-part))
|
|
220 (math-flatten-vector val)))
|
|
221 (setq const-part (if vec
|
|
222 (calcFunc-reducer
|
|
223 (math-calcFunc-to-var func)
|
|
224 (cons 'vec vec))
|
|
225 ident)))
|
|
226 ((or (Math-objectp val) (math-infinitep val))
|
|
227 (setq const-part (if const-part
|
|
228 (funcall func const-part val)
|
|
229 val)))
|
|
230 (t
|
|
231 (setq symb-part (nconc symb-part (list val)))))
|
|
232 (setq vecs (cdr vecs))))
|
|
233 (if const-part
|
|
234 (progn
|
|
235 (setq const-part (math-normalize const-part))
|
|
236 (if symb-part
|
|
237 (funcall func const-part (cons whole-func symb-part))
|
|
238 const-part))
|
|
239 (if symb-part (cons whole-func symb-part) ident)))
|
|
240 )
|
|
241
|
|
242
|
|
243 ;;; Return the number of data elements among the arguments.
|
|
244 (defun calcFunc-vcount (&rest vecs)
|
|
245 (let ((count 0))
|
|
246 (while vecs
|
|
247 (setq count (if (Math-vectorp (car vecs))
|
|
248 (+ count (math-count-elements (car vecs)))
|
|
249 (if (Math-objectp (car vecs))
|
|
250 (1+ count)
|
|
251 (if (and (eq (car-safe (car vecs)) 'var)
|
|
252 (eq (car-safe (calc-var-value
|
|
253 (nth 2 (car vecs))))
|
|
254 'vec))
|
|
255 (+ count (math-count-elements
|
|
256 (symbol-value (nth 2 (car vecs)))))
|
|
257 (math-reject-arg (car vecs) 'numvecp))))
|
|
258 vecs (cdr vecs)))
|
|
259 count)
|
|
260 )
|
|
261
|
|
262 (defun math-count-elements (vec)
|
|
263 (let ((count 0))
|
|
264 (while (setq vec (cdr vec))
|
|
265 (setq count (if (Math-vectorp (car vec))
|
|
266 (+ count (math-count-elements (car vec)))
|
|
267 (1+ count))))
|
|
268 count)
|
|
269 )
|
|
270
|
|
271
|
|
272 (defun math-flatten-many-vecs (vecs)
|
|
273 (let ((p vecs)
|
|
274 (vec (list 'vec)))
|
|
275 (while p
|
|
276 (setq vec (nconc vec
|
|
277 (if (Math-vectorp (car p))
|
|
278 (math-flatten-vector (car p))
|
|
279 (if (Math-objectp (car p))
|
|
280 (list (car p))
|
|
281 (if (and (eq (car-safe (car p)) 'var)
|
|
282 (eq (car-safe (calc-var-value
|
|
283 (nth 2 (car p)))) 'vec))
|
|
284 (math-flatten-vector (symbol-value
|
|
285 (nth 2 (car p))))
|
|
286 (math-reject-arg (car p) 'numvecp)))))
|
|
287 p (cdr p)))
|
|
288 vec)
|
|
289 )
|
|
290
|
|
291 (defun calcFunc-vflat (&rest vecs)
|
|
292 (math-flatten-many-vecs vecs)
|
|
293 )
|
|
294
|
|
295 (defun math-split-sdev-vec (vec zero-ok)
|
|
296 (let ((means (list 'vec))
|
|
297 (wts (list 'vec))
|
|
298 (exact nil)
|
|
299 (p vec))
|
|
300 (while (and (setq p (cdr p))
|
|
301 (not (and (consp (car p))
|
|
302 (eq (car (car p)) 'sdev)))))
|
|
303 (if (null p)
|
|
304 (list vec nil)
|
|
305 (while (setq vec (cdr vec))
|
|
306 (if (and (consp (setq p (car vec)))
|
|
307 (eq (car p) 'sdev))
|
|
308 (or exact
|
|
309 (setq means (cons (nth 1 p) means)
|
|
310 wts (cons (nth 2 p) wts)))
|
|
311 (if zero-ok
|
|
312 (setq means (cons (nth 1 p) means)
|
|
313 wts (cons 0 wts))
|
|
314 (or exact
|
|
315 (setq means (list 'vec)
|
|
316 wts nil
|
|
317 exact t))
|
|
318 (setq means (cons p means)))))
|
|
319 (list (nreverse means)
|
|
320 (and wts (nreverse wts)))))
|
|
321 )
|
|
322
|
|
323
|
|
324 ;;; Return the arithmetic mean of the argument numbers or vectors.
|
|
325 ;;; (If numbers are error forms, computes the weighted mean.)
|
|
326 (defun calcFunc-vmean (&rest vecs)
|
|
327 (let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
|
|
328 (means (car split))
|
|
329 (wts (nth 1 split))
|
|
330 (len (1- (length means))))
|
|
331 (if (= len 0)
|
|
332 (math-reject-arg nil "*Must be at least 1 argument")
|
|
333 (if (and (= len 1) (eq (car-safe (nth 1 means)) 'intv))
|
|
334 (let ((x (math-fix-int-intv (nth 1 means))))
|
|
335 (calcFunc-vmean (nth 2 x) (nth 3 x)))
|
|
336 (math-with-extra-prec 2
|
|
337 (if (and wts (> len 1))
|
|
338 (let* ((sqrwts (calcFunc-map '(var mul var-mul) wts wts))
|
|
339 (suminvsqrwts (calcFunc-reduce
|
|
340 '(var add var-add)
|
|
341 (calcFunc-map '(var div var-div)
|
|
342 1 sqrwts))))
|
|
343 (math-div (calcFunc-reduce '(var add var-add)
|
|
344 (calcFunc-map '(var div var-div)
|
|
345 means sqrwts))
|
|
346 suminvsqrwts))
|
|
347 (math-div (calcFunc-reduce '(var add var-add) means) len))))))
|
|
348 )
|
|
349
|
|
350 (defun math-fix-int-intv (x)
|
|
351 (if (math-floatp x)
|
|
352 x
|
|
353 (list 'intv 3
|
|
354 (if (memq (nth 1 x) '(2 3)) (nth 2 x) (math-add (nth 2 x) 1))
|
|
355 (if (memq (nth 1 x) '(1 3)) (nth 3 x) (math-sub (nth 3 x) 1))))
|
|
356 )
|
|
357
|
|
358 ;;; Compute the mean with an error estimate.
|
|
359 (defun calcFunc-vmeane (&rest vecs)
|
|
360 (let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
|
|
361 (means (car split))
|
|
362 (wts (nth 1 split))
|
|
363 (len (1- (length means))))
|
|
364 (if (= len 0)
|
|
365 (math-reject-arg nil "*Must be at least 1 argument")
|
|
366 (math-with-extra-prec 2
|
|
367 (if wts
|
|
368 (let* ((sqrwts (calcFunc-map '(var mul var-mul) wts wts))
|
|
369 (suminvsqrwts (calcFunc-reduce
|
|
370 '(var add var-add)
|
|
371 (calcFunc-map '(var div var-div)
|
|
372 1 sqrwts))))
|
|
373 (math-make-sdev
|
|
374 (math-div (calcFunc-reduce '(var add var-add)
|
|
375 (calcFunc-map '(var div var-div)
|
|
376 means sqrwts))
|
|
377 suminvsqrwts)
|
|
378 (list 'calcFunc-sqrt (math-div 1 suminvsqrwts))))
|
|
379 (let ((mean (math-div (calcFunc-reduce '(var add var-add) means)
|
|
380 len)))
|
|
381 (math-make-sdev
|
|
382 mean
|
|
383 (list 'calcFunc-sqrt
|
|
384 (math-div (calcFunc-reducer
|
|
385 '(var add var-add)
|
|
386 (calcFunc-map '(var pow var-pow)
|
|
387 (calcFunc-map '(var abs var-abs)
|
|
388 (calcFunc-map
|
|
389 '(var add var-add)
|
|
390 means
|
|
391 (math-neg mean)))
|
|
392 2))
|
|
393 (math-mul len (1- len))))))))))
|
|
394 )
|
|
395
|
|
396
|
|
397 ;;; Compute the median of a list of values.
|
|
398 (defun calcFunc-vmedian (&rest vecs)
|
|
399 (let* ((flat (copy-sequence (cdr (math-flatten-many-vecs vecs))))
|
|
400 (p flat)
|
|
401 (len (length flat))
|
|
402 (hlen (/ len 2)))
|
|
403 (if (= len 0)
|
|
404 (math-reject-arg nil "*Must be at least 1 argument")
|
|
405 (if (and (= len 1) (memq (car-safe (car flat)) '(sdev intv)))
|
|
406 (calcFunc-vmean (car flat))
|
|
407 (while p
|
|
408 (if (eq (car-safe (car p)) 'sdev)
|
|
409 (setcar p (nth 1 (car p))))
|
|
410 (or (Math-anglep (car p))
|
|
411 (math-reject-arg (car p) 'anglep))
|
|
412 (setq p (cdr p)))
|
|
413 (setq flat (sort flat 'math-lessp))
|
|
414 (if (= (% len 2) 0)
|
|
415 (math-div (math-add (nth (1- hlen) flat) (nth hlen flat)) 2)
|
|
416 (nth hlen flat)))))
|
|
417 )
|
|
418
|
|
419
|
|
420 (defun calcFunc-vgmean (&rest vecs)
|
|
421 (let* ((flat (math-flatten-many-vecs vecs))
|
|
422 (len (1- (length flat))))
|
|
423 (if (= len 0)
|
|
424 (math-reject-arg nil "*Must be at least 1 argument")
|
|
425 (math-with-extra-prec 2
|
|
426 (let ((x (calcFunc-reduce '(var mul math-mul) flat)))
|
|
427 (if (= len 2)
|
|
428 (math-sqrt x)
|
|
429 (math-pow x (list 'frac 1 len)))))))
|
|
430 )
|
|
431
|
|
432
|
|
433 (defun calcFunc-agmean (a b)
|
|
434 (cond ((Math-equal a b) a)
|
|
435 ((math-zerop a) a)
|
|
436 ((math-zerop b) b)
|
|
437 (calc-symbolic-mode (math-inexact-result))
|
|
438 ((not (Math-realp a)) (math-reject-arg a 'realp))
|
|
439 ((not (Math-realp b)) (math-reject-arg b 'realp))
|
|
440 (t
|
|
441 (math-with-extra-prec 2
|
|
442 (setq a (math-float (math-abs a))
|
|
443 b (math-float (math-abs b)))
|
|
444 (let (mean)
|
|
445 (while (not (math-nearly-equal-float a b))
|
|
446 (setq mean (math-mul-float (math-add-float a b) '(float 5 -1))
|
|
447 b (math-sqrt-float (math-mul-float a b))
|
|
448 a mean))
|
|
449 a))))
|
|
450 )
|
|
451
|
|
452
|
|
453 (defun calcFunc-vhmean (&rest vecs)
|
|
454 (let* ((flat (math-flatten-many-vecs vecs))
|
|
455 (len (1- (length flat))))
|
|
456 (if (= len 0)
|
|
457 (math-reject-arg nil "*Must be at least 1 argument")
|
|
458 (math-with-extra-prec 2
|
|
459 (math-div len
|
|
460 (calcFunc-reduce '(var add math-add)
|
|
461 (calcFunc-map '(var inv var-inv) flat))))))
|
|
462 )
|
|
463
|
|
464
|
|
465
|
|
466 ;;; Compute the sample variance or standard deviation of numbers or vectors.
|
|
467 ;;; (If the numbers are error forms, only the mean part of them is used.)
|
|
468 (defun calcFunc-vvar (&rest vecs)
|
|
469 (if (and (= (length vecs) 1)
|
|
470 (memq (car-safe (car vecs)) '(sdev intv)))
|
|
471 (if (eq (car-safe (car vecs)) 'intv)
|
|
472 (math-intv-variance (car vecs) nil)
|
|
473 (math-sqr (nth 2 (car vecs))))
|
|
474 (math-covariance vecs nil nil 0))
|
|
475 )
|
|
476
|
|
477 (defun calcFunc-vsdev (&rest vecs)
|
|
478 (if (and (= (length vecs) 1)
|
|
479 (memq (car-safe (car vecs)) '(sdev intv)))
|
|
480 (if (eq (car-safe (car vecs)) 'intv)
|
|
481 (if (math-floatp (car vecs))
|
|
482 (math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
|
|
483 (math-sqrt-12))
|
|
484 (math-sqrt (calcFunc-vvar (car vecs))))
|
|
485 (nth 2 (car vecs)))
|
|
486 (math-sqrt (math-covariance vecs nil nil 0)))
|
|
487 )
|
|
488
|
|
489 ;;; Compute the population variance or std deviation of numbers or vectors.
|
|
490 (defun calcFunc-vpvar (&rest vecs)
|
|
491 (if (and (= (length vecs) 1)
|
|
492 (memq (car-safe (car vecs)) '(sdev intv)))
|
|
493 (if (eq (car-safe (car vecs)) 'intv)
|
|
494 (math-intv-variance (car vecs) t)
|
|
495 (math-sqr (nth 2 (car vecs))))
|
|
496 (math-covariance vecs nil t 0))
|
|
497 )
|
|
498
|
|
499 (defun calcFunc-vpsdev (&rest vecs)
|
|
500 (if (and (= (length vecs) 1)
|
|
501 (memq (car-safe (car vecs)) '(sdev intv)))
|
|
502 (if (eq (car-safe (car vecs)) 'intv)
|
|
503 (if (math-floatp (car vecs))
|
|
504 (math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
|
|
505 (math-sqrt-12))
|
|
506 (math-sqrt (calcFunc-vpvar (car vecs))))
|
|
507 (nth 2 (car vecs)))
|
|
508 (math-sqrt (math-covariance vecs nil t 0)))
|
|
509 )
|
|
510
|
|
511 (defun math-intv-variance (x pop)
|
|
512 (or (math-constp x) (math-reject-arg x 'constp))
|
|
513 (if (math-floatp x)
|
|
514 (math-div (math-sqr (math-sub (nth 3 x) (nth 2 x))) 12)
|
|
515 (let* ((x (math-fix-int-intv x))
|
|
516 (len (math-sub (nth 3 x) (nth 2 x)))
|
|
517 (hlen (math-quotient len 2)))
|
|
518 (math-div (if (math-evenp len)
|
|
519 (calcFunc-sum '(^ (var X var-X) 2) '(var X var-X)
|
|
520 (math-neg hlen) hlen)
|
|
521 (calcFunc-sum '(^ (- (var X var-X) (/ 1 2)) 2)
|
|
522 '(var X var-X)
|
|
523 (math-neg hlen) (math-add hlen 1)))
|
|
524 (if pop (math-add len 1) len))))
|
|
525 )
|
|
526
|
|
527 ;;; Compute the covariance and linear correlation coefficient.
|
|
528 (defun calcFunc-vcov (vec1 &optional vec2)
|
|
529 (math-covariance (list vec1) (list vec2) nil 1)
|
|
530 )
|
|
531
|
|
532 (defun calcFunc-vpcov (vec1 &optional vec2)
|
|
533 (math-covariance (list vec1) (list vec2) t 1)
|
|
534 )
|
|
535
|
|
536 (defun calcFunc-vcorr (vec1 &optional vec2)
|
|
537 (math-covariance (list vec1) (list vec2) nil 2)
|
|
538 )
|
|
539
|
|
540
|
|
541 (defun math-covariance (vec1 vec2 pop mode)
|
|
542 (or (car vec2) (= mode 0)
|
|
543 (progn
|
|
544 (if (and (eq (car-safe (car vec1)) 'var)
|
|
545 (eq (car-safe (calc-var-value (nth 2 (car vec1)))) 'vec))
|
|
546 (setq vec1 (symbol-value (nth 2 (car vec1))))
|
|
547 (setq vec1 (car vec1)))
|
|
548 (or (math-matrixp vec1) (math-dimension-error))
|
|
549 (or (= (length (nth 1 vec1)) 3) (math-dimension-error))
|
|
550 (setq vec2 (list (math-mat-col vec1 2))
|
|
551 vec1 (list (math-mat-col vec1 1)))))
|
|
552 (math-with-extra-prec 2
|
|
553 (let* ((split1 (math-split-sdev-vec (math-flatten-many-vecs vec1) nil))
|
|
554 (means1 (car split1))
|
|
555 (wts1 (nth 1 split1))
|
|
556 split2 means2 (wts2 nil)
|
|
557 (sqrwts nil)
|
|
558 suminvsqrwts
|
|
559 (len (1- (length means1))))
|
|
560 (if (< len (if pop 1 2))
|
|
561 (math-reject-arg nil (if pop
|
|
562 "*Must be at least 1 argument"
|
|
563 "*Must be at least 2 arguments")))
|
|
564 (if (or wts1 wts2)
|
|
565 (setq sqrwts (math-add
|
|
566 (if wts1
|
|
567 (calcFunc-map '(var mul var-mul) wts1 wts1)
|
|
568 0)
|
|
569 (if wts2
|
|
570 (calcFunc-map '(var mul var-mul) wts2 wts2)
|
|
571 0))
|
|
572 suminvsqrwts (calcFunc-reduce
|
|
573 '(var add var-add)
|
|
574 (calcFunc-map '(var div var-div) 1 sqrwts))))
|
|
575 (or (= mode 0)
|
|
576 (progn
|
|
577 (setq split2 (math-split-sdev-vec (math-flatten-many-vecs vec2)
|
|
578 nil)
|
|
579 means2 (car split2)
|
|
580 wts2 (nth 2 split1))
|
|
581 (or (= len (1- (length means2))) (math-dimension-error))))
|
|
582 (let* ((diff1 (calcFunc-map
|
|
583 '(var add var-add)
|
|
584 means1
|
|
585 (if sqrwts
|
|
586 (math-div (calcFunc-reduce
|
|
587 '(var add var-add)
|
|
588 (calcFunc-map '(var div var-div)
|
|
589 means1 sqrwts))
|
|
590 (math-neg suminvsqrwts))
|
|
591 (math-div (calcFunc-reducer '(var add var-add) means1)
|
|
592 (- len)))))
|
|
593 (diff2 (if (= mode 0)
|
|
594 diff1
|
|
595 (calcFunc-map
|
|
596 '(var add var-add)
|
|
597 means2
|
|
598 (if sqrwts
|
|
599 (math-div (calcFunc-reduce
|
|
600 '(var add var-add)
|
|
601 (calcFunc-map '(var div var-div)
|
|
602 means2 sqrwts))
|
|
603 (math-neg suminvsqrwts))
|
|
604 (math-div (calcFunc-reducer '(var add var-add) means2)
|
|
605 (- len))))))
|
|
606 (covar (calcFunc-map '(var mul var-mul) diff1 diff2)))
|
|
607 (if sqrwts
|
|
608 (setq covar (calcFunc-map '(var div var-div) covar sqrwts)))
|
|
609 (math-div
|
|
610 (calcFunc-reducer '(var add var-add) covar)
|
|
611 (if (= mode 2)
|
|
612 (let ((var1 (calcFunc-map '(var mul var-mul) diff1 diff1))
|
|
613 (var2 (calcFunc-map '(var mul var-mul) diff2 diff2)))
|
|
614 (if sqrwts
|
|
615 (setq var1 (calcFunc-map '(var div var-div) var1 sqrwts)
|
|
616 var2 (calcFunc-map '(var div var-div) var2 sqrwts)))
|
|
617 (math-sqrt
|
|
618 (math-mul (calcFunc-reducer '(var add var-add) var1)
|
|
619 (calcFunc-reducer '(var add var-add) var2))))
|
|
620 (if sqrwts
|
|
621 (if pop
|
|
622 suminvsqrwts
|
|
623 (math-div (math-mul suminvsqrwts (1- len)) len))
|
|
624 (if pop len (1- len))))))))
|
|
625 )
|
|
626
|
|
627
|
|
628
|
|
629
|