如何在 scipy 中调试“溢出错误：数学范围错误”

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

它给了我OverflowError: math range errorres = minimize(loglikelihood, (0.01, 0.5,0.6), args = (atimes,))。

我该如何解决？

---

我只是试图将以下最大似然估计 (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。尽管 http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html 我不清楚为什么会这样，但 R 代码必须最大化。

现在我明白了

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，而是来自你自己的函数。