[SciPy-user] computing Bayesian credible intervals

[Neil Martinsen-Burrell]

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