最小化scipy.stats.multivariate_normal.关于协方差的Logpdf &g



我有一个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准则,并将确保矩阵是正定的(因为我们知道方差总是正的,该条件不需要检查),但似乎scipy.optimize.minimize在确定约束是否满足之前评估目标函数(这似乎是一个设计缺陷;在受限域中搜索一个解,而不是搜索所有可能的解,然后确定约束是否满足,这样不是更快吗?不过,我可能搞错了求值的顺序。)

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

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

你做对了。请注意,您的确定性约束减少到优化变量的简单边界,即-∞ <= 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])

最新更新