[SciPy-user] computing Bayesian credible intervals

Johann Cohen-Tanugi cohen@slac.stanford....
Tue May 6 02:45:58 CDT 2008


hi Anne, thanks! I will need some time to digest your suggestions and 
try to implement them. I think that the condition p(a)=p(b) where p is 
the actual pdf and not an integral of it, is the sufficient condition 
for unimodal distributions so that the credible interval is uniquely 
defined as the interval encompassing the mode.
best,
Johann

Anne Archibald wrote:
> 2008/5/4 Johann Cohen-Tanugi <cohen@slac.stanford.edu>:
>
>   
>>  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
>>     
>
> This bites me every time I use quad. Still, it's valuable to check the
> reported error, though it's not always clear what to do with the
> results.
>
>   
>>  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".
>>     
>
> Multidimensional solving, as you are discovering, is a pain. If at all
> possible it's better to turn it into one-dimensional optimization. It
> can also be worth disposing of constraints - at least "open"
> constraints - using funny reparameterizations (log, for example, is
> great for ensuring that positive things stay positive)
>
> Do you really need p(a)=p(b)? I mean, is this the right supplementary
> condition to construct a credible interval? Would it be acceptable to
> choose instead p(x<a)=p(x>b)? This will probably be easier to work
> with, at least if you can get good numerical behaviour out of your
> PDF.
>
> Have you looked into writing down an analytical CDF? It looks
> moderately ghastly - all those s**n exp(k*s) - but possibly something
> that could be done with a bit of sweat. This would improve the
> performance of whatever solution method you use - the roundoff errors
> from quad() follow some irregular, discontinuous pattern. If you must
> use numerical integration, look into using a non-adaptive Gaussian
> quadrature, which won't have that particular problem. (There's one in
> scipy for more or less this reason.) Analytical mean, mode, or median
> could probably be made use of too. Or a different expansion for small
> s. A session with SAGE/MAPLE/Mathematica might yield something.
>
> If you do need p(a)=p(b), perhaps the way to turn the problem
> one-dimensional and reduce your headaches is to use a as the sole
> parameter and write b as a function of a using the fact that your
> distribution is single-peaked (so that p(x)-p(a)==0 has just two
> solutions, one of which you know about already; if you feel like being
> devious you could do root-finding on (p(x)-p(a))/(x-a) ).
>
> Good luck,
> Anne
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>   


More information about the SciPy-user mailing list