[SciPy-user] computing Bayesian credible intervals
Neil Martinsen-Burrell
nmb@wartburg....
Fri May 2 13:02:27 CDT 2008
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
More information about the SciPy-user
mailing list