but not really...
[CommonLispStat.git] / bayes.lsp
blob1a34fec598ec1c2707e5066167dcd4a32fd44b89
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:
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.
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
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*))
33 (in-package :lisp-stat-bayes)
35 ;;; Objects Representing Functions
37 ;; Generic C2 Functions
39 (defproto c2-function-proto '(f h num-derivs))
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))
46 (defmeth c2-function-proto :f (&optional f)
47 (if f (setf (slot-value 'f) f))
48 (slot-value 'f))
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))
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)))))
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))))
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))))
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)
77 (list (first v) (second v) (send self :hessian x h)))
78 (list v (send self :gradient x h) (send self :hessian x h)))))))
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))))
85 ;; Scaled C2 Functions
87 (defproto scaled-c2-function-proto '(scaling) () c2-function-proto)
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))))
112 ;; Tilted C2 Functions
113 ;; **** allow nil values?
114 (defproto tilt-function-proto '(tilt exptilt) () c2-function-proto)
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))
121 (defmeth tilt-function-proto :tilt (&optional tilt)
122 (if tilt (setf (slot-value 'tilt) tilt))
123 (slot-value 'tilt))
125 (defmeth tilt-function-proto :exptilt (&optional (new nil set))
126 (if set (setf (slot-value 'exptilt) new))
127 (slot-value 'exptilt))
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))))))
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))))))
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))))))
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))))))
179 ;; scaled log posterior prototype
181 (defproto scaled-logpost-proto
182 '(tilt-object init-pars) () scaled-c2-function-proto)
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))))
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)))))
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)))))
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)))
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)))
238 ;;;;
239 ;;;; Bayes Model Prototype
240 ;;;;
242 (defproto *bayes-model-proto* '(bayes-internals))
244 ;; initialization methods and constructor function
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)))
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))
270 ;; display method
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 "~%")))
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)))))
291 ;; implementation-dependent access methods
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)))
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)))
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)))
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))))
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)))
324 (defmeth *bayes-model-proto* :parameter-dimension ()
325 (length (select (slot-value 'bayes-internals) 1)))
327 (defmeth *bayes-model-proto* :derivstep ()
328 (select (slot-value 'bayes-internals) 3))
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)))
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)
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))))
352 (defmeth *bayes-model-proto* :mode ()
353 (if (send self :needs-maximizing) (send self :maximize))
354 (first (send self :mode-values)))
356 (defmeth *bayes-model-proto* :new-mode-guess (new)
357 (send self :needs-maximizing t)
358 (send self :mode-values new))
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)))
370 ;;**** need transformed domain here
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)))
387 ;; computing methods
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)))
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.~%"))))))))
428 ;; moments
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))))))
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)))))
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)))))
532 (defmeth *bayes-model-proto* :2ndmoments (&rest args)
533 (apply #'send self :mgfmoments args))
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))
544 ;; margins
546 (defproto laplace-margin-proto '(logpost g x val i j a grad gval lu h))
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))
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))))
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)))
601 (defmeth laplace-margin-proto :log-density (&optional profile)
602 (let ((val (slot-value 'val)))
603 (if profile
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))))))
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)))))))))
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)))))
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)))))
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))))