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