OverFlowError during exponentiation

Question

I am having trouble calculating this. The code works for all values of N up to and including N = 57, but throws an overflow error (34, 'Result too large') for N greater than or equal to 58. Is there a way around this? Thanks.

import numpy as np
import scipy.integrate
import scipy.optimize
import warnings
warnings.filterwarnings('ignore')

R = 6
N = 58
Nb = 4
S_norm = 0.3
Ao = 1/8.02e-5


y = (N-Nb/R)/Ao

def likelihood(psi):

    def func_likelihood(A):
        sum = 0
        for k in range(0, N+1):
            r = (np.math.factorial(N-k+Nb)/(np.math.factorial(k)*np.math.factorial(N-k))*(A*psi*(1+R))**k)
            sum = r+sum
        return sum* (np.exp(-0.5*(np.log(A/Ao)/S_norm)**2)*np.exp(-A*psi)/A)

    return (scipy.integrate.quad(func_likelihood, 0, np.inf, full_output=1,epsabs=1e-300, epsrel=1e-300))[0]

psi = y-y/10

MLE = scipy.optimize.fmin(lambda psi: -likelihood(psi), y, disp=0,full_output=1)[0][0]

normal_factor = 1/likelihood(MLE)

print(normal_factor* likelihood(psi))

Integer division may help. Eg, `(np.math.factorial(N-k+Nb) // (np.math.factorial(k)*np.math.factorial(N-k))` or `(np.math.factorial(N-k+Nb) // (np.math.factorial(k) // np.math.factorial(N-k))` . BTW, you shouldn't use `sum` as a variable name because that shadows the built-in `sum` function. — PM 2Ring, May 31 '17 at 16:23
There should also be an optimized binomial function available in scipy or numpy. — Lutz Lehmann, May 31 '17 at 16:24
@PM2Ring, integer division does not give a correct result. Thanks for your advice regarding using 'sum' as a variable. — DPdl, May 31 '17 at 16:39
Oh well, it was a try. ;) Did you try rearranging the terms, eg `(factorial(N-k+Nb) // factorial(N-k)) / factorial(k)`? Another option is to perform the big multiplication in a loop over `range(N-k+Nb, N-k, -1)`, then find the GCD of that result and `factorial(k)` so you can cancel out the common factor with integer division before attempting the float division. — PM 2Ring, May 31 '17 at 16:45
@PM2Ring, I tried the first option -- still gives me the same error. I will have to try the second option you give. Thanks for your help. — DPdl, May 31 '17 at 16:51
Note that the spelling `np.math.factorial` is a bit misleading, since it has nothing to do with NumPy: it's simply the regular standard library `math.factorial`. (`np.math is math` -> `True`) — Mark Dickinson, May 31 '17 at 18:39

score 2 · Accepted Answer · answered May 31 '17 at 19:04

It is usually more efficient and has less overflow problems to compute the quotient of the terms and use that quotient to update the terms in each step.

Compressed the term of index k is

r[k] = ( (N-k+Nb)! )/( k!*(N-k)! ) * (A*psi*(1+R))**k

so that the quotient to the last term is

r[k+1] / r[k] = ( (N-k) )/( (N-k+Nb)*(k+1) ) * (A*psi*(1+R))

and

r[0] = ( (N+Nb)! )/( N! ) = (N+1)*...*(N+Nb)

which then allows to reformulate the summation function as

 def func_likelihood(A):
      sum = 0
      r = 1
      for k in range(Nb): r *= N+k+1
      sum = r
      for k in range(0, N+1):
          sum += r;
          r *= (A*psi*(1+R)) * (N-k)/((N-k+Nb)*(k+1))
      return sum * (np.exp(-0.5*(np.log(A/Ao)/S_norm)**2)*np.exp(-A*psi)/A)

Prune · Answer 2 · 2017-05-31T16:26:40.543

1

You failed to supply the error message, but I'm betting it's in computing r from the factorials within the for k loop.

There are several ways to get around this. One is to switch the order of operations, such that you eliminate common factors, rather than starting with 58! (which exceeds the default integer limit). This involves a bit more coding, or perhaps calling a combined factorial function, such as one the will do C(n, k) -- @LutzL mentions the binomial function. Granted, you aren't doing quite a binomial, but you could use it and then adjust the numerator as needed for the -k+Nb factors.

Another way is to switch to large integers (see documentation for numpy & scipy, so that you solve the problem in all packages).

edited May 31 '17 at 16:26

answered May 31 '17 at 16:23

Prune

76,765
14
60
81

`np.math.factorial` already uses large integers, and hopefully those factorials can be divided safely. But I agree that using a proper `C(n, k)` is better. – PM 2Ring May 31 '17 at 16:26
Perhaps -- but not handling 58! is a suspicious cut-off. – Prune May 31 '17 at 16:27
1

True, but that's because using `/` on the large integers produces a float result. – PM 2Ring May 31 '17 at 16:28

OverFlowError during exponentiation

2 Answers2