bayes.lsp

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