# Numbers too large for R. How to approximate probability mass function?

Social network data is frequently found in a two-mode form: people vs. events they attend, people vs. classes they attend, countries vs. treaties they sign, etc. A strategy for analyzing this data is to project the rectangular, binary matrix $X$ into a one mode matrix $P = (X * X')$. Then each cell of the new matrix $A_{ij}$ would have the number of times person $i$ co-attended an event with person $j$. But is the number of time they co-attend events more than what would be expected by chance?

I found an interesting paper on the topic and it addresses the issue directly. The author proposes this PMF, where the probability of person $i$ and person $j$ attending exactly $C$ events:

$$\Pr(P_{ij}=C)=\frac{{E \choose C}{E-C \choose P_{ii}-C}{E-P_{ii}\choose P_{jj}-C}}{{E \choose P_{ii}}{E \choose P_{jj}} }$$ In a reasonably small network, there is no difficulty calculating this. But I have a network with thousands of nodes. The numbers in the numerator and denominator are enormous. So large that R just returns Inf and I get a meaningless result.

What I think I should do is find a way to approximate this PMF somehow. I was also considering just writing some code that approximates the distribution through simulation. Is there a better way of approximating this distribution? Is there some known distribution that is very close (read. close enough) to the the theoretical distribution presented?

• Could you use Stirling's approximation and work in log-space? – Nicholas Mancuso Sep 20 '14 at 19:33

Pretty much any decent stats package will provide a log-gamma or log-factorial function.

You mention R; it has:

• lgamma which is the log of the gamma function

• lfactorial which is the log of the factorial function

• lchoose which is the log of the binomial coefficient.

using any of these, you can work out the log of the desired probability. If it's not going to cause underflow to do so, you can always exponentiate it at the end.

See ?gamma

An alternative if you don't have such a function is to keep all the terms from each binomial coefficient in bins (that is, if there's an "11" from expanding a binomial coefficient on a numerator, add '1' to an "11" bin, and if there's an "11" on a denominator, subtract '1'. After you've made your way through all of the coefficients, you can multiply and divide in such an order as to keep the result not too far from 1 (at least until you run out of numerator terms). An advantage of this approach is that you can keep results as exact fractions if you wish. (You can make it more sophisticated by then cancelling common factors before starting the multiplication and division ... but that's not likely to be worthwhile if you just want a numeric answer.)

A third alternative is to generate approximate answers via Stirling's approximation, but this shouldn't be necessary (if I was working it out in my head, I'd do it this way).

• +1. For an example using logs, see the code in the question math.stackexchange.com/questions/465318/… which calculates $\displaystyle\sum_{i=0}^n (-2)^i {n \choose i}\frac{(2n-i)!}{(2n)!}$ for $n= 10^6$ – Henry Sep 20 '14 at 20:18