关于协方差最小化 scipy.stats.multivariate_normal.logpdf

Posted 2023-03-12

【中文标题】关于协方差最小化 scipy.stats.multivariate_normal.logpdf

【英文标题】：Minimizing scipy.stats.multivariate_normal.logpdf with respect to covariance

【发布时间】：2021-11-15 23:32:07

【问题描述】：

我有一个 python 脚本，我使用 scipy 的 multivariate_normal.log_pdf 计算双变量数据样本的正常对数似然函数的值。我假设样本均值和方差的值，只留下变量之间的样本相关性作为未知数，

from scipy.stats import multivariate_normal
from scipy.optimize import minimize

VAR_X = 0.4
VAR_Y = 0.32
MEAN_X = 1
MEAN_Y = 1.2

def log_likelihood_function(x, data):
    log_likelihood = 0
    sigma = [ [VAR_X, x[0]], [x[0], VAR_Y] ]
    mu = [ MEAN_X, MEAN_Y ]
    for point in data:
        log_likelihood += multivariate_normal.logpdf(x=point, mean=mu, cov=sigma)
    return log_likelihood

if __name__ == "__main__":
    some_data = [ [1.1, 2.0], [1.2, 1.9], [0.8, 0.2], [0.7, 1.3] ]
    
    guess = [ 0 ] 

    # maximize log-likelihood by minimizing the negative 
    likelihood = lambda x: (-1)*log_likelihood_function(x, some_data)
    
    result = minimize(fun = likelihood, x0 = guess, options = 'disp': True, method="SLSQP")

    print(result)

无论我的猜测是什么，这个脚本都会可靠地抛出一个 ValueError，

ValueError: the input matrix must be positive semidefinite

现在，据我估计，问题似乎是scipy.optimize.minimize 是在猜测创建一个非正定协方差矩阵的值。所以我需要一种方法来确保最小化算法会丢弃问题域之外的值。我想为最小化调用添加一个约束，

## make the determinant always positive
def positive_definite_constraint(x):
    return VAR_X*VAR_Y - x*x

这基本上是协方差矩阵的Slyvester Criteron，并确保矩阵是正定的（因为我们知道方差总是正的，所以不需要检查该条件）但似乎scipy.optimize.minimize 评估目标在确定是否满足约束之前运行函数（这似乎是一个设计缺陷；在受限域中搜索解决方案不是更快，而不是搜索所有可能的解决方案然后确定是否满足约束吗？我可能不过，评估的顺序是错误的。）

我不确定如何继续。我意识到我在这里通过参数化协方差矩阵然后最小化该参数化来扩展scipy.optimize的目的，我知道有更好的方法来计算正常样本的相关性，但我对此感兴趣问题，因为它泛化到非正态分布。

有什么建议吗？有没有更好的方法来解决这个问题？

【参考方案1】：

你在正确的轨道上。请注意，您的确定性约束简化为优化变量的简单界限，即-∞ <= x[0] <= VAR_X*VAR_Y。变量边界比更一般的约束在内部处理得更好，所以我推荐这样的东西：

bounds = [(None, VAR_X*VAR_Y)]
res = minimize(fun = likelihood, x0 = guess, bounds=bounds, options = 'disp': True, method="SLSQP")

这给了我：

     fun: 6.610504611834715
     jac: array([-0.0063166])
 message: 'Optimization terminated successfully'
    nfev: 9
     nit: 4
    njev: 4
  status: 0
 success: True
       x: array([0.12090069])

【讨论】：

是的，工作就像一个魅力。你知道我对目标函数和约束之间的评估顺序是否正确吗？

@GrantMoore 是的，SLSQP 方法确实如此，请参阅here。