# 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 Jun 29 '14 at 22:38
Control your precision... try e.g. PPP[1000, 2000, 0.4951000, 0.5051000]. – ciao Jun 29 '14 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];

AbsoluteTiming[PPP[1000, 2000, 495/1000, 505/1000];]
(*
{5.040427, Null}
*)

-

By not telling Mathematica otherwise, you have only a few significant digts in the second example, so the result loses accuracy quickly. You can tell Mathematica that a number has more precision than the digits you present, such as

PPP[1000, 2000, 0.4951000, 0.5051000]


See Precision, Accuracy, Numerical Precision, and Controlling the precision and accuracy of numerical results for more information.

-