[SciPy-user] computing Bayesian credible intervals

[Johann Cohen-Tanugi]

mamma mia, shame on me... :(
Of course in my code, it should have read :
 value+=coeff*(s*Ton)**(i-1)/sp.factorial(i-1)
and not
 value=coeff*(s*Ton)**(i-1)/sp.factorial(i-1)
Convergence is much better now with the simple fsolve!

so much for public humiliation. Thanks all for bearing with me :)
I am still interested in testing constrained solver, Dmitrey.

Johann

Robert Kern wrote:
> On Tue, May 6, 2008 at 2:18 AM, Anne Archibald
> <peridot.faceted at gmail.com> wrote:
>   
>>  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.
>>     
>
> It's one of the defining characteristics of the kind of credible
> interval Johann is looking for. "Bayesian credible interval" is a
> somewhat broad designation; it applies to pretty much any interval on
> a Baysian posterior distribution as long as the interval is selected
> according to some rule that statisticians agree has some kind of
> meaning.
>
> In practice, one of the most common such rules is to find the "Highest
> Posterior Density" (HPD) interval, where p(a)=p(b) and P(b)-P(a)=0.95
> or some such chosen credibility level. Imagine the PDF being flooded
> with water up to its peak (let's assume unimodality for now). We
> gradually lower the level of the water such that for both points a and
> b where the water touches the PDF, p(a)=p(b). We lower the water until
> the integral under the "dry" peak is equal to 0.95. Then [a,b] is the
> HPD credible interval for that PDF at the 0.95 credibility level.
>
>