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

ie3xauqp 于 10个月前发布在 Python

关注(0)|答案(1)|浏览(207)

我试图用极大似然估计拟合多元正态分布的参数。

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

型
看来优化算法对任意数组/矩阵都不起作用。我必须将均值数组和协方差矩阵解压缩成一个平面数组来最小化。

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

型
显然，当维数增加时，或者对于一些没有封闭形式解的更复杂的分布，这会很快失控，在这里最好怎么做呢？谢谢。

scipy

来源：https://stackoverflow.com/questions/42889386/maximum-likelihood-estimation-for-multivariate-normal-distribution-in-python

1条答案

按热度按时间

dtcbnfnu1#

假设optimize.minimize * 确实 * 允许你不加修改地传入mean和cov的猜测值--以你想要的任何形状--它会传入相同形状的mean和cov数组到你的负对数似然函数。cov的一些元素会是冗余的，因为它必须是对称的，你需要增加一些约束来保持cov的正定性。
下面的解决方案需要编写特定于分布的函数，用于在决策变量的一维数组和特定于分布的参数之间进行转换，但它可以扩展到任意维数，并且对于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]]

字符串

赞(0）回复(0）举报 10个月前

我来回答

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

1条答案

相关问题

热门标签

最新问答