Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

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.

share|improve this question
3  
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 at 22:38
4  
Control your precision... try e.g. PPP[1000, 2000, 0.495`1000, 0.505`1000]. –  rasher Jun 29 at 23:51
add comment

3 Answers 3

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:

Clear[PPP];
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}]] /
  NIntegrate[
     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}]);

Example:

PPP[1000, 2000, 0.495`16, 0.505`16]
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, .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:

Clear[PPP];
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
add comment

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
add comment

This answer is an elaboration of rasher's.

First, let's get rid of the assumptions; they have no bearing on the basic problem.

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}])

Now let's try inexact arguments of 100 digits of precision.

ppp100 = PPP[1000, 2000, .495`100, .505`100];
Row[{"value: ", Round[#, .00000001], "  acuracy: ", Accuracy[#]}] & @ ppp100

value: 1.915788213096238*10^384 accuracy: -288.211

Even with 100 digits we have lost all accuracy in the result.

So now we try 1000 digits.

ppp1000 = PPP[1000, 2000, .495`1000, .505`1000];
Row[{"value: ", Round[#, .00000001], "  acuracy: ", Accuracy[#]}] & @ ppp1000
value: 0.345408  accuracy: 515.379

This last suggests that we could get away with 500 digits.

ppp500 = PPP[1000, 2000, .495`500, .505`500]; 
Row[{"value: ", Round[#, .00000001], "  accuracy: ", Accuracy[#]}] & @ ppp500
value: 0.345408  accuracy: 15.379

Yes, all is well at 500 digits if we can be happy with a 6-place (or somewhat better) result.

Now let's look at timing

AbsoluteTiming[PPP[1000, 2000, 495/1000, 505/1000];]
{5.720626, Null}
AbsoluteTiming[PPP[1000, 2000, .495`500, .505`500];]
{14.292283, Null}

Moral: In this case, giving exact arguments and using symbolic evaluation not only produces more accurate results, it runs faster.

share|improve this answer
    
Keeping the assumptions made Integrate evaluate faster. See my update. Do you see the same improvement? –  Michael E2 Jun 30 at 5:20
add comment

Your Answer

 
discard

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.