bayes.lsp

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