Multiply a very large integer by a very small probability

Question

I am calculating the probability of rejecting a one-tailed hypothesis test for inequality of two proportions $p_1$ and $p_2$ (i.e., power) using $H_o: p_2 \leq p_1$ vs. the alternative $H_a: p_2 > p_1$ using the exact formulation

$$ Prob=\sum_{i=0}^m \sum_{j=0}^n \binom{m}{i} \binom{n}{j} p_1^i \cdot (1 - p_1)^{m - i} \cdot p_2^j \cdot (1 - p_2)^{n - j} \cdot I(T>z_{1-\alpha}), $$

where $T = \frac{ p_j - p_i}{\sqrt{\bar{p} (1 - \bar{p})}\sqrt{\frac{1}{m}+\frac{1}{n}}}$, $p_i = i / m$, $p_j = j / n$, $\bar{p} = (i + j) / (m + n)$, and typical values of $m$ and $n$ are e.g. $80$ for my use case. $I()$ is a boolean indicator function. The binomial coefficient values $\binom{m}{i}$ take on very large integer values that are $\gg 10^{308}$ - the limit of double precision memory storage for most compilers. Analogously, the probabilities $p_1^i \cdot (1 - p_1)^{m - i} \cdot p_2^j \cdot (1 - p_2)^{n - j}$ can also be e.g. $\ll 10^{-308}$.

Are there any tricks for representing very small probabilities using large integers, or some reliable algorithm that would enable me to perform such multiplications. The limiting issue is that there doesn't seem to be one or more fast methods for memory representation of a real value with $\gg 308$ digit precision.

Answer 1

Do you need $308$ digits of precision or just $10^{308}$ of range? If you need the precision, there is software that can do that, like Mathematica. If you just need the range, you can work with logarithms until the large and small numbers multiply to get in range. Stirling's approximation makes logs of the factorials easy. For the powers of $p$ and $(1-p)$ you can take the logs of $p$ and $1-p$. If one is small and the other is close to $1$ it may be better to use the Taylor series of $\log (1-x)$ for the one that is close to $1$.

Answer 2

This is an answer to the question before its edit.

Given $$\text{Prob}=\sum_{i=1}^n \sum_{j=1}^m \binom{n}{i}\binom{m}{j} p^{(i+j)}(1-p)^{(n+m-i-j)}$$ I should start rewiriting it as $$\text{Prob}=(1-p)^{(n+m)}\, \Gamma (m+1)\, \Gamma (n+1)\, \large \sum$$ with, using $\color{red}{q=\frac{p}{1-p}}$, $$\large \sum=\sum_{i=1}^n \sum_{j=1}^m \frac {q^{i+j}}{\Gamma (i+1) \Gamma (j+1) \Gamma (n-i+1) \Gamma(m-j+1) }$$ The inner sum write $$\frac{q^i \left((q+1)^m-1\right)}{\Gamma (i+1) \Gamma (m+1) \Gamma (n-i+1)}$$ making at the end $$\large \sum=\frac{\left((q+1)^m-1\right) \left((q+1)^n-1\right)}{\Gamma (m+1) \Gamma(n+1)}$$ Recombining and defining $\color{red}{t=\frac{1}{1-p}}$ $$\large\color{blue}{\text{Prob}=(1-t^{-m}) (1-t^{-n} )}$$ Computing $n\log(t)$ and $m\log(t)$ and exponentiating again is not supposed to make major problems.

Now, if $p$ is small, you could use series expansion and obtain $$\text{Prob}=mn\,p^2\Big(1+\sum_{k=1}^r a_k \, p^k \Big)$$ with $$a_1=-\frac 12 (m+n-2)$$ $$a_2=\frac 1{12}\left(2 m^2+3 m n-9 m+2 n^2-9 n+11\right)$$ $$a_3=-\frac 1{24} (m+n-4)\left(m^2+m n-4 m+n^2-4 n+5\right)$$ and so on.

Answer 3

You can calculate this binomial probabilities indirectly using a relationship with the beta distribution, discussed at our sister site here: https://stats.stackexchange.com/questions/4659/relationship-between-binomial-and-beta-distributions That way you avoid the problem with multiplication!

Answer 4

You can directly compute the answer only using bounded numbers. Let $f_p(m,i) := \binom{m}{i} p^i (1-p)^{m-i}$. Because it represents a probability, we have $0 \leq f_p(m,i) \leq 1$ and we have a recurrence:

$$ \begin{align} f_p(0,0) &= 1 \\ f_p(m,i) &= p f_p(m-1,i-1) + (1-p) f_p(m-1,i) & \text{if } 0 \leq i \leq m \\ f_p(m,i) &= 0 & \text{otherwise} \end{align} $$

therefore a dynamic programming algorithm computes the all values $f_p(m,i)$ given $p$ in $O(m^2)$ time.

To compute $\mathrm{Prob}$ in the question, compute the arrays $A_1, A_2$ such that $A_1[i] = f_{p_1}(m,i)$ for $0 \leq i \leq m$, and $A_2[j] = f_{p_2}(n,j)$ for $0 \leq j \leq n$, then

$$ \mathrm{Prob} = \sum_{i=0}^m \sum_{j=0}^n A_1[i] A_2[j] I(i,j) $$

for some function $I(\cdot,\cdot)$. It is total $O(n^2 + m^2 + nm)$ time, and all intermediate values are in $[0,1]$.