I recently had to deal with the same issue when computing conditional probabilities involving numbers on the order of $10^{-10000}$ (because normalizing the probability distribution would have required a great deal of unnecessary calculations). The heart of this difficulty, which comes up repeatedly in statistical computation, concerns taking logarithms of sums which themselves can overflow or underflow the computing platform's form of numeric representation (typically IEEE doubles). In the present case (to abstract away all irrelevant detail) we may write
$$L = \prod_i \left(\sum_j x_{ij}\right)$$
where the products $x_{i1} = \prod_{s=1}^{S_i}L_{is}(y\space|\space \rho_A)\phi$ and $x_{i2} = \prod_{s=1}^{S_i}L_{is}(y\space|\space \rho_B)(1-\phi)$ are themselves best computed using logarithms. Naturally we would like to compute $L$ as
$$L = \exp(\log(L));\quad \log(L) = \sum_i \log \left(\sum_j x_{ij}\right)$$
and in most applications, having $\log(L)$ is all that is really needed. There is the rub: how to compute logarithms of sums? Dropping the (now superfluous) subscript $i$, and assuming the $x_j$ will be computed using logarithms, we are asking how to calculate expressions like
$$\log\sum_j \exp(\log(x_j))$$
without over- or underflow.
(Note, too, that in many statistical applications it suffices to obtain a natural logarithm to an accuracy of only one or two decimal places, so high precision is not usually needed.)
The solution is to keep track of the magnitudes of the arguments. When $x_j \gg x_k$, then (to within a relative error of $x_k/x_j$) it is the case that $x_j + x_k = x_j$.
Often this magnitude-tracking can be done by estimating the order of magnitude of the sum (in some preliminary approximate calculation, for instance) and simply ignoring any term which is much smaller in magnitude. For more generality and flexibility, we can examine the relative orders of magnitude as we go along. This works best when summands have the most similar orders of magnitude, so that as little precision as possible is lost in each partial sum. Provided, then, that we arrange to sum the terms from smallest to largest, we may repeatedly apply this useful approximation and, in the end, lose no meaningful precision.
In many statistical applications (such as computing log likelihoods, which frequently involve Gamma and Beta functions) it is handy to have built-in functions to return the logarithms of factorials $n!$ and binomial coefficients $\binom{n}{k}$. Most statistical and numerical computing platforms have such functions--even Excel does (GAMMALN). Use these wherever possible.
Here is R code to illustrate how this could be implemented. So that it will serve as (executable) pseudocode, only simple Fortran-like expressions are used: they will port easily to almost any computing environment.
log.sum <- function(y, precision=log(10^17)) {
#
# Returns the logarithm of sum(exp(y)) without overflow or underflow.
# Achieves a relative error approximately e^(-precision).
#
# Some precision may be lost in computing `log(1+exp(*))`, depending on
# the computing platform.
#
if (length(y) == 0) return(1)
log.plus <- function(a, b) ifelse(abs(b-a) > precision,
max(b,a), max(b,a) + log(1 + exp(-abs(b-a))))
y <- sort(y)
x <- y[1]
for (z in y[-1]) x <- log.plus(x, z)
return (x)
}
For example, let's compute $\log\sum_{i=0}^{1000}\exp(i) = 1000.45867514538708\ldots$ in R with and without this method:
x <- seq(0, 1000, by=1); print(log(sum(exp(x))), digits=18)
[1] Inf
Summation simply overflows. However,
print(log.sum(x), digits=18)
[1] 1000.45867514538713
is accurate to 17 full decimal places, which is all one can expect of IEEE doubles.
