如何调试scipy中的`OverflowError: math range error`

Question

我有一段测试代码，目的是计算一个最大似然估计（MLE）。

import numpy as np
from scipy.optimize import minimize


def loglikelihood(params, data):
    (mu, alpha, beta) = params
    tlist = np.array(data)
    r = np.zeros(len(tlist)+1)
    for i in xrange(2,len(tlist)):
        r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
    loglik  = data[-1]*mu
    loglik = loglik+alpha/beta*sum(np.exp(-beta*(tlist[-1]-tlist))-1)
    loglik = loglik+np.sum(np.log(mu+alpha*r))
    return -loglik


atimes=[58.98353497,   59.28420225,   59.71571013,   60.06750179,   61.24794134,
   61.70692463,   61.73611983,   62.28593814,   62.51691723,   63.17370423
,   63.20125152,   65.34092403,  214.24934446,  217.0390236,   312.18830525,
  319.38385604,  320.31758188,  323.50201334,  323.76801537,  323.9417007]

print minimize(loglikelihood, (0.01, 0.5,0.6), (atimes,))

在执行 res = minimize(loglikelihood, (0.01, 0.5,0.6), args = (atimes,)) 时，我遇到了 OverflowError: math range error 的错误。

我该如何解决这个问题呢？

---

我只是想把下面的最大似然估计（MLE）R语言代码翻译成Python，并用我放在atimes里的数据进行测试。

neg.loglik <- function(params, data, opt=TRUE) {
  mu <- params[1]
  alpha <- params[2]
  beta <- params[3]
  t <- sort(data)
  r <- rep(0,length(t))
  for(i in 2:length(t)) {
    r[i] <- exp(-beta*(t[i]-t[i-1]))*(1+r[i-1])
  }
  loglik <- -tail(t,1)*mu
  loglik <- loglik+alpha/beta*sum(exp(-beta*(tail(t,1)-t))-1)
  loglik <- loglik+sum(log(mu+alpha*r))
  if(!opt) {
    return(list(negloglik=-loglik, mu=mu, alpha=alpha, beta=beta, t=t,
                r=r))
  }
  else {
    return(-loglik)
  }
}

# insert your values for (mu, alpha, beta) in par
# insert your times for data
opt <- optim(par=c(1,2,3), fn=neg.loglik, data=data)

---

我尝试把Python代码改成返回loglik而不是-loglik。虽然不太清楚为什么R代码是要最大化，但我在http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html上没有找到明确的解释。

现在我得到了

RuntimeWarning: invalid value encountered in log
loglik = loglik+np.sum(np.log(mu+alpha*r))

Answer 1

你对R语言函数的翻译有错误，主要有两个问题：(i) 一是索引从1开始和从0开始的区别，(ii) 二是一些符号错误。

因此，这个函数可能没有下限，也就是说没有最小值。优化器会试着在函数值下降的方向上，使用越来越大的参数值，结果就把很大的数字放进了math.exp，导致溢出。

这个溢出错误不是来自scipy，而是来自你自己写的函数。