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