Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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_] :=
  (Integrate[p^heads (1 - p)^(total - heads), {p, a, b},
     Assumptions -> {0 < a, a < b, b < 1, 0 <= heads, heads <= total}]) /
    (Integrate[p^heads (1 - p)^(total - heads), {p, 0, 1}, 
       Assumptions -> {0 <= heads, heads <= total}])

When I evaluate

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

I get


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


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.

share|improve this question
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
up vote 4 down vote accepted

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:

PPP[heads_?NumericQ, total_?NumericQ, a_?NumericQ, b_?NumericQ] /; 
   Precision[{heads, total, a, b}] < Infinity :=
 NIntegrate[p^heads (1 - p)^(total - heads), {p, a, b}, 
     WorkingPrecision -> Precision[{heads, total, a, b}]] /
     p^heads (1 - p)^(total - heads), {p, 0, 1}, 
     WorkingPrecision -> Precision[{heads, total, a, b}]];
PPP[heads_, total_, a_, b_] :=
 (Integrate[p^heads (1 - p)^(total - heads), {p, a, b}, 
     Assumptions -> {0 < a, a < b, b < 1, 0 <= heads, heads <= total}]) /
  (Integrate[p^heads (1 - p)^(total - heads), {p, 0, 1}, 
     Assumptions -> {0 <= heads, heads <= total}]);


PPP[1000, 2000, 0.495`16, 0.505`16]

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, .495`16, .505`16];]   (* 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:

PPP[heads_, total_, a_, b_] :=
  (Integrate[p^heads (1 - p)^(total - heads), {p, a, b}])/
   (Integrate[p^heads (1 - p)^(total - heads), {p, 0, 1}])

AbsoluteTiming[PPP[1000, 2000, 495/1000, 505/1000];]
  {5.040427, Null}
share|improve this answer

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.495`1000, 0.505`1000]

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

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.