scipy没有优化并返回"由于精度损失而无法实现所需的错误"

Question

我有以下代码尝试最小化对数似然函数.

#!/usr/bin/python
import math
import random
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))
    for i in xrange(1,len(tlist)):
        r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
    loglik  = -tlist[-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 = [ 148.98894201,  149.70253172,  151.13717804,  160.35968355,
        160.98322609,  161.21331798,  163.60755544,  163.68994973,
        164.26131871,  228.79436067]
a= 0.01
alpha = 0.5
beta = 0.6
print loglikelihood((a, alpha, beta), atimes)

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'BFGS',args = (atimes,))
print res

它给了我

28.3136498357
./test.py:17: RuntimeWarning: invalid value encountered in log
  loglik = loglik+np.sum(np.log(mu+alpha*r))
   status: 2
  success: False
     njev: 14
     nfev: 72
 hess_inv: array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]])
      fun: 32.131359359964378
        x: array([ 0.01,  0.1 ,  0.1 ])
  message: 'Desired error not necessarily achieved due to precision loss.'
      jac: array([ -2.8051672 ,  13.06962156, -48.97879982])

请注意,它根本没有设法优化参数,最小值32大于28,这是你得到的a = 0.01,alpha = 0.5,beta = 0.6.通过选择更好的初始猜测可以避免这个问题,但如果是这样,我该如何自动完成？

Answer 1

我复制了你的例子并试了一下.看起来如果你坚持使用BFGS求解器,经过几次迭代后mu+ alpha * r会有一些负数,这就是你获得RuntimeWarning的方法.

我能想到的最简单的解决方法是切换到Nelder Mead求解器.

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'Nelder-Mead',args = (atimes,))

它会给你这个结果:

28.3136498357
  status: 0
    nfev: 159
 success: True
     fun: 27.982451280648817
       x: array([ 0.01410906,  0.68346023,  0.90837568])
 message: 'Optimization terminated successfully.'
     nit: 92

Answer 2

另一个解决方案（对我有用）是将您的函数（和梯度）缩放到接近 0 的值。例如，当我必须评估 60k 点的对数似然时，我的问题就出现了。这意味着我的对数似然是一个非常大的数字。从概念上讲，对数似然是一个非常非常尖峰的函数。

梯度开始很大（攀登这座尖峰山），然后变得适度小，但绝不会小于gtolBGFS 例程中的默认参数（这是所有梯度必须低于终止的阈值）。此外，此时我基本上已经得出了正确的值（我使用了生成的数据，所以我知道了真实值）。

发生的事情是我的梯度大约是。60k * average individual gradient value，即使average individual gradient value很小，也可以说小于 1e-8， 60k * 1e-8 > gtol。因此，即使我已经找到了解决方案，我也从未满足过阈值。

从概念上讲，由于这座山峰非常陡峭，算法正在迈出小步，但跨越了真正的最小值并且从未达到过average individual gradient << 1e-8，这意味着我的梯度从未下降过gtol。

两种解决方案：

1）按一个因子缩放您的对数似然和梯度，例如样本数量1/n在哪里n。

2）扩展你的gtol：例如"gtol": 1e-7 * n

Answer 3

留意 log() 函数的负值，解决它们并通过添加惩罚来告诉优化器它们是坏的：

#!/usr/bin/python
import math
import random
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))
    for i in xrange(1,len(tlist)):
        r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
    loglik = -tlist[-1]*mu
    loglik += alpha/beta*sum(np.exp(-beta*(tlist[-1]-tlist))-1)
    argument = mu + alpha * r
    limit = 1e-6
    if np.min(argument) < limit:
        # add a penalty for too small argument of log
        loglik += np.sum(np.minimum(0.0, argument - limit)) / limit
        # keep argument of log above the limit
        argument = np.maximum(argument, limit)
    loglik += np.sum(np.log(argument))
    return -loglik

atimes = [ 148.98894201,  149.70253172,  151.13717804,  160.35968355,
        160.98322609,  161.21331798,  163.60755544,  163.68994973,
        164.26131871,  228.79436067]
a= 0.01
alpha = 0.5
beta = 0.6
print loglikelihood((a, alpha, beta), atimes)

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'BFGS',args = (atimes,))
print res

Answer 4

面对同样的警告，我通过重写对数似然函数来解决它，以获取log(params)和log(data)作为参数，而不是参数和数据。

因此，np.log()如果可能，我避免在似然函数或雅可比矩阵中使用。