Group Abstract Group Abstract

Message Boards Message Boards

NIntergrate for values very close to zero?

Posted 9 years ago

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?

POSTED BY: Niall Hughes
3 Replies
Posted 9 years ago

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.

POSTED BY: Niall Hughes

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]

enter image description here

Note that

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

gives 1.

Cheers,

Marco

POSTED BY: Marco Thiel

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}]
POSTED BY: Gianluca Gorni
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard