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

[Bruce Southey]

Hi,
I think that you need an alternative approach here because:

1) The Poisson is a discrete distribution not continuous (so can not use 
the solvers).
2) The Poisson is also a skewed distribution so finding points such that 
p(a)=p(b) is impossible unless the Poisson parameter is sufficiently large.
3) Depending on the parameters, it is impossible to get exact 
probabilities like 5% or 95% from discrete distributions.

These issues probably are reflected in the reported problems.

Depending on your parameter and hence how accurate do you want for your 
interval, Normal/Gaussian can provide a suitable approximation.

If you must stay with the Poisson, you need to solve it by brute force 
and allow for p(a) != p(b).

Regards
Bruce

dmitrey wrote:
> 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 at scipy.org
>>> http://projects.scipy.org/mailman/listinfo/scipy-user
>>>   
>>>       
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> SciPy-user mailing list
>> SciPy-user at scipy.org
>> http://projects.scipy.org/mailman/listinfo/scipy-user
>>     
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
>