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