; Ranking functions required by mwu and wilcoxon ranking tests.
(defun rank (data0)
  (let* ((starter "someCraZYsymBOL")
	 (old starter)
	 (n (length data0))
	 (data (sort data0 #'<))
	 (start 1)
	 (sum 0)
	 skipping
	 r
	 (ranks (make-hash-table)))
    (do ((i 1 (1+ i))) ((> i n))
      (setf skipping (or
		      (equal old starter)
		      (equal (nth i data) old)))
      (if skipping
	  (incf sum i)
	  (progn
	    (setf r (/ sum (- i start)))
	    (do ((j start (1+ j))) ((= j i))
	      (setf (gethash (nth j data) ranks) r))
	    (setf start i
		  sum i)))
      (setf old (nth i data)))
    (if skipping
	(setf (gethash (nth n data) ranks) (/ sum (- n start)))
	(if (equal (gethash (nth (- n 1) data) ranks "NotHere") "NotHere")
	    (setf (gethash (nth (- n 1) data) ranks) (+ r 1))))
    ranks))

; Slight bug somewhere. will fix soon.

;;; Wilcoxon statistical test - How it works.
;;; http://faculty.vassar.edu/lowry/ch12a.html
(defun wilcoxon (population1 population2 critical)
  (let* ((n 0)
	 (w 0)
	 sigma
	 z
	 ranks
	 (diff (make-list (length population1) :initial-element 0))
	 (absDiff (make-list (length population1) :initial-element 0)))
    (dotimes (i (length population1))
      (let* ((delta (- (nth i population1) (nth i population2))))
	(print delta)
	(if (not (= delta 0))
	    (progn
	      (incf n)
	      (setf (nth i diff) delta
		    (nth i absDiff) (abs delta))))))
  (setf ranks (rank absDiff))
  (dotimes (i (1- (length absDiff)))
    (let* ((w0 (gethash (nth i absDiff) ranks 0)))
      (print w0)
      (incf w (if (< (nth i diff) 0)
		  (* -1 w0)
		  w0))))
  (print "n")
  (print n)
  (setf sigma (sqrt (/ (* n
			  (+ n 1)
			  (+ (* 2 n) 1))
		       6)))
  (setf z (if (= (floor sigma) 0.0)
	      (/ (- w 0.5)
		 sigma)
	      0))
  (if (and 
       (>= z 0)
       (<= z (critical-value critical)))
      0
      1)))

;; Mann-Whitney Statistical test - How it works:
;; http://faculty.vassar.edu/lowry/ch11a.html
(defun mann-whitney(population1 population2 critical)
  (let* ((data (make-list (+ (length population1) (length population2))))
	 (ranks1 (make-list (length population1)))
	 (ranks2 (make-list (length population2)))
	 (n 0)
	 (n1 0)
	 (n2 0)
	 (sum1 0)
	 (sum2 0))
    (dolist (i population1)
      (setf (nth n data) i)
      (incf n))
    (dolist (i population2)
      (setf (nth n data) i)
      (incf n))
    (let* ((ranks (rank data)))
      (dotimes (i (length population1))
	(progn
	  (incf n1)
	  (setf (nth i ranks1) (gethash (nth i population1) ranks 0.001))
	  (incf sum1 (nth i ranks1))))
      (dotimes (i (length population2))
	(progn
	  (incf n2)
	  (setf (nth i ranks2) (gethash (nth i population2) ranks 0.001))
	  (incf sum2 (nth i ranks2))))
      (let* ((meanU (/ (* n1
			  (+ n1 n2 1))
		       2))
	     (sdU (expt (/ (* n1 n2
			      (+ n1 n2 1))
			   12)
			0.5))
	     (correction (if (> sum1 meanU) -0.5 0.5))
	     (z (abs (/ (+ (- sum1
			      meanU)
			   correction)
			sdU))))
	(if (and 
	     (>= z 0)
	     (<= z (critical-value critical)))
	    0
	    1)))))


	
(defun test-wilcoxon ()
  (let* ((population (list 1 4.6 2 4.7 3 4.9 4 5.1 5 5.2 6 5.5 7 5.8 8 6.1 9 6.5 10 6.5 11 7.2))
	 (population2 (map 'list (lambda(x) (* 2 x)) population)))
    (wilcoxon population population 95)))
  

(defun test-mwu ()
    (let* ((population (list 1 4.6 2 4.7 3 4.9 4 5.1 5 5.2 6 5.5 7 5.8 8 6.1 9 6.5 10 6.5 11 7.2))
	   (population2 (map 'list (lambda(x) (* 0.5 x)) population)))
      (mann-whitney population population2 95)))


; Picks a critical value for the test. the default case is for 0.95
(defun critical-value (conf)
  (cond 
    ((= conf 99) 2.326)
    ((= conf 90) 1.646)
    (t 1.960)))