[SciPy-user] computing Bayesian credible intervals

[Damian Eads]

Johann Cohen-Tanugi wrote:
> Hello,
> The question is a bit more general than that : How can I use SciPy to 
> compute the values a and b defined by :
> 1) integral from a to b of a pdf p(x) (or any function for that matter) 
> = a given value q
> 2) p(a)=p(b)
> 
> thanks in advance,
> Johann

Computing the area of the posterior over a credible interval involves 
integration. If your prior and likelihood are conjugate, it might be 
easier to use a conjugate distribution parameterized with the posterior 
hyperparameters and then compute CDF(b)-CDF(a). See 
http://en.wikipedia.org/wiki/Conjugate_prior for a list of conjugate 
priors with the hyperparameters worked out.

Now, I guess what your really asking is how to do the inverse of that 
(question #1), i.e. how do you find the end points of the interval if 
you know the area? Try the inverse CDF or inverse survival function. In 
Scipy, some distributions have an isf member function.
    b=post.isf(q)
    a=0.0

Now onto question #2, let's assume your posterior distribution is 
symmetric, you can try the inverse CDF or the inverse survival function. 
  For example, if q=0.7 (70%) and the posterior is symmetric, then
     L=[post.isf(0.85), post.isf(0.15)]
     a=min(L)
     b=max(L)
Note, post.pdf(a) should be equal to post.pdf(b) because post is symmetric.

Cheers,

Damian Eads