(defun ndup (e n)
  (if (= n 0)    
    nil      
    (cons e (ndup e (- n 1)))))

(defun sum (l)
  (if (eq l nil)
    0
    (+ (first l) (sum (rest l)))))

(defun mean (l)
  (/ (sum l) (length l)))

(defun square (x)
  (* x x))

(defun variance (l)
  (/ (sum (mapcar #'square (mapcar #'- l (ndup (mean l) (length l))))) (- (length l) 1)))

(defun sd (l)
  (sqrt (variance l)))

(defun rand ()
  (float (/ (random 100000000) 100000000)))

(defun bernoulli (p)
  (if (< (rand) p) 1 0))

(defun binomial (n p)
  (if (= n 0) 0 (+ (bernoulli p) (binomial (- n 1) p))))

(defun uniform (_min _max)
  (+ _min (* (rand) (- _max _min))))

(defun exponential (rate)
  (/ (* -1 (log (rand))) rate))

(defun accept_reject ()
  (let* ((Y1 (exponential 1))
	 (Y2 (exponential 1))
	 (Z (if (< (rand) 0.5) Y1 (* -1 Y1))))
    (if (> (- Y2 (/ (expt (- Y1 1) 2) 2)) 0) Z (accept_reject))))

(defun normal (mu sigma)
  (+ mu (* (accept_reject) sigma)))