# [SciPy-user] computing Bayesian credible intervals

Johann Cohen-Tanugi cohen@slac.stanford....
Tue May 6 01:24:52 CDT 2008

```Damian, Bruce, Dmitrey,
thank you very much for your answers. I feel I need to be more precise
as to what I am doing. Look at eqs 5.13 and 5.14 of the following article
http://www.astro.cornell.edu/staff/loredo/bayes/promise.pdf
These 2 eqs define the posterior probability I am trying to work with.
Bruce, it is a continuous function of the count rate s, though of course
the underlying Poisson process is discrete. Damian, I will definitely
look into your suggestions to understand better what is available in
that respect, but I am not expecting this pdf to be implemented for cdf
or isf computing. Dmitrey could you send me your calls to the openopt
optimizers, that would be helpful for me to get a headstart on using
your very nice framework and play with other possibilities then fsolve.

As eqs 5.13 and 5.14 show, the resulting pdf is by construction
normalized to 1. In my code, I tried integrate.quad and did not get a
result close to 1. So I think that I have a numerical issue in my python
formula, which would presumably also explain why it fails to converge
(the pdf seems to go down to 0 at low s value much too quickly).

thanks again,
Johann

> 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)
>>
>> 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,
>