NIntergrate for values very close to zero?

Questions
, ,
1

I am trying to run the following. I am interested in b (the posterior distribution of x given p)

t = 0.8 - (x/4000);
p= Binomial[10000, 5000] t^5000 (1 - t)^(5000);
intp = NIntegrate[p, {x, 0, 1}];
b= p/intp;
Plot[{b}, {x,
   0, 1}, {WorkingPrecision -> 200, PlotLegends -> Automatic} ]

p is a very small number because t is very flat. But, p does still vary in x, and should have most mass around x=1.

when I integrate p , mathematica struggles with the values of p being very close to zero for the whole range of x. It cannot evaluate precisely enough…returning a value of zero.

Then b=p\intp gives an error message as mathematica is trying to divide by zero.

Any ideas?

2

This is a case where you can find the exact value for the integral, if you are willing to wait a few minutes:

t = 8/10 - (x/4000);
p = Binomial[10000, 5000] t^5000 (1 - t)^5000;
intp = Integrate[Expand[p], {x, 0, 1}];
intp // N
Plot[p/intp, {x, 0, 1}]
3

I also favour the analytical/exact solution, but this might work, too, and evaluates practically instantaneously.

t = 0.8 - (x/4000); 
p = Binomial[10000, 5000] t^5000 (1 - t)^(5000); 
intp = NIntegrate[p, {x, 0, 1}, WorkingPrecision -> 10]; 
b = p/intp; Plot[{b}, {x, 0, 1}, PlotLegends -> Automatic]

ScreenShot2016-08-26at13.51.01.png
ScreenShot2016-08-26at13.51.01.png760×496 11.7 KB

Note that

NIntegrate[b, {x, 0, 1}]

gives 1.

Cheers,

Marco

4

thanks a lot to both of you. The first one seems to work but takes quite a while. I’m had pretty good success so far with the second approach.

cheers.