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