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?