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