使用scipy minimize时log函数出现RuntimeWarning，该如何修复？

问题描述 #

使用NumPy结合scipy的minimize函数进行参数拟合时触发RuntimeWarning: invalid value encountered in log，对应代码为计算model的语句。流程是先生成模拟数据，定义对数似然函数log_likelihood，通过最小化负对数似然拟合参数，但运行后minimize返回初始参数值，同时出现上述警告。

已检查初始参数下的log函数参数数组1-s*np.sin(2*np.pi*t/p_orb-phi_orb)，确认无零或负值，但警告仍存在。

相关代码 #

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)

检查代码及结果 #

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.")

输出：

Array does not contain any zero.
Array does not contain any negative value.

修复方案 #

1. 给参数添加边界约束 #

minimize在迭代过程中会尝试不同的参数值，可能出现s的绝对值大于1的情况，导致1-s*sin(...)变成非正值。针对参数的物理意义设置边界：

from scipy.optimize import Bounds

# 设置边界：r>0, s∈[-1,1], D>0
bounds = Bounds([0, -1, 0], [np.inf, 1, np.inf])
soln = minimize(nll, initial, args=(t, data, sigmas), bounds=bounds)

L-BFGS-B算法（默认支持边界的算法）会自动限制参数在合法范围内迭代。

2. 提升log计算的数值稳定性 #

即使有边界约束，浮点精度问题可能导致输入接近0，触发警告。可以用np.maximum强制log输入大于极小值：

def log_likelihood(theta , t , data , sigmas):
    r , s , D = theta
    # 确保log输入始终大于等于1e-12
    log_arg = np.maximum(1-s*np.sin(2*np.pi*t/p_orb-phi_orb), 1e-12)
    model = ((-2*r) * np.log(log_arg)) + (-(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))

3. 调试迭代过程中的参数 #

在似然函数中加入调试代码，定位迭代时是否出现非法参数：

def log_likelihood(theta , t , data , sigmas):
    r , s , D = theta
    log_arg = 1-s*np.sin(2*np.pi*t/p_orb-phi_orb)
    # 打印当前参数和log输入的极值
    print(f"Params: r={r:.6e}, s={s:.6f}, D={D:.0f}")
    print(f"log_arg min: {np.min(log_arg):.6e}, max: {np.max(log_arg):.6e}")
    if np.any(log_arg <= 0):
        print("⚠️ Invalid log argument detected!")
    log_arg = np.maximum(log_arg, 1e-12)
    model = ((-2*r) * np.log(log_arg)) + (-(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))

通过输出可以确认约束是否生效，以及问题出现的具体参数组合。

4. 指定合适的优化算法 #

明确指定支持边界约束的算法，避免默认算法切换带来的问题：

soln = minimize(nll, initial, args=(t, data, sigmas), bounds=bounds, method='L-BFGS-B')

内容的提问来源于stack exchange，提问作者Tihomb