[SciPy-user] computing Bayesian credible intervals

[Johann Cohen-Tanugi]

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

Damian Eads wrote:
> 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
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>