最大似然估计伪代码

Posted 2023-02-19

【中文标题】最大似然估计伪代码

【英文标题】：Maximum Likelihood Estimate pseudocode

【发布时间】：2011-12-04 19:19:41

【问题描述】：

我需要编写一个最大似然估计器来估计一些玩具数据的均值和方差。我有一个包含 100 个样本的向量，使用 numpy.random.randn(100) 创建。数据应具有零均值和单位方差高斯分布。

我查看了***和一些额外的资源，但我有点困惑，因为我没有统计背景。

是否有最大似然估计器的伪代码？我得到了 MLE 的直觉，但我不知道从哪里开始编码。

Wiki 说采用对数似然的 argmax。我的理解是：我需要通过使用不同的参数来计算对数似然，然后我将采用给出最大概率的参数。我没有得到的是：我首先在哪里可以找到参数？如果我随机尝试不同的均值和方差以获得高概率，我应该什么时候停止尝试？

【问题讨论】：

如果你有“数据”，那么均值 = 数据，方差 = 0.0

很抱歉，数据是一个包含 100 个样本的向量。

【参考方案1】：

如果您进行最大似然计算，您需要采取的第一步如下：假设分布取决于某些参数。既然你generate 你的数据（你甚至知道你的参数），你“告诉”你的程序假设高斯分布。但是，您不会告诉程序您的参数（0 和 1），而是让它们先验未知，然后再计算它们。

现在，您有了样本向量（我们称之为x，它的元素是x[0] 到x[100]），您必须对其进行处理。为此，您必须计算以下内容（f 表示 probability density function of the Gaussian distribution）：

f(x[0]) * ... * f(x[100])

正如您在我给定的链接中看到的那样，f 使用两个参数（希腊字母 µ 和 σ）。 您现在必须以f(x[0]) * ... * f(x[100]) 取最大可能值的方式计算 µ 和 σ 的值。

完成后，µ 是平均值的最大似然值，σ 是标准差的最大似然值。

请注意，我没有明确告诉您如何计算 µ 和 σ 的值，因为这是一个我手头没有的非常数学的过程（我可能不明白它）;我只是告诉你获取值的技术，它也可以应用于任何其他分布。

由于您想最大化原始项，您可以“简单地”最大化原始项的对数 - 这样您就无需处理所有这些乘积，并将原始项转换为带有一些加法的总和。

如果你真的想计算它，你可以做一些简化，得出以下术语（希望我没有搞砸任何东西）：

现在，您必须找到 µ 和 σ 的值，以使上述野兽最大。这样做是一项非常重要的任务，称为非线性优化。

您可以尝试的一种简化如下：修复一个参数并尝试计算另一个参数。这使您不必同时处理两个变量。

【讨论】：

感谢您的回答。我的理解是：如果我保持一个参数固定并计算另一个参数，反之亦然，我实际上会做期望最大化算法，对吗？ en.wikipedia.org/wiki/…

我认为可能是这种情况（但我对此不确定）。我认为从平均值（平均值）作为 µ 的起点（将 µ 固定到平均值）然后最大化 σ 可能是一个好的开始...

@Kyle: 也许en.wikipedia.org/wiki/… 是你感兴趣的……

@Kyle FYI 高斯的 MLE 都可以通过分析获得。它们是样本均值和样本方差，尽管后者对于小样本量略有偏差，因此通常除以 n-1 而不是 n。更一般地，您需要学习牛顿法，也许还需要学习 EM（期望最大化）。

@joran 实际上，对于高斯分布，如果我取样本均值和样本方差，我将获得数据集的 MLE。实际上，我需要对数据集应用有偏和无偏的 MLE 估计。那么在这种情况下，您知道仅取样本均值和样本方差是否有效？

【参考方案2】：

我刚遇到这个，我知道它很旧，但我希望其他人能从中受益。尽管之前的 cmets 对什么是 ML 优化给出了很好的描述，但没有人给出实现它的伪代码。 Python 在 Scipy 中有一个最小化器可以做到这一点。这是线性回归的伪代码。

# import the packages
import numpy as np
from scipy.optimize import minimize
import scipy.stats as stats
import time

# Set up your x values
x = np.linspace(0, 100, num=100)

# Set up your observed y values with a known slope (2.4), intercept (5), and sd (4)
yObs = 5 + 2.4*x + np.random.normal(0, 4, 100)

# Define the likelihood function where params is a list of initial parameter estimates
def regressLL(params):
    # Resave the initial parameter guesses
    b0 = params[0]
    b1 = params[1]
    sd = params[2]

    # Calculate the predicted values from the initial parameter guesses
    yPred = b0 + b1*x

    # Calculate the negative log-likelihood as the negative sum of the log of a normal
    # PDF where the observed values are normally distributed around the mean (yPred)
    # with a standard deviation of sd
    logLik = -np.sum( stats.norm.logpdf(yObs, loc=yPred, scale=sd) )

    # Tell the function to return the NLL (this is what will be minimized)
    return(logLik)

# Make a list of initial parameter guesses (b0, b1, sd)    
initParams = [1, 1, 1]

# Run the minimizer
results = minimize(regressLL, initParams, method='nelder-mead')

# Print the results. They should be really close to your actual values
print results.x

这对我很有用。当然，这只是基础知识。它没有描述或给出参数估计的 CI，但它是一个开始。您还可以使用 ML 技术来查找例如 ODE 和其他模型的估计值，正如我所描述的 here。

我知道这个问题很老了，希望你从那时起就明白了，但希望其他人会从中受益。

【讨论】：

不是斜率和y-int分别(1.0101010101010102, 1.0)吗？

您在此处发布的链接已被删除。

【参考方案3】：

正如 joran 所说，正态分布的最大似然估计可以通过分析计算。答案是通过找到对数似然函数关于参数的偏导数，将每个参数设置为零，然后同时求解两个方程来找到的。

在正态分布的情况下，您将导出关于均值 (mu) 的对数似然，然后导出关于方差 (sigma^2) 的两个方程均等于零。求解 mu 和 sigma^2 的方程后，您将得到样本均值和样本方差作为答案。

有关详细信息，请参阅wikipedia page。

【参考方案4】：

您需要一个数值优化程序。不确定是否在 Python 中实现了任何东西，但如果是，那么它将在 numpy 或 scipy 和朋友中实现。

寻找“Nelder-Mead 算法”或“BFGS”之类的内容。如果一切都失败了，请使用 Rpy 并调用 R 函数 'optim()'。

这些函数通过搜索函数空间并尝试找出最大值在哪里来工作。想象一下在雾中寻找山顶。您可能只是尝试始终以最陡峭的方式前进。或者你可以派一些朋友带着收音机和 GPS 装置出去做一些调查。任何一种方法都可能导致您进入错误的顶峰，因此您通常需要从不同的点开始执行几次。否则你可能会认为南峰是最高的，而北峰却是巨大的。

【讨论】：

如果概率密度函数有闭式解，则不必使用数值优化。例如，多变量高斯的参数可以通过 w.r.t. 的导数来评估。 mu 和 sigma 并将其等同于 0。最佳参数对应于数据的 mu 和 sigma。