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))

Answer 1

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)