使用scipy中的minimize函数时出现日志函数运行时错误，如何解决？

Question

为了让大家更明白，我在用NumPy这个库的时候，遇到了一个警告信息。

    RuntimeWarning: invalid value encountered in log
  model = ((-2*r) * np.log(1-s*np.sin(2*np.pi*t/p_orb-phi_orb))) + (-(R**2)/(c*D*9.461E15)*(np.cos(Beta)**2)*np.cos((4*np.pi*t)/(3.154E7)-2*Lambda))

这是我用的代码：

r = 9.85E-6
i = np.pi/2
s = np.sin(i)
N = 157788
t = np.linspace(0 , 9.461E7 , N)
p_orb = 2400
phi_orb = np.pi/3
R = 1.496E11
c = 299752498
D = 1000
Beta = 50*np.pi/180
Lambda = np.pi/4
sigmas = np.ones(N)*0.000001
data = ((-2*r) * np.log(1-s*np.sin(2*np.pi*t/p_orb-phi_orb))) + (-(R**2)/(c*D*9.461E15)*(np.cos(Beta)**2)*np.cos((4*np.pi*t)/(3.154E7)-2*Lambda)) + sigmas

def log_likelihood(theta , t , data , sigmas):
    r , s , D = theta
    model = ((-2*r) * np.log(1-s*np.sin(2*np.pi*t/p_orb-phi_orb))) + (-(R**2)/(c*D*9.461E15)*(np.cos(Beta)**2)*np.cos((4*np.pi*t)/(3.154E7)-2*Lambda))
    return -0.5 * sum(((data - model)**2)/(sigmas**2))

nll = lambda *args: -log_likelihood(*args)
initial = np.array([r , s , D])
soln = minimize(nll , initial , args = (t , data , sigmas))
r_ml, s_ml, D_ml = soln.x
print(r_ml , s_ml , D_ml)

在我的Jupyter笔记本里运行这段代码时，除了我给最小化命令的初始值外，还出现了这个警告信息。我用下面的代码检查了一下有没有零或者负数：

array = 1-s*np.sin(2*np.pi*t/p_orb-phi_orb)
# Check for zeros
zeros_indices = array == 0
zeros_exist = np.any(zeros_indices)

# Check for negative values
negative_indices = array < 0
negatives_exist = np.any(negative_indices)

if zeros_exist:
    print("Array contains zero(s).")
else:
    print("Array does not contain any zero.")

if negatives_exist:
    print("Array contains negative value(s).")
else:
    print("Array does not contain any negative value.")

这里的数组是对数函数的参数。看起来没有负数或零。我不知道还应该检查什么。

Answer 1

问题在于，s是你优化的变量之一，所以minimize()这个函数可以随意尝试不同的值来计算像log(1-s.sin(at-b))这样的表达式。它几乎肯定会尝试比预期答案1更大的或更小的值。

一旦尝试的s值超过1，那么几乎可以肯定（这取决于t[]数组的内容），对数函数的参数会变成负数（这是不合法的）。

你可以通过使用bounds=选项来限制s的范围，从而在运行代码时停止警告。因此，下面的代码会“有效”，也就是说它能恢复你原来的数据参数（包括s=1）：

soln = minimize(nll, initial, args=(t,data,sigmas), bounds=((None,None),(0.0,1.0),(None,None)) )

这会生成(r,S,D)，结果是

9.93687622297525e-06 1.0 1000.0

不过，这更多是运气而不是判断。你的离散t值并没有产生一个sin()等于1的值，所以log(1-s.sin())从未失败。

我建议你不要这样尝试最小化，特别是使用从s=sin(pi/2)=1生成的模型数据，并且一定要给你的优化参数设置正式的范围，以防止像log()这样的函数出现无效的参数。