How do I implement a numerically stable weighted logaddexp?

Question

What is the most numerically stable way of calculating:

log[(wx * exp(x) + wy * exp_y)/(wx + wy)]

where the weights wx, wy > 0?

Without the weights, this function is logaddexp and could be implemented in Python with NumPy as:

tmp = x - y
return np.where(tmp > 0,
                x + np.log1p(np.exp(-tmp)),
                y + np.log1p(np.exp(tmp)))

How should I generalize this to the weighted version?

Answer 1

You could use the original logaddexp function for thus purpose, if you rewrite the weighted expression as,

This is equivalent to,

logaddexp( x + log(w_x), y + log(w_y) ) - log(w_x + w_y)

which should be as numerically stable as the original logaddexp implementation.

Note: I'm referring to the numpy.logaddexp function that takes in x and y, not x and exp_y, as you mention in the question.