GENEALOGY-DNA-L Archives

Archiver > GENEALOGY-DNA > 2002-07 > 1026177000


From: "John F. Chandler" <>
Subject: Re: [DNA] Poisson vs. others
Date: Mon, 8 Jul 2002 21:10 EDT
In-Reply-To: gbonner@smhrimsx.mhri.med.umich.edu message <2449F3F867FBD211936600C04F68569EA3CB05@smhrimsx.mhri.med.umich.edu> of Sat, 6 Jul 2002 16:56:12 -0600


Gregg wrote:
> I don't know any other way. In my case, I can calculate 432! on other
> computers which are readily available to me. Nevertheless, I would be
> interested to learn how to determine exact binomial without having to
> determine the (n) factorial.

The point is that the binomial coefficient "N take M" is N!/M!/(N-M)!,
and that expression is all you need. You don't need the individual
factors. Therefore, if you find yourself in the regime where M is close
to N, you note that N!/M! = Nx(N-1)x(N-2)x...x(M+1). What's more, you
can do the whole thing in purely integer arithmetic by alternating:
N x (N-1) / 2 x (N-2) / 3 x ... x (M+1) / (N-M). (This is because of
a theorem that shows successively that Nx(N-1) is always divisible
exactly by 2, and Nx(N-1)x(N-2) is always divisible exactly by 2x3,
and so on.)

> I would like to see my exact example in sensible implementation without
> using the factorial. For a case of 36 transmission events, and 12 markers
> per event (viz., 432 opportunities to mutate), how would I calculate the
> probabilities of finding the range of number of mutations numbering zero to
> 10 without using as a determined intermediate number something that has an
> exponent with a magnitude of > +100?

If the rate of mutation is R, then the probability of 0 mutations is
(1-R)^432. Take that result and follow the chain of multiplication and
division for the successive binomial coefficients, and you get all the
probabilities you want as intermediate results:

P(0) x 432 x R/(1-R) = P(1)
P(1) x 431 / 2 x R/(1-R) = P(2)
P(2) x 430 / 3 x R/(1-R) = P(3)
P(3) x 429 / 4 x R/(1-R) = P(4)
P(4) x 428 / 5 x R/(1-R) = P(5)
P(5) x 427 / 6 x R/(1-R) = P(6)

and so on.

John Chandler

This thread: