On 06/29/11 13:05, Invisible wrote:
> Suppose you flip a coin 1000 times. What is the probability of getting
> 402 heads?
> 
> The answer is of course that the head count follows a binomial
> distribution. That means that you need to compute the probability mass
> function. Unfortunately, in this instance the PMF in question is
> 
>   Pr(402 heads) = (1000 choose 402) / 2^1000
> 
> Now this is a problem. 2^1000 is roughly the square of the number of
> atoms in the observable universe, while 1000 choose 402 is defined as
> 
>   1000 choose 402 = 1000! / (402! * 598!)
> 
> In principle it looks plausible that you could compute the factorial of
> a big number. It might just be a bit slow, eh? Well... no. It turns out
> that 1000! is over 3 million bits in length, for a start.
> 
> Even if you could compute 1000 choose 402, it's again almost as huge as
> 2^1000. And taking the quotient of two vast numbers is unlikely to be
> very numerically stable.

Division by a power of 2 is numerically very stable for binary
arithmetics: It is just subtraction of the exponent. :-)

Moreover: If you feed the formula above into a symbolic calculator it is
able to

  a) compute the binomial more efficient than by division of three
     factorials.
  b) eliminate each factor of 2 in the numerator by one in the
     denominator so the division does not operates on this high numbers.


So the program 'maxima' prints after a very short time:


(%i1) binomial(1000,402)/(2^1000);
(%o1) 138375989102164508869118744263820047102478631240107164548096353\
437104485769496863239251132205776942412106967230343812796549218291211\
502370211070909495252274618853753553269750995170002584212749043301928\
705163749324086671038996613530725409460218145152701829040905136436576\
746556205530054031875 / 133938575898283415118553131132500226320175601\
463191700930468798546293881390617015311649797351961982265949334114694\
143353148393160711539255449807219683732185049182097185302887317763432\
563279639273474427276913080937294774265842484594489569299325963286432\
1399559710817770957553728956578048354650708508672

(I put spaces around the / to make it clearer)

To get a floating point approximation of the quotient:

(%i2) float(binomial(1000,402)/(2^1000));
(%o2)                       1.0331302104275842E-10


So where is the problem in calculating such expressions?


			Lars R.

Post a reply to this message