# Function given exact arguments returns hugely different value than it returns when given equivalent inexact arguments

I was trying to compute the probability that a coin is from a particular underlying distribution given that a particular set of tosses was observed. (I know this can be done in a different way, but I was just thinking of using Bayes in this context.)

I wrote this function:

PPP[heads_, total_, a_, b_] :=
Assumptions -> {0 < a, a < b, b < 1, 0 <= heads, heads <= total}]) /


When I evaluate

N[PPP[1000, 2000, 495/1000, 505/1000]]


I get

0.345408,


which seems reasonable. (The probability that your coin's distribution is between 0.495 and 0.505 given that you've witnessed 1000 heads in 2000 tosses is 0.345408.)

but, if I do

PPP[1000, 2000, 0.495, 0.505]


I get

3.936310242253867*10^468.


This makes no sense just on the face of it. The integral at the top is a subset of the integral at the bottom, and the integrand is strictly positive, so there's no way the numerator can be greater than the denominator.

Why is this happening? (I'm using Mathematica 9.0.0). Here's the notebook with those equations, if you don't want to cut and paste: ppp.nb.

• I have a feeling this is a numerical issue to do with the fact you're raising very small numbers to very big powers (leading to even smaller numbers). The first output that uses purely symbolic calculations is unaffected by it.
– Aky
Commented Jun 29, 2014 at 22:38
• Control your precision... try e.g. PPP[1000, 2000, 0.4951000, 0.5051000].
– ciao
Commented Jun 29, 2014 at 23:51

You could use NIntegrate on inexact numbers. It monitors precision and you won't get such catastrophic round-off error; and it will probably be faster as well. Something like this will set the WorkingPrecision to the minimum precision of the inputs, if the arguments are inexact numbers; otherwise it uses the OP's exact code:

Clear[PPP];
Precision[{heads, total, a, b}] < Infinity :=
WorkingPrecision -> Precision[{heads, total, a, b}]] /
NIntegrate[
WorkingPrecision -> Precision[{heads, total, a, b}]];
Assumptions -> {0 < a, a < b, b < 1, 0 <= heads, heads <= total}]) /


Example:

PPP[1000, 2000, 0.49516, 0.50516]
Precision[%]
(*
0.3454082823427782
15.699
*)


Warning: Each integral vanishes at MachinePrecision:

PPP[1000, 2000, 0.495, 0.505]


NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>
...
Power::infy: Infinite expression 1/0. encountered. >>
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

(* Indeterminate *)


Edit update: timing

To elaborate on the speed of NIntegrate:

AbsoluteTiming[PPP[1000, 2000, 495/1000, 505/1000];] (* Integrate *)
AbsoluteTiming[PPP[1000, 2000, .49516, .50516];]   (* NIntegrate, WorkingPrecision -> 16 *)
(*
{1.016703, Null}
{0.334477, Null}
*)


Interestingly, the assumptions with Integrate above make the integral almost five times faster. (On the face of it, m_goldberg seemed correct: for numeric arguments the assumptions are mathematically redundant.) From m_goldberg's code:

Clear[PPP];

PPP[1000, 2000, 0.4951000, 0.5051000]