On 29/06/2011 12:05 PM, 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!)

Just for giggles, I tried it.

It turns out that the GMP can easily handle factorials of this size. I 
was even able to convert this to a double-precision floating point 
number and divide by 2^1000 to get the correct answer.

However, if the parameters become sufficiently large, it eventually 
stops working. Apparently the "choose" function yields its largest 
result when k = n/2. 1028 choose 514 works OK. However, 1030 choose 515 
is so huge that it registers as "infinity". The rest of the calculation 
then fails to work as a result.

Performing the division by 2^n in infinite-precision arithmetic and then 
converting to double-precision allowed me to get the correct probability 
value for 1030 flips with 515 heads. However, if we now look at the 
lowest probability, around 1060 flips the probability reported for zero 
heads began to get progressively less accurate. (Denormal numbers?) At 
1070 flips 0 heads the probability actually came out as "zero". (It 
should in fact be 7.7 * 10^326.) Still, the fraction representation was 
quite correct; it's merely that double-precision cannot handle such tiny 
numbers.

In short, it seems there are limits to what can be done without some 
very clever numerical programming. Still, I got a lot further than I 
ever expected to get...

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message