假设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]]

字符串