(defun mann-whitney (a b &optional (conf 95))
  "The 'method 1' of http://faculty.vassar.edu/lowry/ch11a.html."
  (labels ((as-ranks (l r) (mapcar #'(lambda (x) (gethash x r)) l))
	   (sum      (l)   (let ((s 0)) (dolist (x l s)(incf s x))))
	   (median   (l n) (let ((sorted (sort l #'<))
				 (midv (floor (/ n 2))))
			     (if (oddp n)
				 (nth midv sorted)
				 (/ (+ (nth (1- midv) sorted) 
				       (nth  midv     sorted)) 2)))))
    (let* ((all    (sort (copy-list (append a  b)) #'<))
	   (ranks  (rank all))
	   (ranksa (as-ranks a ranks))
	   (ranksb (as-ranks b ranks))
	   (na     (length a))
	   (nb     (length b))
	   (n      (+ na nb))	   
           (tcrict (tcritical n conf))           
	   (suma   (* 1.0 (sum ranksa)))
           (ta     (/ (* na (+ n 1)) 2.0))
	   (sigma  (sqrt (/  (* na nb (+ n 1)) 12.0)))
	   (za     (/ (+ (- suma ta) 0.5) sigma)))
                                        ;     (o suma ta sigma tcrict) ;debug line
      (cond  ((< (abs za) tcrict)                        0)
	     ((< (median ranksa na) (median ranksb nb)) -1)
	     (t                                          1))
      )))

(defun mann-whitney-demo-small ()
  (mann-whitney '(4.6   4.7   4.9  5.1   5.2   5.5  
		  5.8   6.1   6.5  6.5   7.2)
		'(5.2   5.3   5.4  5.6   6.2   6.3  
		  6.8   7.7   8.0  8.1)))
	
(defun mann-whitney-demo-big (&optional (fudge 1)) 
  "Generate two lists of 10000 random floats. Multiple the
  second  list by 'fudge'. Check if the median ranks of list one
  is the same as list two (return '0'), smaller than list two
  (return '-1'), or larger than list two (return '1')."
  (labels ((big (n s) (let (out) 
			(dotimes (i n out) 
			  (push (* s (round (random 100))) out)))))
    (let ((one (big 10000 1))
	  (two (big 10000 fudge)))
      (mann-whitney one two))))

(deftest test-mann-whitney ()
  (check
    (= 0 (mann-whitney-demo-small))
    (= 0 (mann-whitney-demo-big))
    (= 1 (mann-whitney-demo-big -2))
    (= -1 (mann-whitney-demo-big 2))
  ))