[SciPy-user] computing Bayesian credible intervals : help on constrained optimisation schemes?

dmitrey dmitrey.kroshko@scipy....
Mon May 5 09:33:31 CDT 2008

I have tried to solve your problem via openopt.NSLP solver nssolve & 
openopt.GLP solver galileo, both give maxResidual=0.885, that is far 
from desired zero, so your system (with required non-negative solution) 
seems to have no solution.

As for using fsolve or other smooth-funcs intended tools, it (AFAIK) is 
senseless wrt non-smooth funcs, like your numerical integration yields.

Regards, D.

Johann Cohen-Tanugi wrote:
> Hello,
> I am attaching my script. I followed Neil's suggestion but it fails to 
> converge, due seemingly to issues with the fact that a and b must be 
> positive, which I cannot enforce with fsolve, AFAIK.
> I am ready to dive into constrained schemes, especially in openOpt, 
> but I was hoping for some suggestions beforehand as to which path to 
> follow to solve this problem now.
> I remind that I am trying to find a and b so that :
> integral from a to b of p(x) = Q
> and p(a)=p(b)
> where Q is given (0.95 in my script) and p is a Poisson posterior pdf 
> for ON/OFF source experiments. a,b, and x are source rates, and as 
> such are positive.
> People will have recognized the computation of a Bayesian credible 
> interval here!!
> thanks a lot in advance,
> 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
> ------------------------------------------------------------------------
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user

More information about the SciPy-user mailing list