I want to use numpy's logsumexp()
in python 2.7.
The formula I need to solve looks like this:
log ( 1 + e^a1 + e^a2 + e^a3 + ... e^an - e^ax )
The last term which is a negative number just has to be appended on. Excluding this last term, I would do the following:
myarray = numpy.array([0, a1, a2, a3, ..., an])
That way, with the first element being 0
, then e^0 = 1
and so I have my first term, which is 1
. Then I would just use
result = numpy.logsumexp(myarray)
and I would get the correct result.
But now I have to append a -e^ax
, and because it's negative, I can't simply append ax
to the end of myarray
. I also can't append -ax
because that's just wrong, it would mean that I'm adding 1/e^ax
, instead of -e^ax
.
Is there any direct way to append this so that I can still use logsumexp()
? The only reason I'm insisting on using logsumexp()
rather than separately using numpy.exp()
and numpy.sum()
and numpy.log()
is because I have the impression that logsumexp
also incorporates stability within it in order to prevent underflows (correct me if I'm wrong). However if there's no other way around then I guess I have no choice.
According to scipy.misc.logsumexp
documentation:
scipy.misc.logsumexp(a, axis=None, b=None) Parameters: b: array-like, optional Scaling factor for exp(a). Must be of the same shape as a or broadcastable to a. New in version 0.12.0.
So, you could add list of factors like this:
In [2]: a = [0, 1, 3, 2]
In [3]: logsumexp(a, b=[1] * (len(a) - 1) + [-1])
Out[3]: 2.7981810916785101