home *** CD-ROM | disk | FTP | other *** search
/ Fish 'n' More 2 / fishmore-publicdomainlibraryvol.ii1991xetec.iso / fish / applications / xlispstat / lisp.lzh / XLisp-Stat / bayes.lsp < prev    next >
Text File  |  1990-10-03  |  27KB  |  658 lines
  1. ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
  2. ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
  3. ;;;; You may give out copies of this software; for conditions see the file
  4. ;;;; COPYING included with this distribution.
  5.  
  6. (require "maximize")
  7. (provide "bayes")
  8.  
  9. ;;;;
  10. ;;;; Objects Representing Functions
  11. ;;;;
  12.  
  13. ;; Generic C2 Functions
  14.  
  15. (defproto c2-function-proto '(f h num-derivs))
  16.  
  17. (defmeth c2-function-proto :isnew (f &optional (h .001) (num-derivs 0))
  18.   (setf (slot-value 'f) f)
  19.   (setf (slot-value 'h) (if (numberp h) (list h h) h))
  20.   (setf (slot-value 'num-derivs) num-derivs))
  21.  
  22. (defmeth c2-function-proto :f (&optional f)
  23.   (if f (setf (slot-value 'f) f))
  24.   (slot-value 'f))
  25.  
  26. (defmeth c2-function-proto :grad-h () (first (slot-value 'h)))
  27. (defmeth c2-function-proto :hess-h () (second (slot-value 'h)))
  28. (defmeth c2-function-proto :num-derivs () (slot-value 'num-derivs))
  29.  
  30. (defmeth c2-function-proto :value (x)
  31.   (let ((f (send self :f)))
  32.     (if (objectp f) 
  33.         (send f :value x) 
  34.         (let ((v (funcall f x)))
  35.           (if (consp v) (first v) v)))))
  36.  
  37. (defmeth c2-function-proto :gradient (x &optional (h (send self :grad-h)))
  38.   (let ((f (send self :f)))
  39.     (if (objectp f) (send f :gradient x h) (numgrad f x nil h))))
  40.  
  41. (defmeth c2-function-proto :hessian (x &optional (h (send self :hess-h)))
  42.   (let ((f (send self :f)))
  43.     (if (objectp f) (send f :hessian x h) (numhess f x nil h))))
  44.  
  45. (defmeth c2-function-proto :vals (x &optional (h (send self :hess-h)))
  46.   (let ((f (send self :f)))
  47.     (if (objectp f)
  48.         (send f :vals x h)
  49.         (let ((v (funcall f x)))
  50.           (if (consp v)
  51.               (if (= (length v) 3)
  52.                   v  
  53.                   (list (first v) (second v) (send self :hessian x h)))
  54.               (list v (send self :gradient x h) (send self :hessian x h)))))))
  55.  
  56. (defmeth c2-function-proto :vals (x &optional (h (send self :hess-h)))
  57.   (let ((f (send self :f)))
  58.     (if (objectp f) (send f :vals x h) (numhess f x nil h t))))
  59.  
  60.  
  61. ;; Scaled C2 Functions
  62.  
  63. (defproto scaled-c2-function-proto '(scaling) () c2-function-proto)
  64.  
  65. ;;**** allow function objects?
  66. (defmeth scaled-c2-function-proto :isnew (f &optional 
  67.                                             theta
  68.                                             sigma
  69.                                             (center 0)
  70.                                             (scale 1)
  71.                                             (h 0.001))
  72.   (let* ((value (funcall f theta))
  73.          (num-derivs (if (consp value) (- (length value) 1) -1))
  74.          (sigma-t (if (< 0 num-derivs) (transpose sigma))))
  75.     (labels ((scale (v)
  76.                (if v
  77.                    (case num-derivs
  78.                      (-1 (/ (- v center) scale))
  79.                      (0 (/ (- (first v) center) scale))
  80.                      (1 (list (/ (- (first v) center) scale)
  81.                               (matmult sigma-t (/ (second v) scale))))
  82.                      (2 (list (/ (- (first v) center) scale)
  83.                               (matmult sigma-t (/ (second v) scale))
  84.                               (matmult sigma-t (/ (third v) scale) sigma))))))
  85.              (sf (x) (scale (funcall f (ax+y sigma x theta t)))))
  86.       (call-next-method #'sf h num-derivs))))
  87.  
  88. ;; Tilted C2 Functions
  89. ;; **** allow nil values?
  90. (defproto tilt-function-proto '(tilt exptilt) () c2-function-proto)
  91.  
  92. (defmeth tilt-function-proto :isnew (&optional f (tilt .1) (h .001))
  93.   (call-next-method f h)
  94.   (setf (slot-value 'exptilt) t)
  95.   (setf (slot-value 'tilt) tilt))
  96.  
  97. (defmeth tilt-function-proto :tilt (&optional tilt)
  98.   (if tilt (setf (slot-value 'tilt) tilt))
  99.   (slot-value 'tilt))
  100.  
  101. (defmeth tilt-function-proto :exptilt (&optional (new nil set))
  102.   (if set (setf (slot-value 'exptilt) new))
  103.   (slot-value 'exptilt))
  104.  
  105. (defmeth tilt-function-proto :value (x)
  106.   (let ((f (send self :f))
  107.         (tilt (send self :tilt))
  108.         (exptilt (send self :exptilt)))
  109.     (flet ((value (f)
  110.              (let ((v (send f :value x)))
  111.                (if exptilt v (log v)))))
  112.       (* tilt (if (consp f) (apply #'+ (mapcar #'value f)) (value f))))))
  113.  
  114. (defmeth tilt-function-proto :gradient (x &optional (h (send self :grad-h)))
  115.   (let ((f (send self :f))
  116.         (tilt (send self :tilt))
  117.         (exptilt (send self :exptilt)))
  118.     (flet ((gradient (f)
  119.              (if exptilt
  120.                  (send f :gradient x h)
  121.                  (let ((v (send f :value x))
  122.                        (grad (send f :gradient x h)))
  123.                    (/ grad v)))))
  124.       (* tilt 
  125.          (if (consp f) (apply #'+ (mapcar #'gradient f)) (gradient f))))))
  126.  
  127. (defmeth tilt-function-proto :hessian (x &optional (h (send self :hess-h)))
  128.   (let ((f (send self :f))
  129.         (tilt (send self :tilt))
  130.         (exptilt (send self :exptilt)))
  131.     (flet ((hessian (f)
  132.              (let* ((vals (send f :vals x h))
  133.                     (v (first vals))
  134.                     (grad (if exptilt (second vals) (/ (second vals) v)))
  135.                     (hess (if exptilt (third vals) (/ (third vals) v))))
  136.                (if exptilt hess (- hess (outer-product grad grad))))))
  137.       (* tilt (if (consp f) (apply #'+ (mapcar #'hessian f)) (hessian f))))))
  138.  
  139. (defmeth tilt-function-proto :vals (x &optional (h (send self :hess-h)))
  140.   (let ((f (send self :f))
  141.         (tilt (send self :tilt))
  142.         (exptilt (send self :exptilt)))
  143.     (flet ((vals (f)
  144.              (let ((vals (send f :vals x h)))
  145.                (if exptilt
  146.                    vals
  147.                    (let* ((v (first vals))
  148.                           (grad (/ (second vals) v))
  149.                           (hess (- (/ (third vals) v) 
  150.                                    (outer-product grad grad))))
  151.                      (list (log v) grad hess))))))
  152.       (let ((v (if (consp f) (mapcar #'vals f) (vals f))))
  153.         (* tilt (if (consp f) (apply #'+ v) v))))))
  154.         
  155. ;; scaled log posterior prototype
  156.  
  157. (defproto scaled-logpost-proto 
  158.   '(tilt-object init-pars) () scaled-c2-function-proto)
  159.  
  160. (defmeth scaled-logpost-proto :isnew (f &optional 
  161.                                         theta sigma
  162.                                         (center 0) (scale 1) (h .001))
  163.   (let* ((n (length theta))
  164.          (m (repeat 0 n))
  165.          (m-grad (repeat 0 n))
  166.          (m-hess (- (identity-matrix n)))
  167.          (pars (list m m-grad m-hess)))
  168.     (call-next-method f theta sigma center scale h)
  169.     (setf (slot-value 'init-pars) pars)
  170.     (setf (slot-value 'tilt-object) (send tilt-function-proto :new))))
  171.  
  172. (defmeth scaled-logpost-proto :log-laplace (g &optional
  173.                                               (count-limit 2) det-only (h .1))
  174.   (let* ((x (send self :tilt-newton g count-limit))
  175.          (vals (send self :vals x h))
  176.          (gvals (if g (send g :vals x h)))
  177.          (hess (if g (+ (third vals) (third gvals)) (third vals)))
  178.          (det (- (sum (log (diagonal (first (chol-decomp (- hess)))))))))
  179.     (if det-only 
  180.         det 
  181.         (if g (+ (first vals) (first gvals) det) (+ (first vals) det)))))
  182.  
  183. (defmeth scaled-logpost-proto :tilt-newton (tilt &optional (count-limit 2))
  184.   (let* ((pars (slot-value 'init-pars))
  185.          (mode (first pars))
  186.          (mode-grad (second pars))
  187.          (mode-hess (third pars)))
  188.     (flet ((gradhess (x initial) 
  189.              (let ((gh (if (and initial mode-grad mode-hess)
  190.                                  (list mode-grad mode-hess)
  191.                                  (rest (send self :vals x)))))
  192.                (if tilt (+ gh (rest (send tilt :vals x))) gh)))
  193.            (newton-step (x gh) (- x (solve (second gh) (first gh)))))
  194.       (do* ((count 1 (+ count 1))
  195.             (gradhess (gradhess mode t) (gradhess x nil))
  196.             (x (newton-step mode gradhess) (newton-step x gradhess)))
  197.            ((>= count count-limit) x)))))
  198.  
  199. (defmeth scaled-logpost-proto :tilt-laplace (g tilt &optional 
  200.                                                (exptilt t) maxiter det-only h)
  201.   (let ((tilt-object (slot-value 'tilt-object)))
  202.     (send tilt-object :exptilt exptilt)
  203.     (send tilt-object :f g)
  204.     (send tilt-object :tilt tilt)
  205.     (send self :log-laplace tilt-object maxiter det-only h)))
  206.  
  207. (defmeth scaled-logpost-proto :tilt-mode (g tilt &key (exptilt t) (maxiter 2))
  208.   (let ((tilt-object (slot-value 'tilt-object)))
  209.     (send tilt-object :exptilt exptilt)
  210.     (send tilt-object :f g)
  211.     (send tilt-object :tilt tilt)
  212.     (send self :tilt-newton tilt-object maxiter)))
  213.  
  214. ;;;;
  215. ;;;; Bayes Model Prototype
  216. ;;;;
  217.  
  218. (defproto bayes-model-proto '(bayes-internals))
  219.  
  220. ;; initialization methods and constructor function
  221.  
  222. (defmeth bayes-model-proto :isnew (logpost mode &key
  223.                                            scale 
  224.                                            (derivstep .001)
  225.                                            (verbose t)
  226.                                            (maximize t)
  227.                                            domain)
  228.   (send self :set-bayes-internals 
  229.         logpost mode scale derivstep nil nil t domain)
  230.   (if maximize (send self :maximize verbose)))
  231.  
  232. (defun bayes-model (logpost mode &rest args &key (quick t) (print t))
  233. "Args: (logpost mode &key scale derivstep (verbose t)
  234.        (quick t) (print t)))
  235. LOGPOST computes the logposterior density. It should return the
  236. function, or a list of the function value and gradient, or a list of
  237. the function value, gradient and Hessian. MODE is an initial guess for
  238. the mode. SCALE and DERIVSTEP are used for numerical derivatives and
  239. scaling. VERBOSE controls printing of iteration information during
  240. optimization, PRINT controls printing of summary information. If QUICK
  241. is T the summary is based on first order approximations."
  242.   (let ((m (apply #'send bayes-model-proto :new logpost mode args)))
  243.     (if print (send m :display :quick quick))
  244.     m))
  245.  
  246. ;; display method
  247.  
  248. (defmeth bayes-model-proto :display (&key (quick t))
  249.   (let* ((moments (send self (if quick :1stmoments :moments)))
  250.          (means (first moments))
  251.          (stdevs (second moments))
  252.          (p-names (send self :parameter-names)))
  253.     (if quick
  254.         (format t "~2%First Order Approximations to Posterior Moments:~2%")
  255.         (format t "~2%Approximate Posterior Moments:~2%"))
  256.     (mapcar #'(lambda (name mu sd)
  257.                 (format t "~22a  ~10g (~a)~%" name mu sd))
  258.             p-names
  259.             means
  260.             stdevs)
  261.     (format t "~%")))
  262.  
  263. (defmeth bayes-model-proto :parameter-names ()
  264.   (let ((n (length (send self :mode))))
  265.     (mapcar #'(lambda (x) (format nil "Parameter ~d" x)) (iseq 0 (- n 1)))))
  266.  
  267. ;; implementation-dependent access methods
  268.  
  269. (defmeth bayes-model-proto :set-bayes-internals (lp m s h mval ch max dom)
  270.   (setf (slot-value 'bayes-internals) 
  271.         (vector lp m s h mval ch max dom)))
  272.  
  273. (defmeth bayes-model-proto :logpost (&optional new)
  274.   (let ((internals (slot-value 'bayes-internals)))
  275.     (when new
  276.           (setf (select internals 0) new)
  277.           (send self :needs-maximizing t))
  278.     (select internals 0)))
  279.  
  280. (defmeth bayes-model-proto :domain (&optional new)
  281.   (let ((internals (slot-value 'bayes-internals)))
  282.     (if new (setf (select internals 7) new))
  283.     (select internals 7)))
  284.  
  285. (defmeth bayes-model-proto :mode-values (&optional mode mval ch)
  286.   (let ((internals (slot-value 'bayes-internals)))
  287.     (when mode
  288.           (setf (select internals 1) mode)
  289.           (setf (select internals 4) mval)
  290.           (setf (select internals 5) ch))
  291.     (list (select internals 1)
  292.           (select internals 4)
  293.           (select internals 5))))
  294.  
  295. (defmeth bayes-model-proto :parameter-scale (&optional new)
  296.   (let ((internals (slot-value 'bayes-internals)))
  297.     (if new (setf (select internals 2) new))
  298.     (select internals 2)))
  299.  
  300. (defmeth bayes-model-proto :parameter-dimension ()
  301.   (length (select (slot-value 'bayes-internals) 1)))
  302.  
  303. (defmeth bayes-model-proto :derivstep ()
  304.   (select (slot-value 'bayes-internals) 3))
  305.  
  306. (defmeth bayes-model-proto :needs-maximizing (&optional (new nil set))
  307.   (let ((internals (slot-value 'bayes-internals)))
  308.     (if set (setf (select internals 6) new))
  309.     (select internals 6)))
  310.  
  311. ;; Transformation-Related Methods
  312. ;; (These should be the only ones needing to be changed to handle 
  313. ;; an internal parameter transformation; perhaps also :logpost)
  314.  
  315. ;; **** fix to be more careful about use of functionp
  316. (defun function-list (g &optional n)
  317.   (cond
  318.     ((or (functionp g) (objectp g)) (list g))
  319.     ((integerp g) 
  320.      (if (null n)
  321.          (list #'(lambda (x) (elt x g)))
  322.          (let ((grad (make-array n :initial-element 0))
  323.                (hess (make-array (list n n) :initial-element 0)))
  324.            (setf (aref grad g) 1)
  325.            (list #'(lambda (x) (list (elt x g) grad hess))))))
  326.     (t (mapcar #'(lambda (x) (car (function-list x n))) g))))
  327.  
  328. (defmeth bayes-model-proto :mode ()
  329.   (if (send self :needs-maximizing) (send self :maximize))
  330.   (first (send self :mode-values)))
  331.  
  332. (defmeth bayes-model-proto :new-mode-guess (new)
  333.   (send self :needs-maximizing t)
  334.   (send self :mode-values new))
  335.  
  336. (defmeth bayes-model-proto :transformed-logpost ()
  337.   (if (send self :needs-maximizing) (send self :maximize))
  338.   (let* ((m-values (send self :mode-values))
  339.          (mode (first m-values))
  340.          (mval (second m-values))
  341.          (ch (third m-values))
  342.          (h (send self :derivstep))
  343.          (f (send self :logpost)))
  344.     (send scaled-logpost-proto :new f mode ch mval 1 h)))
  345.  
  346. ;;**** need transformed domain here
  347.  
  348. (defmeth bayes-model-proto :transformed-functions (&optional g (c 0) (s 1))
  349.   (if (send self :needs-maximizing) (send self :maximize))
  350.   (let* ((m-values (send self :mode-values))
  351.          (mode (first m-values))
  352.          (mval (second m-values))
  353.          (ch (third m-values))
  354.          (h (send self :derivstep))
  355.          (n (length mode))
  356.          (g (function-list (if g g (iseq n)) n))
  357.          (c (if (numberp c) (repeat c (length g)) c))
  358.          (s (if (numberp s) (repeat s (length g)) s)))
  359.     (mapcar #'(lambda (g c s) 
  360.                 (send scaled-c2-function-proto :new g mode ch c s h))
  361.             g c s)))
  362.  
  363. ;; computing methods
  364.  
  365. (defmeth bayes-model-proto :maximize (&optional (verbose 0))
  366.   (let* ((lp (send self :logpost))
  367.          (x (first (send self :mode-values)))
  368.          (scale (send self :parameter-scale))
  369.          (h (send self :derivstep))
  370.          (minfo (newtonmax lp x 
  371.                            :scale scale
  372.                            :derivstep h
  373.                            :verbose verbose
  374.                            :return-derivs t))
  375.          (mode (first minfo))
  376.          (mval (second minfo))
  377.          (ch (first (chol-decomp (inverse (- (fourth minfo)))))))
  378.     (send self :mode-values mode mval ch)
  379.     (send self :needs-maximizing nil)
  380.     (send self :check-derivatives verbose)))
  381.  
  382. (defmeth bayes-model-proto :check-derivatives (&optional 
  383.                                                (verbose 0) 
  384.                                                (epsilon .00001))
  385.   (let* ((verbose (if (numberp verbose) (< 0 verbose) verbose))
  386.          (n (send self :parameter-dimension))
  387.          (tlp (send self :transformed-logpost))
  388.          (hess (send tlp :hessian (repeat 0 n)))
  389.          (needs-max (send self :needs-maximizing)))
  390.     (when (> (max (abs (+ hess (identity-matrix n)))) epsilon)
  391.           (if verbose (format t "Adjusting derivatives...~%"))
  392.           (let* ((ch (first (chol-decomp (- (inverse hess)))))
  393.                  (mvals (send self :mode-values))
  394.                  (m (matmult (third mvals) ch)))
  395.             (send self :mode-values (first mvals) (second mvals) m)
  396.             (if (not needs-max) (send self :needs-maximizing nil))
  397.             (if verbose
  398.                 (let* ((tlp (send self :transformed-logpost))
  399.                        (hess (send tlp :hessian (repeat 0 n))))
  400.                   (if (> (max (abs (+ hess (identity-matrix n)))) epsilon)
  401.                       (format t 
  402.                               "Derivatives may not be well-behaved.~%"))))))))
  403.  
  404. ;; moments
  405.  
  406. (defmeth bayes-model-proto :1stmoments (&optional gfuns &key covar)
  407. "Args: (&optional gfuns &key covar) 
  408. Computes first order approximations to posterior moments. GFUNS can be
  409. a parameter index, list of indices, a function of the parameters or a
  410. list of such functions. Returns a the list of first order approximate
  411. means and standard deviations if COVAR is NIL. If COVAR is T the
  412. covaraince is appended to the end of the result as well."
  413.   (if (send self :needs-maximizing) (send self :maximize))
  414.   (let* ((n (send self :parameter-dimension))
  415.          (x (repeat 0 n))
  416.          (g (send self :transformed-functions gfuns 0 1))
  417.          (grads (apply #'bind-columns 
  418.                        (mapcar #'(lambda (g) (send g :gradient x)) g)))
  419.          (mean (mapcar #'(lambda (g) (send g :value x)) g))
  420.          (cov (matmult  (transpose grads) grads)))
  421.     (if covar
  422.         (list mean (sqrt (diagonal cov)) cov)
  423.         (list mean (sqrt (diagonal cov))))))
  424.  
  425. (defmeth bayes-model-proto :mgfmoments (&optional g &key covar 
  426.                                                   (mgfdel .1)
  427.                                                   ((:derivstep h) .1)
  428.                                                   (maxiter 2))
  429.   (let* ((moms1 (send self :1stmoments g :covar covar))
  430.          (mean1 (first moms1))
  431.          (stdev1 (second moms1))
  432.          (cov1 (if covar (third moms1)))
  433.          (l-object (send self :transformed-logpost))
  434.          (g-objects (send self :transformed-functions g mean1 stdev1))
  435.          (ldet0 (send l-object :log-laplace nil maxiter t h)))
  436.     (labels ((lapdet (g tilt)
  437.                (- (send l-object :tilt-laplace g tilt t maxiter t h) ldet0))
  438.              (moms2 (m s g)
  439.                (let ((ldet1 (lapdet g mgfdel))
  440.                      (ldet2 (lapdet g (- mgfdel))))
  441.                  (list (+ m (* s (/  (- ldet1 ldet2) (* 2 mgfdel))))
  442.                        (* s (sqrt (+ 1 (/ (+ ldet1 ldet2) (^ mgfdel 2))))))))
  443.              (covar (g mean-sd)
  444.                (let* ((mu (first mean-sd))
  445.                       (sd (second mean-sd))
  446.                       (cov (diagonal (^ sd 2)))
  447.                       (var1 (^ stdev1 2))
  448.                       (var (^ sd 2))
  449.                       (rvdiff (/ (- var var1) var))
  450.                       (tilt mgfdel)
  451.                       (2tilt2 (* 2 (^ tilt 2)))
  452.                       (negtilt (- tilt)))
  453.                  (dotimes (i (length g) cov)
  454.                    (dotimes (j i)
  455.                      (let* ((g (select g (list i j)))
  456.                             (rvdi (select rvdiff i))
  457.                             (rvdj (select rvdiff j))
  458.                             (sdi (select sd i))
  459.                             (sdj (select sd j))
  460.                             (ldt1 (lapdet g tilt))
  461.                             (ldt2 (lapdet g negtilt))
  462.                             (del2 (/ (+ ldt1 ldt2) 2tilt2))
  463.                             (d (- del2 (* 0.5 rvdi) (* 0.5 rvdj)))
  464.                             (c (+ (aref cov1 i j) (* d sdi sdj))))
  465.                        (setf (aref cov i j) c)
  466.                        (setf (aref cov j i) c)))))))
  467.       (let ((mean-sd (transpose (mapcar #'moms2 mean1 stdev1 g-objects))))
  468.         (if covar 
  469.             (append mean-sd (list (covar g-objects mean-sd)))
  470.             mean-sd)))))
  471.  
  472. (defmeth bayes-model-proto :fullmoments (&optional g &key covar
  473.                                                    ((:derivstep h) .1)
  474.                                                    (maxiter 2))
  475.   (let* ((moms1 (send self :1stmoments g))
  476.          (mean1 (first moms1))
  477.          (stdev1 (second moms1))
  478.          (l-object (send self :transformed-logpost))
  479.          (g-objects (send self :transformed-functions g 0 mean1))
  480.          (loglap0 (send l-object :log-laplace nil maxiter nil h)))
  481.     (labels ((loglap (g tilt)
  482.                (- (send l-object :tilt-laplace g tilt nil maxiter nil h)
  483.                   loglap0))
  484.              (moms2 (g mu)
  485.                (let ((mu1 (exp (loglap g 1.0)))
  486.                      (mu2 (exp (loglap g 2.0))))
  487.                  (* mu (list mu1 (sqrt (- mu2 (^ mu1 2)))))))
  488.              (covar (g mean-sd)
  489.                (let* ((mu (/ (first mean-sd) mean1))
  490.                       (sd (second mean-sd))
  491.                       (cov (diagonal (^ sd 2))))
  492.                  (dotimes (i (length g) cov)
  493.                    (dotimes (j i)
  494.                      (let* ((g (select g (list i j)))
  495.                             (muij (exp (loglap g 1.0)))
  496.                             (mui (select mu i))
  497.                             (muj (select mu j))
  498.                             (mu1i (select mean1 i))
  499.                             (mu1j (select mean1 j))
  500.                             (c (* (- muij (* mui muj)) mu1i mu1j)))
  501.                        (setf (aref cov i j) c)
  502.                        (setf (aref cov j i) c)))))))
  503.       (let ((mean-sd (transpose (mapcar #'moms2 g-objects mean1))))
  504.         (if covar 
  505.             (append mean-sd (list (covar g-objects mean-sd)))
  506.             mean-sd)))))
  507.  
  508. (defmeth bayes-model-proto :2ndmoments (&rest args)
  509.   (apply #'send self :mgfmoments args))
  510.  
  511. (defmeth bayes-model-proto :moments (&rest args)
  512. "Args: (&optional gfuns &key covar) 
  513. Computes second order approximations to posterior moments. GFUNS can be
  514. a parameter index, list of indices, a function of the parameters or a
  515. list of such functions. Returns a the list of second order approximate
  516. means and standard deviations if COVAR is NIL. If COVAR is T the
  517. covaraince is appended to the end of the result as well."
  518.   (apply #'send self :2ndmoments args))
  519.  
  520. ;; margins
  521.  
  522. (defproto laplace-margin-proto '(logpost g x val i j a grad gval lu h))
  523.  
  524. (defmeth laplace-margin-proto :isnew (logpost g n k h)
  525.   (setf (slot-value 'logpost) logpost)
  526.   (setf (slot-value 'g) g)
  527.   (setf (slot-value 'x) (repeat 0 (+ n k)))
  528.   (setf (slot-value 'i) (iseq n))
  529.   (setf (slot-value 'j) (+ n (iseq k)))
  530.   (setf (slot-value 'a)
  531.         (make-array (list (+ n k) (+ n k)) :initial-element 0))
  532.   (setf (slot-value 'h) h)
  533.   (send self :adjust-internals t))
  534.  
  535. (defmeth laplace-margin-proto :adjust-internals (&optional initial)
  536.   (let* ((logpost (slot-value 'logpost))
  537.          (g (slot-value 'g))
  538.          (i (slot-value 'i))
  539.          (j (slot-value 'j))
  540.          (x (slot-value 'x))
  541.          (a (slot-value 'a))
  542.          (h (slot-value 'h))
  543.          (y (select x i))
  544.          (lambda (select x j))
  545.          (n (length y))
  546.          (vals (if initial 
  547.                    (list 0 (repeat 0 n) (- (identity-matrix n)))
  548.                    (send logpost :vals y h)))
  549.          (val (first vals))
  550.          (grad (second vals))
  551.          (hess (third vals))
  552.          (gvals (mapcar #'(lambda (x) (send x :vals y h)) g))
  553.          (gval (mapcar #'first gvals))
  554.          (ggrad (mapcar #'second gvals))
  555.          (ghess (mapcar #'third gvals))
  556.          (ggradmat (apply #' bind-columns ggrad)))
  557.     (setf (slot-value 'val) val)
  558.     (setf (slot-value 'grad) (apply #'+ grad (* lambda ggrad)))
  559.     (setf (slot-value 'gval) gval)
  560.     (setf (select a i i) (apply #'+ hess (* lambda ghess)))
  561.     (setf (select a i j) ggradmat)
  562.     (setf (select a j i) (transpose ggradmat))
  563.     (setf (slot-value 'lu) (lu-decomp a))))
  564.  
  565. ;; **** test for nonsingularity?
  566. (defmeth laplace-margin-proto :move-to (target)
  567.   (let* ((x (slot-value 'x))
  568.          (grad (slot-value 'grad))
  569.          (gval (slot-value 'gval))
  570.          (lu (slot-value 'lu))
  571.          (next-x (- x (lu-solve lu (combine grad (- gval target))))))
  572.     (setf (slot-value 'x) next-x)
  573.     (send self :adjust-internals)))
  574.  
  575. (defmeth laplace-margin-proto :log-density (&optional profile)
  576.   (let ((val (slot-value 'val)))
  577.     (if profile 
  578.         val
  579.         (let* ((lu (slot-value 'lu))
  580.                (nonsing (null (fourth lu))))
  581.           (if nonsing
  582.               (+ (* -0.5 (sum (log (abs (diagonal (first lu))))))
  583.                  val))))))
  584.  
  585. ;; ***** fix step choice
  586. ;; ***** Cut off at first nil?
  587. (defmeth bayes-model-proto :log-margin1 (g x &key 
  588.                                            ((:derivstep h) .05)
  589.                                            (spline t)
  590.                                            profile)
  591.   (let* ((moms1 (send self :1stmoments g))
  592.          (mean1 (select (first moms1) 0))
  593.          (stdev1 (select (second moms1) 0))
  594.          (n (send self :parameter-dimension))
  595.          (l-ob (send self :transformed-logpost))
  596.          (g-obs (send self :transformed-functions g mean1 stdev1))
  597.          (xs (/ (- x mean1) stdev1))
  598.          (xlow (coerce (sort-data (select xs (which (<= xs 0)))) 'list))
  599.          (xhigh (coerce (sort-data (select xs (which (> xs 0)))) 'list)))
  600.     (flet ((margin (x)
  601.              (let ((margin (send laplace-margin-proto :new l-ob g-obs n 1 h)))
  602.                (flet ((nextmargin (x)
  603.                         (send margin :move-to x)
  604.                         (send margin :log-density profile)))
  605.                  (mapcar #'nextmargin x)))))
  606.       (let* ((ylow (reverse (margin (reverse xlow))))
  607.              (yhigh (margin xhigh))
  608.              (x (append xlow xhigh))
  609.              (y (append ylow yhigh))
  610.              (i (which (mapcar #'numberp y)))
  611.              (xi (select x i))
  612.              (yi (select y i))
  613.              (xy (if spline (spline xi yi) (list xi yi))))
  614.         (list (+ mean1 (* stdev1 (first xy)))
  615.               (- (second xy) (log stdev1) (* 0.5 (log (* 2 pi)))))))))
  616.  
  617. (defmeth bayes-model-proto :margin1 (g x &key 
  618.                                        (derivstep .05)
  619.                                        (spline t)
  620.                                        profile)
  621. "Args: (g x &key (:derivstep .05) (spline t) profile)
  622. Computes Laplace approximation to marginal posterior density of G at
  623. points X. G can be an index or a function of the parameter vector. X
  624. is a sequence that should include the modal value of G. If SPLINE is
  625. true the log density is splined. If PROFILE is true, a profile of the
  626. posterior is returned."
  627.   (let* ((logmar (send self :log-margin1 g x 
  628.                        :derivstep derivstep 
  629.                        :spline spline
  630.                        :profile profile)))
  631.     (list (first logmar) (exp (second logmar)))))
  632.     
  633. ;;**** allow domain test function
  634. (defmeth bayes-model-proto :impsample (&optional g &key (n 100) (df 2))
  635.   (let* ((l-ob (send self :transformed-logpost))
  636.          (g-obs (send self :transformed-functions g))
  637.          (k (send self :parameter-dimension))
  638.          (v (chisq-rand n df))
  639.          (z (* (normal-rand (repeat k n)) (sqrt (/ df v))))
  640.          (c (- (log-gamma (/ (+ k df) 2)) 
  641.                (log-gamma (/ df 2)) 
  642.                (* (/ k 2) (log (/ df 2))))))
  643.     (flet ((w (z)
  644.               (let ((lp (send l-ob :value z))
  645.                     (lt (* -0.5 (+ k df) (log (+ 1 (/ (sum (* z z)) df))))))
  646.                 (if (realp lp) (exp (- lp lt c)) 0)))
  647.            (gvals (z) (mapcar #'(lambda (g) (send g :value z)) g-obs)))
  648.       (list (mapcar #'gvals z) (mapcar #'w z)))))
  649.  
  650. (defmeth bayes-model-proto :impmoments (&optional g &key (n 100) (df 2))
  651.   (let* ((impsample (send self :impsample g :n n :df df))
  652.          (means (/ (reduce #'+ (* (first impsample) (second impsample))) 
  653.                    (reduce #'+ (second impsample))))
  654.          (x (mapcar #'(lambda (z) (^ (- z means) 2)) (first impsample)))
  655.          (vars (/ (reduce #'+ (* x (second impsample))) 
  656.                   (reduce #'+ (second impsample)))))
  657.     (list means (sqrt vars))))
  658.