title: OverFlowError during exponentiation
body: <p>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.</p>
  
  <pre><code>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))
  </code></pre>