; An example of running XLISPSTAT and finding the Laplace approximation
; for the marginal of $\theta_0$.
; xlispstat 

(def Y  (list 65  156 100 134 16 108 121 4 39 143 56 26 22 1 1 5 65))


(def wbc (list 2300 750 4300 2600 6000 10500 10000 17000 5400 7000 9400 
	 32000 35000 100000 100000 52000 100000))

(def  x (log (/ wbc 10000)))

(defun loglik (theta)
  (let* ((theta0 (first theta))
	 (theta1 (second theta))
	 (t1x (* theta1 x)))
    (- (sum t1x) 
       (* (length x) (log theta0))
       (/ (sum (* y (exp t1x)))
	  theta0))))


(defun logprior (theta)
   (let ((theta0 (first theta))
	 (theta1 (second theta)))
     (+ (log (* (/ 1 5) (exp (- (/ theta1 5)))))
	(log (* (/ 1 theta0) 
		(exp (- (/ (^ (/ (- (log theta0) (log 52)) (log 2))  2) 2))))
	     ))
     )
)

(defun logpost (theta) (+ (loglik theta) (logprior theta)))

(def post.obj (bayes-model #'logpost (list (exp 3.5) .8)))

(send post.obj :moments)  


(send post.obj :help :margin1)


(plot-lines (send post.obj :margin1 0 (rseq 1 100 100)))


; list xy pairs to copy into a file and read into S

(def xy (transpose (send post.obj :margin1 0 (rseq 1 100 100))))
(dotimes (i 30) (apply #'format t "~g ~g ~%" (select xy i)))