python computing likelihood causing exp overflow

Question

I am using numpy to compute the likelihood of a variable $Z$ using numpy. $Z$ is a Bernoulli random variable which has two outcomes $[0,1]$. I compute the log likelihood of observing $Z$ given the parameter is $x=[ -3146,-1821]$. These numbers are too large, np.exp(x) returns [0,0], and np.exp(x)/np.exp(x).sum() returns [nan, nan]. I know the reason is that these two numbers are both smaller than np.log(sys.float_info.min)(the value is $-708.3964185322641$ in my PC). Is there any method to accurately calculate the likelihood from log likelihood if the numbers are too small?

I tried to add a large number on $x$: $x = x+1821$, but the first number is still too small.

Update

The context is a state observation problem. There are some (unknown states) $\{z_1,z_2,\dots,z_m\} \in Z$ whose true values are unknown, we only know that it can take two possible values $[0,1]$. $Z$ can only be observed by some sensors $\{s_1,s_2,\dots,s_n\} \in S$, each sensor can provide values to multiple unknown states in $Z$. The observed sensor value from $s_i$ about $z_j$ is $x_{ij}$. $e_i$ is used to model the error of a sensor $s_i$, larger $e_i$ implies sensor $s_i$ is less accurate and $e_i$ is unknown as well.

I model $Z$ as the latent variable, and $E$ as the model parameters where $E = \{e_i,e_2,\dots,e_n\}$. I am trying to use EM algorithm to work out $Z$ and $E$ given the observations $X$. The E-step is given below:

$p(z_j|X_j,E_j) \propto p(z_j)\prod\limits_i p(x_{ij}|z_j,e_i)$, where $X_j$ is the set of sensors values for $z_j$, $E_j$ is the set of errors of sensors that observe $z_j$, $p(z_j)$ is the prior and $p(x_{ij}|z_j,e_i)$ is the likelihood. Then we have:

$p(z_j = k|X_j,E_j) = \frac{p(z_j = k|X_j,E_j)}{\sum_{k'=0}^1 p(z_j = k'|X_j,E_j)}$

M-step is not related to the question here, not going to write it here.

Back to the question, I am computing the value of $\log p(z_j = k|X_j,E_j) = \log p(z_j=k) - \log\sum\limits_i p(z_j = k|X_j,E_j)$. Since there are a lot of sensors, the log likelihood is very small: $[-3146,-1821]$. Using np.exp([-3146,-1821]) to convert it back to likelihood would raise some warnings because the numbers are too small to be handled by np.exp.

Update 2

I understand that $\exp(-3146)$ is tiny and could be ignored if compared with $\exp(-1821)$. It is attempted to approx the result as $[0,1]$. But I think in terms of probability, $\exp(-3146)$ is indeed very very small, it still has a very tiny probability that $z_j=0$. If the approximated result $\exp([-3146,-1821])\approx[0,1]$ is used, it would make the probability $p(z_j = 0|X_j,E_j) = 0$, i.e., $z_j=0$ is impossible.

Answer 1

so if I understand what you are doing, you want to calculate $$\log \frac{\exp(x_i)}{\sum exp(x_j)}$$ where $x$ is large negative number. So what you need to do is write as $$\log \frac{\exp(x_i)}{\exp(x_m) (1 + \sum_{j\ne m} \exp(x_j - x_m))}=x_j - x_m - \log (1 + \sum_{j\ne m} \exp(x_j - x_m))\approx x_j-x_m$$

Where $x_m$ is the largest of your negative numbers (ie -1821 in your example) Hopefully this is sufficient for your purposes otherwise you would have to do Taylor series of $\log(1 + x)=x ... $

Answer 2

The cause of the problem is that numpy and python's math module cannot calculate $\exp(x)$ if $x$ is greater than sys.float_info.max of smaller than sys.float_info.min. I find some packages, e.g., mpmath and sympy, can compute precise exponent, but the efficiency is significantly worse than numpy, especially in the case of vectorized operations.

If programming in Python, the softmax function in scipy.special module can solve this problem. For example, softmax([-2000,-2005]) returns array([0.99330715, 0.00669285]).