Numerically stable calculation of multinomial probabilities

103 Views Asked by At

I'm looking for a numerically stable method to compute expressions of the form $$\frac{(a+b+c+d)!}{a!b!c!d!}\left(\frac{1}{4}\right)^{a+b+c+d}$$

So far I've been using a compensated sum algorithm to compute

$$\log \left(\frac{1}{1}\frac{1}{4}\right) + \ldots + \log \left(\frac{a+1}{1}\frac{1}{4}\right) + \ldots + \log \left(\frac{a+b+c+d}{d}\frac{1}{4}\right)$$

and then take the exponential. It works, but there's probably a faster method that doesnt require computing the log and the exponential.

1

There are 1 best solutions below

0
On

From a coding point of view, I should first generate a table of $x_i=\log(i)$ for all values of $i$ from $1$ to $(a+b+c+d)$. So, $$A=\frac{(a+b+c+d)!}{a!\,b!\,c!\,d!}\left(\frac{1}{4}\right)^{a+b+c+d}$$ Taking logarithms $$\log(n!)=\sum_{i=1}^n \log(i)=\sum_{i=1}^n x_i$$ Doing so, the logarithms are only computed once.

We could do better if we set the numbers such that $a<b<c<d$. In such a case, compute first $$\log(a!)=\sum_{i=1}^a x_i$$ $$\log(b!)=\log(a!)+\sum_{i=a+1}^b x_i$$ $$\log(c!)=\log(b!)+\sum_{i=b+1}^c x_i$$ $$\log(d!)=\log(c!)+\sum_{i=c+1}^d x_i$$ $$\log\big((a+b+c+d)!\big)=\log(d!)+\sum_{i=d+1}^{a+b+c+d} x_i$$