Python中多元正态分布的最大似然估计

Question

我正在尝试使用MLE拟合多元正态分布的参数。

import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize
from scipy.stats import multivariate_normal as mnorm

def estimation(obs,fun,init,method='Nelder-Mead'):
    mle = lambda param: -np.sum(fun(*[obs,param])) ## negate since we will minimize
    result = minimize(mle,init,method=method)
    return result.x

拟合单变量正态分布是可以的：

obs = np.random.normal(1,4,50000)
ini = [0,1]
print(estimation(obs,lambda ob,p:norm.logpdf(ob,p[0],p[1]),ini))

但在多元方面遇到了一些问题（分配数组给变量时出错）：

obs_m = np.random.multivariate_normal([0,0],[[1,0],[0,100]],50000)
ini_m = [[0,0],[[1,0],[0,100]]]
print(estimation(obs_m,lambda ob,p:mnorm.logpdf(ob,p[0],p[1],ini_m))

看起来优化算法不适用于任意数组/矩阵。我必须将均值数组和协方差矩阵解包成一个平坦的数组，然后传递给 minimize 函数。

ini_m = [0,0,1,0,0,100]
print(estimation(obs_m,lambda ob,p:mnorm.logpdf(ob,[p[0],p[1]],[[p[2],p[3]],[p[4],p[5]]]),ini_m))

显然，当维度增加或某些更复杂的分布没有封闭形式解时，这种情况很快就会失控。在这里应该怎么办呢？谢谢。

- iwbabn

Answer 1

假设optimize.minimize允许您在不进行修改的情况下传入mean和cov的猜测值，无论其形状如何，并且它将传入与您的负对数似然函数相同形状的mean和cov数组。由于cov必须是对称的，因此一些cov的元素将是多余的，并且您需要以某种方式添加约束条件以保持cov正定。

下面的解决方案需要编写特定于分布的函数，用于在决策变量的1-D数组和特定于分布的参数之间进行转换，但它适用于任意维数，并且对于multivariate_normal，它解决了冗余决策变量和确保正定性的问题。

import numpy as np
from scipy import stats, optimize
rng = np.random.default_rng(747458538356)

def mle(obs, dist, params0, method='slsqp'):
    # Distribution-independent MLE function
    dim = obs.shape[-1]

    def nllf(x):
        params = x2params(x, dim)
        logpdf = dist.logpdf(obs, *params)
        return -np.sum(logpdf)

    x0 = params2x(params0)
    res = optimize.minimize(nllf, x0, method=method)
    return x2params(res.x, dim)


dist = stats.multivariate_normal

def x2params(x, dim):
    # multivariate_normal-specific function to convert decision variable
    # array to distribution parameters
    mean = x[:dim]

    # Decision variables are elements of Cholesky decomposition
    # to ensure positive definite covariance matrix
    L = np.zeros((dim, dim))
    indices = np.tril_indices(dim)
    L[indices] = x[dim:]
    cov = L @ L.T

    return mean, cov

def params2x(params):
    # multivariate_normal-specific function to convert distribution parameters
    # to decision variable array
    mean, cov = params
    dim = len(mean)
    indices = np.tril_indices(dim)
    L = np.linalg.cholesky(cov)
    return np.concatenate((np.ravel(mean), L[indices]))


# generate observations
size = 500  # number of observations
mean = [0, 0]
cov = np.array([[5, 2], [2, 10]])
params0 = (mean, cov)
obs = dist.rvs(mean, cov, size=size, random_state=rng)

# compare our function against reference
res = mle(obs, dist, params0)
ref = stats.multivariate_normal.fit(obs)
np.testing.assert_allclose(res[0], ref[0], rtol=1e-3)
np.testing.assert_allclose(res[1], ref[1], rtol=1e-3)

mean, cov = res

print(mean)
# [0.09533408 0.04011049]

print(cov)
# [[ 4.91358864  2.25704525]
#  [ 2.25704525 10.59167997]]

- Matt Haberland