Have you ever thought about adding Welch's t-test as a primitive to newLISP? Could be another "mode" of the current t-test primitive or you could come up with another name.
Here's a newLISP implementation -- for fun, but also I had to implement it for a stat analysis I recently performed because I couldn't assume that the two population variances were equal, and people better than I at this inform me that this is the better test under this condition.
Code: Select all
(define (welch-t-test xs ys)
(letn (nx (length xs) ny (length ys)
nx-1 (- nx 1) ny-1 (- ny 1)
mx (mean xs) my (mean ys)
s2x (s2 xs) s2y (s2 ys)
s2x/nx (div s2x nx) s2y/ny (div s2y ny)
df (div (sq (add s2x/nx s2y/ny))
(add (div (sq s2x/nx) nx-1)
(div (sq s2y/ny) ny-1)))
t (div (sub mx my)
(sqrt (add s2x/nx s2y/ny))))
(list mx my (ssd xs) (ssd ys)
t df (mul 2 (prob-t (abs t) df)))))
Code: Select all
;; Arithmetic mean
(define (mean xs) (div (apply add xs) (length xs)))
;; Square function
(define (sq x) (mul x x))
;; Sum of squares
(define (sumsq xs) (apply add (map sq xs)))
;; The sample variance, just the numerator.
(define (s2-num xs)
(let (n (length xs))
(sub (sumsq xs) (div (sq (apply add xs)) n))))
;; Sample variance
(define (s2 xs) (div (s2-num xs) (- (length xs) 1)))
;; Sample standard deviation
(define (ssd xs) (sqrt (s2 xs)))
Code: Select all
> (t-test '(10 4 7 1 1 6 1 8 2 4) '(4 6 9 4 6 8 9 3 9 7))
(4.4 6.5 3.238655414 2.273030283 -1.678359739 18 0.1105533782)
> (welch-t-test '(10 4 7 1 1 6 1 8 2 4) '(4 6 9 4 6 8 9 3 9 7))
(4.4 6.5 3.238655414 2.273030283 -1.678359739 16.13523413 0.1126979729)
> ;; Now, try unequal sample sizes (which will yield different t values).
> (t-test '(10 4 7 1 1 6 1 8 2 4) '(4 6 9 4 6 8 9 3))
(4.4 6.125 3.238655414 2.356601669 -1.260037546 16 0.2257247063)
> (welch-t-test '(10 4 7 1 1 6 1 8 2 4) '(4 6 9 4 6 8 9 3))
(4.4 6.125 3.238655414 2.356601669 -1.30656126 15.90049882 0.2110406063)