[SciPy-user] computing Bayesian credible intervals

[Johann Cohen-Tanugi]

hi Neil, thanks for your answer and sorry I was not clear enough. Of 
course I require the 2 conditions. 1) defines *a* credible interval if p 
is a posterior pdf; and 2) sets a constraint that for common situation 
yield *the* standard Bayesian credible interval. I will have a look at 
brentq, I do not know what it refers to.
best,
Johann


Neil Martinsen-Burrell wrote:
> Johann Cohen-Tanugi <cohen <at> slac.stanford.edu> writes:
>
>   
>> 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)
>>     
>
> Problem 1 is under-specified.  You have one equation for two unknowns and
> generally will not be able to solve uniquely for both of the endpoints.  It may
> be common in your field to require some sort of symmetry between a and b, e.g. a
> = mu + E, b = mu - E where mu is some fixed quantity and now you solve for E
> rather than a and b separately.  I (perhaps naively) would do this in scipy
> using scipy.optimize.brentq such as
>
> form numpy import exp
> from scipy.optimize import brentq
>
> mu = 0
>
> def density(x):
>     return exp(-x)
>
> def integral_function(E, mu):
>     return scipy.integrate(density, mu-E, mu+E) - q
>
> margin = brentq(integral_function, 0, 1000)
> a,b = mu - margin, mu+margin
>
> For problem 2, I'm not sure what the function p represents, but a suitable
> adaptation of the above should work.
>
> -Neil
>
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>