;;;;Peter Santiago, Started June 16th 2008
;;;;compiled under:
;;;;slime 1:20070927-3
;;;;SBCL  1:1.0.11.0-1
;;;;Wilcoxon algorithm and information shamelessly stolen from Richard Lowry:
;;;;http://faculty.vassar.edu/lowry/ch12a.html

(load "utests.lisp") ;loads unit-testing file
  
(deftest wilcoxon-test ()
  "runs wilcoxon on data with known results."
  (check
   (= 1 (wilcoxon '(72 90 40 84 22 78 50 64 50 30 52 64 45 64 24 78)
                  '(32 58 20 68 36 68 40 56 44 25 56 68 48 62 24 78) 95)
            )))

(defun wilcoxon (x.a x.b &optional (signum 95))
  "Runs a wilcoxon signed-rank test"
  (labels ((combine (a &rest elements) `(,a ,@ elements))
           (as-ranks (l r) (mapcar #'(lambda (x) (gethash x r)) l))
           (simple-quad (c) (/ (+ -1 (sqrt (- 1 (* -1 8 c)))) 2)))
    
    (let* ((x.dif (mapcar #'- x.a x.b)) ;x.a-x.b
           (x.abs (mapcar #'abs x.dif)) ;absolute value of x.dif
           (x (sort (mapcar #'combine x.abs x.dif x.b x.a)
                    #'< :key #'first)) ;sorts by absolute value
           (ranks (as-ranks (mapcar #'first x)
                            (rank (mapcar #'first x)))) ;gets the ranks
           (x (mapcar #'cons ranks x)) ;adds absolute ranks to x
           (x (mapcar #'cons
                      (mapcar #'*
                              (mapcar #'first x)
                              (mapcar #'signum
                                      (mapcar #'third x))) x))
                                        ;adds signed ranks
           (n (simple-quad
               (loop for item in ranks
                     summing item))) ;number of ranks used
           (w (loop for item in (mapcar #'first x)
                    summing item)) ;total weight of all signed ranks
           (z (/ (- w .5) (sqrt (/ (* n (+ n 1) (+ (* 2 n) 1)) 6))))) ;z-value
      
      (test-significance signum n w z))))

(defun test-significance (signum total.ranks weight z)
  "Tests for significance: 0 if not, 1 if"
  (let ((val 0))
    (cond ((< total.ranks 5))
          
          ((< total.ranks 10)
           (< weight (get-wcrit signum total.ranks))
           (if (> weight (get-wcrit signum total.ranks))
               (setf val 1)))
          
          ((< 10 total.ranks)
           (< z (get-zcrit signum))
           (if (> z (get-zcrit signum))
               (setf val 1))))
    val))

(defun get-zcrit (signum)
  "for tests with greater than 10 ranks, use this table"
  (let ((table '((90   1.645)
                 (95   1.960)
                 (98   2.326)
                 (99   2.576)
                 (99.9 3.291))))
    (dotimes (i (list-length table))
      (if (= (first (nth i table)) signum)
          (return (second (nth i table)))))))

(defun get-wcrit (signum total.ranks)
  "for tests run with 5-9 ranks, you must use this table"
  (let* ((inf most-positive-fixnum)
         (table `((90   (5 15)    (6 17)    (7 22)   (8 26)  (9 29))
                  (95   (5 ,inf)  (6 21)    (7 24)   (8 30)  (9 35))
                  (98   (5 ,inf)  (6 ,inf)  (7 28)   (8 34)  (9 39))
                  (99   (5 ,inf)  (6 ,inf)  (7 ,inf) (8 36)  (9 43))
                  (99.9 (5 ,inf)  (6 ,inf)  (7 ,inf) (8 ,inf)(9 ,inf)))))
    (dotimes (i (list-length table))
      (if (= (first (nth i table)) signum)
          (dotimes (j 5)
            (if (= total.ranks (first (nth (+ j 1) (nth i table))))
                (return-from get-wcrit
                  (second (nth (+ j 1) (nth i table))))))))))

(defun rank (l &optional (ranks (make-hash-table)) (n 0))
  "Returns a hash of the ranks in a sorted list.  All numbers in a run of
   repeated entries get the average rank of that run."
  (if (null l)
      ranks
      (let (repeats sum now)
	(labels ((walk () (incf n) (pop l))
		 (new  () (setf repeats 1) (setf sum n))
		 (same () (incf sum n) (incf repeats))
		 (spin () (when (eql now (car l))
			    (walk) (same) (spin))))
	  (setf now (walk))
	  (new)
	  (spin)
          (cond ((= now 0)
                 (setf repeats 1 sum 0 now 0 n 0)))
	  (setf (gethash now ranks) (/ sum repeats))
	  (rank l ranks n)))))