Update copyright year to 2014 by running admin/update-copyright.
[emacs.git] / lisp / calc / calc-stat.el
blob2b63821a6ba55a602c213c59a1f1d47e56a70d80
1 ;;; calc-stat.el --- statistical functions for Calc
3 ;; Copyright (C) 1990-1993, 2001-2014 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is free software: you can redistribute it and/or modify
11 ;; it under the terms of the GNU General Public License as published by
12 ;; the Free Software Foundation, either version 3 of the License, or
13 ;; (at your option) any later version.
15 ;; GNU Emacs is distributed in the hope that it will be useful,
16 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ;; GNU General Public License for more details.
20 ;; You should have received a copy of the GNU General Public License
21 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
23 ;;; Commentary:
25 ;;; Code:
27 ;; This file is autoloaded from calc-ext.el.
29 (require 'calc-ext)
30 (require 'calc-macs)
32 ;;; Statistical operations on vectors.
34 (defun calc-vector-count (arg)
35 (interactive "P")
36 (calc-slow-wrapper
37 (calc-vector-op "coun" 'calcFunc-vcount arg)))
39 (defun calc-vector-sum (arg)
40 (interactive "P")
41 (calc-slow-wrapper
42 (if (calc-is-hyperbolic)
43 (calc-vector-op "vprd" 'calcFunc-vprod arg)
44 (calc-vector-op "vsum" 'calcFunc-vsum arg))))
46 (defun calc-vector-product (arg)
47 (interactive "P")
48 (calc-hyperbolic-func)
49 (calc-vector-sum arg))
51 (defun calc-vector-max (arg)
52 (interactive "P")
53 (calc-slow-wrapper
54 (if (calc-is-inverse)
55 (calc-vector-op "vmin" 'calcFunc-vmin arg)
56 (calc-vector-op "vmax" 'calcFunc-vmax arg))))
58 (defun calc-vector-min (arg)
59 (interactive "P")
60 (calc-invert-func)
61 (calc-vector-max arg))
63 (defun calc-vector-mean (arg)
64 (interactive "P")
65 (calc-slow-wrapper
66 (if (calc-is-hyperbolic)
67 (if (calc-is-inverse)
68 (calc-vector-op "harm" 'calcFunc-vhmean arg)
69 (calc-vector-op "medn" 'calcFunc-vmedian arg))
70 (if (calc-is-inverse)
71 (calc-vector-op "meae" 'calcFunc-vmeane arg)
72 (calc-vector-op "mean" 'calcFunc-vmean arg)))))
74 (defun calc-vector-mean-error (arg)
75 (interactive "P")
76 (calc-invert-func)
77 (calc-vector-mean arg))
79 (defun calc-vector-median (arg)
80 (interactive "P")
81 (calc-hyperbolic-func)
82 (calc-vector-mean arg))
84 (defun calc-vector-harmonic-mean (arg)
85 (interactive "P")
86 (calc-invert-func)
87 (calc-hyperbolic-func)
88 (calc-vector-mean arg))
90 (defun calc-vector-geometric-mean (arg)
91 (interactive "P")
92 (calc-slow-wrapper
93 (if (calc-is-hyperbolic)
94 (calc-binary-op "geom" 'calcFunc-agmean arg)
95 (calc-vector-op "geom" 'calcFunc-vgmean arg))))
97 (defun calc-vector-sdev (arg)
98 (interactive "P")
99 (calc-slow-wrapper
100 (if (calc-is-hyperbolic)
101 (if (calc-is-inverse)
102 (calc-vector-op "pvar" 'calcFunc-vpvar arg)
103 (calc-vector-op "var" 'calcFunc-vvar arg))
104 (if (calc-is-inverse)
105 (calc-vector-op "psdv" 'calcFunc-vpsdev arg)
106 (calc-vector-op "sdev" 'calcFunc-vsdev arg)))))
108 (defun calc-vector-pop-sdev (arg)
109 (interactive "P")
110 (calc-invert-func)
111 (calc-vector-sdev arg))
113 (defun calc-vector-variance (arg)
114 (interactive "P")
115 (calc-hyperbolic-func)
116 (calc-vector-sdev arg))
118 (defun calc-vector-pop-variance (arg)
119 (interactive "P")
120 (calc-invert-func)
121 (calc-hyperbolic-func)
122 (calc-vector-sdev arg))
124 (defun calc-vector-covariance (arg)
125 (interactive "P")
126 (calc-slow-wrapper
127 (let ((n (if (eq arg 1) 1 2)))
128 (if (calc-is-hyperbolic)
129 (calc-enter-result n "corr" (cons 'calcFunc-vcorr
130 (calc-top-list-n n)))
131 (if (calc-is-inverse)
132 (calc-enter-result n "pcov" (cons 'calcFunc-vpcov
133 (calc-top-list-n n)))
134 (calc-enter-result n "cov" (cons 'calcFunc-vcov
135 (calc-top-list-n n))))))))
137 (defun calc-vector-pop-covariance (arg)
138 (interactive "P")
139 (calc-invert-func)
140 (calc-vector-covariance arg))
142 (defun calc-vector-correlation (arg)
143 (interactive "P")
144 (calc-hyperbolic-func)
145 (calc-vector-covariance arg))
147 (defun calc-vector-op (name func arg)
148 (setq calc-aborted-prefix name
149 arg (prefix-numeric-value arg))
150 (if (< arg 0)
151 (error "Negative arguments not allowed"))
152 (calc-enter-result arg name (cons func (calc-top-list-n arg))))
157 ;;; Useful statistical functions
159 ;;; Sum, product, etc., of one or more values or vectors.
160 ;;; Each argument must be either a number or a vector. Vectors
161 ;;; are flattened, but variables inside are assumed to represent
162 ;;; non-vectors.
164 (defun calcFunc-vsum (&rest vecs)
165 (math-reduce-many-vecs 'calcFunc-add 'calcFunc-vsum vecs 0))
167 (defun calcFunc-vprod (&rest vecs)
168 (math-reduce-many-vecs 'calcFunc-mul 'calcFunc-vprod vecs 1))
170 (defun calcFunc-vmax (&rest vecs)
171 (if (eq (car-safe (car vecs)) 'sdev)
172 '(var inf var-inf)
173 (if (eq (car-safe (car vecs)) 'intv)
174 (nth 3 (math-fix-int-intv (car vecs)))
175 (math-reduce-many-vecs 'calcFunc-max 'calcFunc-vmax vecs
176 '(neg (var inf var-inf))))))
178 (defun calcFunc-vmin (&rest vecs)
179 (if (eq (car-safe (car vecs)) 'sdev)
180 '(neg (var inf var-inf))
181 (if (eq (car-safe (car vecs)) 'intv)
182 (nth 2 (math-fix-int-intv (car vecs)))
183 (math-reduce-many-vecs 'calcFunc-min 'calcFunc-vmin vecs
184 '(var inf var-inf)))))
186 (defun math-reduce-many-vecs (func whole-func vecs ident)
187 (let ((const-part nil)
188 (symb-part nil)
189 val vec)
190 (let ((calc-internal-prec (+ calc-internal-prec 2)))
191 (while vecs
192 (setq val (car vecs))
193 (and (eq (car-safe val) 'var)
194 (eq (car-safe (calc-var-value (nth 2 val))) 'vec)
195 (setq val (symbol-value (nth 2 val))))
196 (cond ((Math-vectorp val)
197 (setq vec (append (and const-part (list const-part))
198 (math-flatten-vector val)))
199 (setq const-part (if vec
200 (calcFunc-reducer
201 (math-calcFunc-to-var func)
202 (cons 'vec vec))
203 ident)))
204 ((or (Math-objectp val) (math-infinitep val))
205 (setq const-part (if const-part
206 (funcall func const-part val)
207 val)))
209 (setq symb-part (nconc symb-part (list val)))))
210 (setq vecs (cdr vecs))))
211 (if const-part
212 (progn
213 (setq const-part (math-normalize const-part))
214 (if symb-part
215 (funcall func const-part (cons whole-func symb-part))
216 const-part))
217 (if symb-part (cons whole-func symb-part) ident))))
220 ;;; Return the number of data elements among the arguments.
221 (defun calcFunc-vcount (&rest vecs)
222 (let ((count 0))
223 (while vecs
224 (setq count (if (Math-vectorp (car vecs))
225 (+ count (math-count-elements (car vecs)))
226 (if (Math-objectp (car vecs))
227 (1+ count)
228 (if (and (eq (car-safe (car vecs)) 'var)
229 (eq (car-safe (calc-var-value
230 (nth 2 (car vecs))))
231 'vec))
232 (+ count (math-count-elements
233 (symbol-value (nth 2 (car vecs)))))
234 (math-reject-arg (car vecs) 'numvecp))))
235 vecs (cdr vecs)))
236 count))
238 (defun math-count-elements (vec)
239 (let ((count 0))
240 (while (setq vec (cdr vec))
241 (setq count (if (Math-vectorp (car vec))
242 (+ count (math-count-elements (car vec)))
243 (1+ count))))
244 count))
247 (defun math-flatten-many-vecs (vecs)
248 (let ((p vecs)
249 (vec (list 'vec)))
250 (while p
251 (setq vec (nconc vec
252 (if (Math-vectorp (car p))
253 (math-flatten-vector (car p))
254 (if (Math-objectp (car p))
255 (list (car p))
256 (if (and (eq (car-safe (car p)) 'var)
257 (eq (car-safe (calc-var-value
258 (nth 2 (car p)))) 'vec))
259 (math-flatten-vector (symbol-value
260 (nth 2 (car p))))
261 (math-reject-arg (car p) 'numvecp)))))
262 p (cdr p)))
263 vec))
265 (defun calcFunc-vflat (&rest vecs)
266 (math-flatten-many-vecs vecs))
268 (defun math-split-sdev-vec (vec zero-ok)
269 (let ((means (list 'vec))
270 (wts (list 'vec))
271 (exact nil)
272 (p vec))
273 (while (and (setq p (cdr p))
274 (not (and (consp (car p))
275 (eq (car (car p)) 'sdev)))))
276 (if (null p)
277 (list vec nil)
278 (while (setq vec (cdr vec))
279 (if (and (consp (setq p (car vec)))
280 (eq (car p) 'sdev))
281 (or exact
282 (setq means (cons (nth 1 p) means)
283 wts (cons (nth 2 p) wts)))
284 (if zero-ok
285 (setq means (cons (nth 1 p) means)
286 wts (cons 0 wts))
287 (or exact
288 (setq means (list 'vec)
289 wts nil
290 exact t))
291 (setq means (cons p means)))))
292 (list (nreverse means)
293 (and wts (nreverse wts))))))
296 ;;; Return the arithmetic mean of the argument numbers or vectors.
297 ;;; (If numbers are error forms, computes the weighted mean.)
298 (defun calcFunc-vmean (&rest vecs)
299 (let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
300 (means (car split))
301 (wts (nth 1 split))
302 (len (1- (length means))))
303 (if (= len 0)
304 (math-reject-arg nil "*Must be at least 1 argument")
305 (if (and (= len 1) (eq (car-safe (nth 1 means)) 'intv))
306 (let ((x (math-fix-int-intv (nth 1 means))))
307 (calcFunc-vmean (nth 2 x) (nth 3 x)))
308 (math-with-extra-prec 2
309 (if (and wts (> len 1))
310 (let* ((sqrwts (calcFunc-map '(var mul var-mul) wts wts))
311 (suminvsqrwts (calcFunc-reduce
312 '(var add var-add)
313 (calcFunc-map '(var div var-div)
314 1 sqrwts))))
315 (math-div (calcFunc-reduce '(var add var-add)
316 (calcFunc-map '(var div var-div)
317 means sqrwts))
318 suminvsqrwts))
319 (math-div (calcFunc-reduce '(var add var-add) means) len)))))))
321 (defun math-fix-int-intv (x)
322 (if (math-floatp x)
324 (list 'intv 3
325 (if (memq (nth 1 x) '(2 3)) (nth 2 x) (math-add (nth 2 x) 1))
326 (if (memq (nth 1 x) '(1 3)) (nth 3 x) (math-sub (nth 3 x) 1)))))
328 ;;; Compute the mean with an error estimate.
329 (defun calcFunc-vmeane (&rest vecs)
330 (let* ((split (math-split-sdev-vec (math-flatten-many-vecs vecs) nil))
331 (means (car split))
332 (wts (nth 1 split))
333 (len (1- (length means))))
334 (if (= len 0)
335 (math-reject-arg nil "*Must be at least 1 argument")
336 (math-with-extra-prec 2
337 (if wts
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-make-sdev
344 (math-div (calcFunc-reduce '(var add var-add)
345 (calcFunc-map '(var div var-div)
346 means sqrwts))
347 suminvsqrwts)
348 (list 'calcFunc-sqrt (math-div 1 suminvsqrwts))))
349 (let ((mean (math-div (calcFunc-reduce '(var add var-add) means)
350 len)))
351 (math-make-sdev
352 mean
353 (list 'calcFunc-sqrt
354 (math-div (calcFunc-reducer
355 '(var add var-add)
356 (calcFunc-map '(var pow var-pow)
357 (calcFunc-map '(var abs var-abs)
358 (calcFunc-map
359 '(var add var-add)
360 means
361 (math-neg mean)))
363 (math-mul len (1- len)))))))))))
366 ;;; Compute the median of a list of values.
367 (defun calcFunc-vmedian (&rest vecs)
368 (let* ((flat (copy-sequence (cdr (math-flatten-many-vecs vecs))))
369 (p flat)
370 (len (length flat))
371 (hlen (/ len 2)))
372 (if (= len 0)
373 (math-reject-arg nil "*Must be at least 1 argument")
374 (if (and (= len 1) (memq (car-safe (car flat)) '(sdev intv)))
375 (calcFunc-vmean (car flat))
376 (while p
377 (if (eq (car-safe (car p)) 'sdev)
378 (setcar p (nth 1 (car p))))
379 (or (Math-anglep (car p))
380 (math-reject-arg (car p) 'anglep))
381 (setq p (cdr p)))
382 (setq flat (sort flat 'math-lessp))
383 (if (= (% len 2) 0)
384 (math-div (math-add (nth (1- hlen) flat) (nth hlen flat)) 2)
385 (nth hlen flat))))))
388 (defun calcFunc-vgmean (&rest vecs)
389 (let* ((flat (math-flatten-many-vecs vecs))
390 (len (1- (length flat))))
391 (if (= len 0)
392 (math-reject-arg nil "*Must be at least 1 argument")
393 (math-with-extra-prec 2
394 (let ((x (calcFunc-reduce '(var mul math-mul) flat)))
395 (if (= len 2)
396 (math-sqrt x)
397 (math-pow x (list 'frac 1 len))))))))
400 (defun calcFunc-agmean (a b)
401 (cond ((Math-equal a b) a)
402 ((math-zerop a) a)
403 ((math-zerop b) b)
404 (calc-symbolic-mode (math-inexact-result))
405 ((not (Math-realp a)) (math-reject-arg a 'realp))
406 ((not (Math-realp b)) (math-reject-arg b 'realp))
408 (math-with-extra-prec 2
409 (setq a (math-float (math-abs a))
410 b (math-float (math-abs b)))
411 (let (mean)
412 (while (not (math-nearly-equal-float a b))
413 (setq mean (math-mul-float (math-add-float a b) '(float 5 -1))
414 b (math-sqrt-float (math-mul-float a b))
415 a mean))
416 a)))))
419 (defun calcFunc-vhmean (&rest vecs)
420 (let* ((flat (math-flatten-many-vecs vecs))
421 (len (1- (length flat))))
422 (if (= len 0)
423 (math-reject-arg nil "*Must be at least 1 argument")
424 (math-with-extra-prec 2
425 (math-div len
426 (calcFunc-reduce '(var add math-add)
427 (calcFunc-map '(var inv var-inv) flat)))))))
431 ;;; Compute the sample variance or standard deviation of numbers or vectors.
432 ;;; (If the numbers are error forms, only the mean part of them is used.)
433 (defun calcFunc-vvar (&rest vecs)
434 (if (and (= (length vecs) 1)
435 (memq (car-safe (car vecs)) '(sdev intv)))
436 (if (eq (car-safe (car vecs)) 'intv)
437 (math-intv-variance (car vecs) nil)
438 (math-sqr (nth 2 (car vecs))))
439 (math-covariance vecs nil nil 0)))
441 (defun calcFunc-vsdev (&rest vecs)
442 (if (and (= (length vecs) 1)
443 (memq (car-safe (car vecs)) '(sdev intv)))
444 (if (eq (car-safe (car vecs)) 'intv)
445 (if (math-floatp (car vecs))
446 (math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
447 (math-sqrt-12))
448 (math-sqrt (calcFunc-vvar (car vecs))))
449 (nth 2 (car vecs)))
450 (math-sqrt (math-covariance vecs nil nil 0))))
452 ;;; Compute the population variance or std deviation of numbers or vectors.
453 (defun calcFunc-vpvar (&rest vecs)
454 (if (and (= (length vecs) 1)
455 (memq (car-safe (car vecs)) '(sdev intv)))
456 (if (eq (car-safe (car vecs)) 'intv)
457 (math-intv-variance (car vecs) t)
458 (math-sqr (nth 2 (car vecs))))
459 (math-covariance vecs nil t 0)))
461 (defun calcFunc-vpsdev (&rest vecs)
462 (if (and (= (length vecs) 1)
463 (memq (car-safe (car vecs)) '(sdev intv)))
464 (if (eq (car-safe (car vecs)) 'intv)
465 (if (math-floatp (car vecs))
466 (math-div (math-sub (nth 3 (car vecs)) (nth 2 (car vecs)))
467 (math-sqrt-12))
468 (math-sqrt (calcFunc-vpvar (car vecs))))
469 (nth 2 (car vecs)))
470 (math-sqrt (math-covariance vecs nil t 0))))
472 (defun math-intv-variance (x pop)
473 (or (math-constp x) (math-reject-arg x 'constp))
474 (if (math-floatp x)
475 (math-div (math-sqr (math-sub (nth 3 x) (nth 2 x))) 12)
476 (let* ((x (math-fix-int-intv x))
477 (len (math-sub (nth 3 x) (nth 2 x)))
478 (hlen (math-quotient len 2)))
479 (math-div (if (math-evenp len)
480 (calcFunc-sum '(^ (var X var-X) 2) '(var X var-X)
481 (math-neg hlen) hlen)
482 (calcFunc-sum '(^ (- (var X var-X) (/ 1 2)) 2)
483 '(var X var-X)
484 (math-neg hlen) (math-add hlen 1)))
485 (if pop (math-add len 1) len)))))
487 ;;; Compute the covariance and linear correlation coefficient.
488 (defun calcFunc-vcov (vec1 &optional vec2)
489 (math-covariance (list vec1) (list vec2) nil 1))
491 (defun calcFunc-vpcov (vec1 &optional vec2)
492 (math-covariance (list vec1) (list vec2) t 1))
494 (defun calcFunc-vcorr (vec1 &optional vec2)
495 (math-covariance (list vec1) (list vec2) nil 2))
498 (defun math-covariance (vec1 vec2 pop mode)
499 (or (car vec2) (= mode 0)
500 (progn
501 (if (and (eq (car-safe (car vec1)) 'var)
502 (eq (car-safe (calc-var-value (nth 2 (car vec1)))) 'vec))
503 (setq vec1 (symbol-value (nth 2 (car vec1))))
504 (setq vec1 (car vec1)))
505 (or (math-matrixp vec1) (math-dimension-error))
506 (or (= (length (nth 1 vec1)) 3) (math-dimension-error))
507 (setq vec2 (list (math-mat-col vec1 2))
508 vec1 (list (math-mat-col vec1 1)))))
509 (math-with-extra-prec 2
510 (let* ((split1 (math-split-sdev-vec (math-flatten-many-vecs vec1) nil))
511 (means1 (car split1))
512 (wts1 (nth 1 split1))
513 split2 means2 (wts2 nil)
514 (sqrwts nil)
515 suminvsqrwts
516 (len (1- (length means1))))
517 (if (< len (if pop 1 2))
518 (math-reject-arg nil (if pop
519 "*Must be at least 1 argument"
520 "*Must be at least 2 arguments")))
521 (if (or wts1 wts2)
522 (setq sqrwts (math-add
523 (if wts1
524 (calcFunc-map '(var mul var-mul) wts1 wts1)
526 (if wts2
527 (calcFunc-map '(var mul var-mul) wts2 wts2)
529 suminvsqrwts (calcFunc-reduce
530 '(var add var-add)
531 (calcFunc-map '(var div var-div) 1 sqrwts))))
532 (or (= mode 0)
533 (progn
534 (setq split2 (math-split-sdev-vec (math-flatten-many-vecs vec2)
535 nil)
536 means2 (car split2)
537 wts2 (nth 2 split1))
538 (or (= len (1- (length means2))) (math-dimension-error))))
539 (let* ((diff1 (calcFunc-map
540 '(var add var-add)
541 means1
542 (if sqrwts
543 (math-div (calcFunc-reduce
544 '(var add var-add)
545 (calcFunc-map '(var div var-div)
546 means1 sqrwts))
547 (math-neg suminvsqrwts))
548 (math-div (calcFunc-reducer '(var add var-add) means1)
549 (- len)))))
550 (diff2 (if (= mode 0)
551 diff1
552 (calcFunc-map
553 '(var add var-add)
554 means2
555 (if sqrwts
556 (math-div (calcFunc-reduce
557 '(var add var-add)
558 (calcFunc-map '(var div var-div)
559 means2 sqrwts))
560 (math-neg suminvsqrwts))
561 (math-div (calcFunc-reducer '(var add var-add) means2)
562 (- len))))))
563 (covar (calcFunc-map '(var mul var-mul) diff1 diff2)))
564 (if sqrwts
565 (setq covar (calcFunc-map '(var div var-div) covar sqrwts)))
566 (math-div
567 (calcFunc-reducer '(var add var-add) covar)
568 (if (= mode 2)
569 (let ((var1 (calcFunc-map '(var mul var-mul) diff1 diff1))
570 (var2 (calcFunc-map '(var mul var-mul) diff2 diff2)))
571 (if sqrwts
572 (setq var1 (calcFunc-map '(var div var-div) var1 sqrwts)
573 var2 (calcFunc-map '(var div var-div) var2 sqrwts)))
574 (math-sqrt
575 (math-mul (calcFunc-reducer '(var add var-add) var1)
576 (calcFunc-reducer '(var add var-add) var2))))
577 (if sqrwts
578 (if pop
579 suminvsqrwts
580 (math-div (math-mul suminvsqrwts (1- len)) len))
581 (if pop len (1- len)))))))))
583 (provide 'calc-stat)
585 ;;; calc-stat.el ends here