;; gibbs.lsp

(defun norm-rand (mu sigma)
  (+ (* (first (normal-rand 1)) sigma) mu)
  )

(defun gibbs (n rho)
  (let (
	(x 0)
	(y 0)
	(sig (sqrt (- 1 (^ rho 2))))
	(xvec (iseq n))
	(yvec (iseq n))
	)
    (setf (select xvec 0) x)
    (setf (select yvec 0) y)
    (dolist (i (iseq 1 (- n 1)))
      (setf x (norm-rand (* rho y) sig))
      (setf y (norm-rand (* rho x) sig))
      (setf (select xvec i) x)
      (setf (select yvec i) y)
      )
    (list xvec yvec)
    )
  )

(def bvn (gibbs 10000 0.98))
(plot-points (select bvn 0) (select bvn 1))
(plot-lines (select bvn 0) (select bvn 1))
(histogram (select bvn 0))
(histogram (select bvn 1))






;; end