Python 中 MLE（非线性模型）的多变量梯度下降

Question

我正在尝试执行梯度下降来计算三个 MLE（从头开始）。我有数据 $x_i=s_i+w_i$ 其中 $s_i=A(nu_i/nu_0)^{alpha}(nu_i/nu_0+1)^{-4alpha}$ 我已分析计算一阶导数：

import numpy as np
import matplotlib.pyplot as plt

#model function s_i
def signal(A,nu_0,alpha,nu):
    return A*(nu/nu_0)**alpha*(1+nu/nu_0)**(-4*alpha)

#partial derivatives of log likelihood function
def MLE_A(A,nu_0,alpha,nu,x_i):
    nu=np.array(nu)
    x_i=np.array(x_i)
    return -np.sum(((nu/nu_0)**alpha*((A*(nu/nu_0)**alpha)/(nu/nu_0+1)**(4*alpha)-x_i))/(nu/nu_0+1)**(4*alpha))

def MLE_alpha(A,nu_0,alpha,nu,x_i):
    nu=np.array(nu)
    x_i=np.array(x_i)
    
    return -np.sum((A*(nu/nu_0)**alpha*(4*np.log(nu/nu_0+1)-np.log(nu/nu_0))*(x_i*(nu/nu_0+1)**(4*alpha)-A*(nu/nu_0)**alpha))/(nu/nu_0+1)**(8*alpha))

def MLE_nu_0(A,nu_0,alpha,nu,x_i):
    nu=np.array(nu)
    x_i=np.array(x_i)
    
    return -np.sum((A*alpha*(nu/nu_0)**(alpha)*(nu_0-3*nu)*((x_i*((nu)/nu_0+1)**(4*alpha))-A*(nu/nu_0)**alpha))/(nu_0*(nu+nu_0)*((nu)/nu_0+1)**(8*alpha)))

 
def gradient_descent(A_init,nu_0_init,alpha_init,nu,x_i,iterations=1000,learning_rate=0.01):
    
    A=A_init
    nu_0=nu_0_init
    alpha=alpha_init
    theta=np.array([A_init,nu_0_init,alpha_init])
    updated_theta=[theta]
    for i in range(iterations):
        new_theta = theta - learning_rate * np.array([MLE_A(A,nu_0,alpha,nu,x_i), MLE_nu_0(A,nu_0,alpha,nu,x_i), MLE_alpha(A,nu_0,alpha,nu,x_i)])
        theta=new_theta
        updated_theta.append(theta)
        A,nu_0,alpha = new_theta[0],new_theta[1],new_theta[2]
    return(updated_theta)

A=6
nu_0=2
alpha=1
nu=np.linspace(0.05,1.0,200)
x_i=signal(A,nu_0,alpha,nu)+np.random.normal(0,0.05,len(nu))


params= gradient_descent(A,nu_0,alpha,nu,x_i,iterations=10000,learning_rate=0.01)

print(params[-1])
A_fit=params[-1][0]
nu_0_fit=params[-1][1]
alpha_fit=params[-1][2]

plt.plot(nu,x_i)
plt.plot(nu,signal(A_fit,nu_0_fit,alpha_fit,nu))

plt.show()

有时我会收到类似 RuntimeWarning: Overflow returned in power 和 RuntimeWarning: invalid value returned in true_divide 的错误，有时我会得到严重偏离值的错误。我对学习率使用了不同的值，但它没有解决这个问题。我已经手动并使用符号软件检查了这些功能。另外，我使用 Latexify 来查看我确实正确输入了它们，所以我假设这是我对梯度下降的实现在某种程度上关闭了。

非常感谢任何帮助。

Answer 1

观察结果

在评估导数时，您的 MCVE 似乎很难处理某些

xi

值：它会遇到溢出或除以 0。因此，您的参数向量包含

NaN

并且算法失败。

对于某些设置，问题不存在，例如将随机种子固定为：

np.random.seed(1234567)

生成一个防止代码失败的数据集，然后我们得到以下结果：

# [5.66701368 2.18827948 5.48118174]

这可能表明最后的导数不准确。

MCVE

通过对梯度进行数值评估来执行 MSE 梯度下降效果很好（对于其他种子也是如此），并且与模型无关（用户不必提供导数）：

import numpy as np
import numdifftools as nd
from scipy import optimize
import matplotlib.pyplot as plt

def model(x, A, nu0, alpha):
    return A * np.power((x/nu0), 1. * alpha) * np.power((1 + x/nu0), (-4. * alpha))

def gradient_factory(model, x, y):
    def wrapped(p):
        return 0.5 * np.sum(np.power(y - model(x, *p), 2.))
    return nd.Gradient(wrapped)

np.random.seed(1234567)
p0 = (6, 2, 1)
nu = np.linspace(0.05, 1.0, 200)
xi = model(nu, *p0) + 0.05 * np.random.normal(size=nu.size)

def gradient_descent(model, x, y, p0, tol=1e-16, maxiter=500, rate=0.001, atol=1e-10, rtol=1e-8):

    gradient = gradient_factory(model, x, y)
    
    p = np.array(p0)
    dp = gradient(p)

    for _ in range(maxiter):

        # Update gradient descent:
        p_ = p - rate * dp
        dp_ = gradient(p_)

        # Update rate:
        Dp_ = p_ - p
        Ddp_ = dp_ - dp
        rate = np.abs(Dp_.T @ Ddp_) / np.power(np.linalg.norm(Ddp_), 2)

        # Break when precision is reached:
        if np.allclose(p, p_, atol=atol, rtol=rtol):
            break

        # Next iteration:
        p = p_
        dp = dp_

    else:
        raise RuntimeError("Max Iteration (maxiter=%d) reached" % maxiter)
        
    return p

p = gradient_descent(model, nu, xi, [1., 1., 1.])
# array([5.82733367, 2.06411618, 0.98227514])

并且收敛到与

curve_fit

相同的解决方案：

popt, pcov = optimize.curve_fit(model, nu, xi)
# (array([5.82905928, 2.06399904, 0.98240413]),
#  array([[ 0.56149129, -0.03771281,  0.04195576],
#         [-0.03771281,  0.00493833, -0.00287965],
#         [ 0.04195576, -0.00287965,  0.00314415]]))

您可以调整此 MCVE 来执行 MLE，而不是最小化 MSE。您还可以提供自己的梯度定义来驱动下降，而不是用数字来评估它。