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