[SciPy-user] computing Bayesian credible intervals
Johann Cohen-Tanugi
cohen@slac.stanford....
Sun May 4 14:41:38 CDT 2008
hi Neil,
thanks a lot. I needed to modify :
integral_error = quad(density, a , b) - q
into
integral_error = quad(density, a , b)[0] - q
as quad returns a tuple with the result and the estimated error
More interesting, as a test I used :
q=0.96
def density(s):
if s<0.5:
return 4*s
else :
return 4-4*s
For which the solution is a=0.1, and b=0.9
If I keep
results=optimize.fsolve(solve_me, [0, 1])
fsolve fails because it wanders in directions a<0 and b>1.
With
results=optimize.fsolve(solve_me, [0.5, 0.5])
fsolve converges without trouble.... Is there a simple way to get it to
work with constraints like a>0 and b<1? I am afraid that the answer is
no, if I understood well the context of another recent thread entitled
"constrained optimization".
Thanks a lot,
Johann
Neil Martinsen-Burrell wrote:
> Johann Cohen-Tanugi <cohen <at> slac.stanford.edu> writes:
>
>
>> 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.
>>
>
> scipy.optimize.brentq is Brent's method for finding a root of a given scalar
> equation. Since you are looking for two values, a and b, with two conditions,
> then Brent's method is not appropriate (barring some symmetry-based reduction to
> one variable). I like to use scipy.optimize.fsolve to find roots of
> multivariable equations, such as
>
> def solve_me(x): # x is an array of the values you are solving for
> a,b = x
> integral_error = quad(density, a , b) - q
> prob_difference = density(b) - density(a)
> return np.array([integral_error, prob_difference])
>
> fsolve(solve_me, [0.0, 1.0]) # initial guess is a = 0, b = 1
>
>
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-user
mailing list