Stata：最大似然估计(MLE)简易教程

发布时间：2022-12-10 阅读 1870

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

课程详情： https://gitee.com/lianxh/Course

课程主页： https://gitee.com/lianxh/Course

Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

PDF下载 - 推文合集

编译：王萃芳 (东北财经大学)
邮箱：wcfsophia@26.com

编者按：本文主要整理自 Maximum Likelihood Estimation (MLE)，特此致谢！

本文主要是介绍最大似然估计 (MLE) 的一些基础知识，并演示如何在 Stata 软件中进行最大似然估计。

1. MLE 简介

MLE 的本质是参数估计，最合理的参数估计量应该使得从模型中抽取 $n$ 组样本观测值的概率最大，也就是概率分布函数或者说是似然函数最大。还有一些非参数 MLE，其分布函数来自训练数据集或与我们的数据非常相似的其他数据。例如，如果有一百万张猫的图片，我们可以使用它们对我们自己图片中的猫进行分类。

2.最大似然估计的简单例子

让我们从抛硬币的经典例子开始。抛硬币有二种结果：正面 ( $H$ ) 和反面 ( $T$ )。每次抛硬币的初始条件相同，即硬币不变，每次抛硬币不相互影响。在统计学中，这些属性意味着抛硬币是独立同分布。

假设得到正面 ( $H$ ) 的概率等于某个未知值 $p$ 或者 $P r (H) = p$ 。那么，得到反面的概率，即没有得到正面的概率 ( $H^{'}$ 或者 $1 - H$ ) 为 $1 - p$ ，公式为 $P r (T) = P r (1 - H) = 1 - p$ 。

如何利用 MLE 计算 $p$ 的值呢？我们抛 5 次硬币，得到如下序列：H T T H H。那么从我们观察到的数据中得到正面 H 的可能性是 ( $L$ ) 是多少，即 $L (p)$ 是多少？

假设概率是通过我们观察到的 5 次回合中正面的联合概率分布得出。从统计学上我们知道联合概率分布等于结果 1 的概率乘以结果 2 的概率乘以所有其他结果的概率。因此，我们的问题的似然函数可以写成：

L (p) = p \cdot (1 - p) \cdot (1 - p) \cdot p \cdot p

或者是上述情况可以简化为：

L (p) = p^{3} (1 - p)^{2}

其中，3 是 $H$ 出现的次数，2 是 $T$ 出现的次数。简单来说，如果我们总共有 $N$ 次抛掷，其中 $n$ 次是正面，那么我们将上述问题简化为以下通用函数形式：

L (p) = p^{n} (1 - p)^{N - n}

这是一个有未知数 $p$ 的函数，我们将 $L$ 关于 $p$ 微分并令其等于 0 以找到 $p$ ( $L (p) = d L / d p = 0$ ) 的最优解。因此，我们会得到如下表达式：

\frac{d L (p)}{d p} = n p^{n - 1} (1 - p)^{N - n} - (N - n) (1 - p)^{N - n - 1} p^{n}

接下来，对等式两边同时取对数，将似然函数转化为对数似然函数 (LL)。

L L (p) = n \ln (p) + (N - n) \ln (1 - p)

对 $p$ 求导，并求得 $p$ ：

\begin{aligned} \frac{d L L (p)}{d p} & = \frac{n}{p} - \frac{N - n}{1 - p} = 0 \\ p & = \frac{n}{N} \end{aligned}

实际上，从我们的示例中，由于我们知道 $n = 3$ (3 次正面) 和 $N = 5$ (5 次抛掷)，因此概率 $p$ 等于 $p = 3 / 5 = 0.6$ 。如果我们开始增加投掷次数 ( $N$ 变大)，那么如果硬币是公平的，概率将收敛到 0.5。

3.MLE 的 Stata 实现

3.1 范例 1：抛硬币例子的 MLE 估计

我们先生成一个变量表示投掷硬币的情况。由于我们知道对数似然 (LL) 函数等于 $L L (p) = 3 \times l o g (p) + 2 \times l o g (1 - p)$ ，因此可以绘制下图：

. clear all
. set obs 5
. gen n = _n
. gen x = .
. replace x = 1 in 1
. replace x = 0 in 2
. replace x = 0 in 3
. replace x = 1 in 4
. replace x = 1 in 5
. twoway (function 3 * log(x) + 2 * log(1-x), range(0 1)), ytitle("LL")

可以看到 LL 在 0.6 附近有一个最大值，我们可以利用 Mata 中的 optimize 函数，来找到最大值。

mata
void myfunc(todo, x, y, g, H) {
     y = 3 * log(x) + 2 * log(1-x)
}
maxval = optimize_init()
optimize_init_which(maxval, "max")
optimize_init_evaluator(maxval, &myfunc())
optimize_init_params(maxval, 0.2)
xmax = optimize(maxval)
xmax
end

我们在这定义了一个最大化的函数，并以任意值 0.2 开始搜索，我们可以看到它很接近图中最大值点。对于这样一个简单的优化问题，起始值是无关紧要的。通过以下代码也可以绘制上图：

mata
void myfunc(todo, x, y, g, H) {
 y = 3 * log(x) + 2 * log(1-x)
 }
maxval = optimize_init()
optimize_init_which(maxval, "max")
optimize_init_evaluator(maxval, &myfunc())
optimize_init_params(maxval, 0)
xmax = optimize(maxval)
xmax
ymax = 3 * log(xmax) + 2 * log(1-xmax)
st_local("maxx", strofreal(xmax))
st_local("maxy", strofreal(ymax))
end

twoway (function 3 * log(x) + 2 * log(1-x), range(0 1)),      ///
     yline(`maxy') xline(`maxx') ytitle("LL") xlabel(0(0.2)1)

从中我们得到 LL 的最大值为 -3.365。但是由于 Stata 具有内置的 MLE 函数 (请参阅参考资料help ml) 和其他学习文档，因此我们不需要费劲地在 Mata 中定义优化问题。我们可以简单地使用该 ml 功能。这个命令的使用有点复杂，为了使用它，我们需要编写小程序，在其中定义需要优化的 LL 函数。在抛硬币的情况下，我们已经知道形式上是线性的 LL 函数，我们可以使用 mlexp 选项 (请参阅参考资料help mlexp)。

. mlexp ((x * log({p})) + ((1 - x) * log(1 - {p})))


Maximum likelihood estimation
Log likelihood = -3.3650583                                  Number of obs = 5
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
          /p |      0.600      0.219     2.74   0.006        0.171       1.029
------------------------------------------------------------------------------

在上面的命令中，x 是我们生成的 Stata 变量，{p} 是我们需要最大化的值，该值应始终写在大括号中。也可以有多个参数，每个参数都放在自己的大括号中。 $p$ 值被正确地估计为 0.6。标准误差是由有限样本的渐近方差导出的，在我们的例子中，它是非常小的。如果我们增加投掷硬币的次数，那么 $p$ 应该接近 0.5，标准误差应该接近 0，这是总体分布的真实值。我们可以做如下检查：

. clear all
. set obs 1000000         // or higher or lower
. gen x = uniform() < .5  // a 0/1 variable
. mlexp ((x * log({p})) + ((1 - x) * log(1 - {p})))

Maximum likelihood estimation
Log likelihood = -693146.23                          Number of obs = 1,000,000
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
          /p |      0.499      0.000   998.62   0.000        0.498       0.500
------------------------------------------------------------------------------

mlexp 命令还有几个估计选项，用于执行各种测试、参数的线性组合、创建边距图等。有关详细信息，请参阅 help mlexp postestimate。

3.2 范例 2：Logit 模型的 MLE 估计

MLE 的另一个常见示例是 Logit 或 Probit 回归。这是处理二分变量的标准工具。Logit 的可能性来自逻辑分布，它采用以下函数形式：

L (z) = \prod_{n = 1}^{N} {(\frac{e^{z_{n}}}{1 + e^{z_{n}}})}^{y_{n}} {(\frac{1}{1 + e^{z_{n}}})}^{1 - y_{n}}

虽然这个公式看起来很复杂，但是他遵循与投掷硬币相同的逻辑。这里我们计算 $y = 1$ 的概率：

P r (y = 1) = e x p (z) / (1 + e x p (z))

或者 $y = 0$ 的概率：

P r (y = 0) = 1 - e x p (z) / (1 + e x p (z)) = 1 / (1 + e x p (z))

上面显示的似然公式是二分变量 $y$ 的联合概率分布。Probit 模型的似然性来自正态分布，如下所示：

L (z) = \prod_{n = 1}^{N} {(\int_{- \infty}^{z_{n}} \frac{1}{\sqrt{2 π}} e^{- 0.5 u^{2}} d u)}^{y_{n}} {(1 - \int_{- \infty}^{z_{n}} \frac{1}{\sqrt{2 π}} e^{- 0.5 u^{2}} d u)}^{1 - y_{n}}

在这两个似然分布中， $z = X β$ 是解释变量的矩阵， $n$ 是观察值的数量，从 1 到 N。Logit 模型更适用，因为它更容易处理。我们生成五个观察结果：

. clear all
. set obs 5
. gen y = .
. gen x = .
. replace y = 1 in 1
. replace y = 0 in 2
. replace y = 0 in 3
. replace y = 1 in 4
. replace y = 1 in 5
. replace x = -1 in 1
. replace x =  1 in 2
. replace x = -4 in 3
. replace x =  8 in 4
. replace x = -5 in 5

其中 $y$ 是结果变量， $x$ 是解释变量。可以有更多的解释变量，但由于我们想要绘制似然函数，因此只使用一个解释变量。在模型中，我们还包括一个截距， $y$ 被解释为 $z = α + β x$ 。对于上述观察结果，似然函数可以推导为：

L (z) = \frac{e^{3 α + 12 β}}{(e^{α} + e^{β}) (e^{α} + e^{4 β}) (e^{α} + e^{5 β}) (1 + e^{α + β}) (1 + e^{α + 8 β})}

由于对 $z = α, β$ 的微分和求解很困难，因此我们使用对数似然 (LL) 表达式：

L L (z) = \ln (\frac{1}{1 + e^{α - 4 β}}) + \ln (\frac{e^{α}}{e^{α} + e^{β}}) + \ln (\frac{e^{α}}{e^{α} + e^{5 β}}) + \ln (\frac{1}{1 + e^{α + β}}) + \ln (\frac{e^{α + 8 β}}{1 + e^{α + 8 β}})

这样通用 Logit LL 可以写成：

L L (z) = \sum_{n = 1}^{N} (y_{n}) \ln (\frac{e^{z_{n}}}{1 + e^{z_{n}}}) + (1 - y_{n}) \ln (\frac{1}{1 + e^{z_{n}}})

为了在 Stata 中解决这个问题，我们使用了标准 ml 函数 (详细资料参考 help ml)。对于 Logit 情况，我们需要定义一个如下的程序：

cap prog drop mllogit
program define mllogit
     args lnf Xb
     quietly replace `lnf' = ln(    invlogit(`Xb')) if $ML_y1==1
     quietly replace `lnf' = ln(1 - invlogit(`Xb')) if $ML_y1==0
end

根据上面的代码，ml 参数描述如下。程序第三行用于声明输入项 args，其后紧跟对数似然函数值 lnf 和参数名称 Xb。我们分别指定因变量 $y = 0$ 和 $y = 1$ 的条件。通用术语 $ML_y1 表示解释变量 $y$ 。分布 invlogit 是一个内置函数，其中 invlogit(x) = exp(x)/(1 + exp(x))。另一种指定上述条件的方法是：

program define mllogit2
     args lnf Xb
quietly replace `lnf' =       - ln(1+exp(-`Xb')) if $ML_y1==1
quietly replace `lnf' = -`Xb' - ln(1+exp(-`Xb')) if $ML_y1==0
end

下一步我们需要在 ml 实例中调用这个程序：

. ml model lf mllogit (y = x)
. ml maximize  // or minimize depending on the setup

                                                        Number of obs =      5
                                                        Wald chi2(1)  =   0.25
Log likelihood = -3.2265299                             Prob > chi2   = 0.6169
------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
           x |      0.111      0.221     0.50   0.617       -0.323       0.545
       _cons |      0.460      0.956     0.48   0.630       -1.414       2.334
------------------------------------------------------------------------------

还可以对照标准的 Logit 命令检查这一点：

. logit y x

Logistic regression                                     Number of obs =      5
                                                        LR chi2(1)    =   0.28
                                                        Prob > chi2   = 0.5986
Log likelihood = -3.2265299                             Pseudo R2     = 0.0412
------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
           x |      0.111      0.221     0.50   0.617       -0.323       0.545
       _cons |      0.460      0.956     0.48   0.630       -1.414       2.334
------------------------------------------------------------------------------

3.3 范例 3：Probit 模型的 MLE 估计

除了 Logit 函数，我们还可以使用从正态分布导出的 Probit 函数。Probit 的 LL 就是我们之前定义的可能性的对数。在 Stata 中，可以这样定义 Probit 的 LL：

cap prog drop mlprobit
program mlprobit
args lnf xb
quietly replace `lnf' = ln(    normal(`xb')) if $ML_y1==1
quietly replace `lnf' = ln(1 - normal(`xb')) if $ML_y1==0
end

. ml model lf mlprobit (y = x)
. ml maximize

                                                        Number of obs =      5
                                                        Wald chi2(1)  =   0.27
Log likelihood = -3.2193281                             Prob > chi2   = 0.6036
------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
           x |      0.072      0.139     0.52   0.604       -0.200       0.344
       _cons |      0.292      0.589     0.50   0.620       -0.862       1.446
------------------------------------------------------------------------------

同样也可以与标准 Stata 命令比较：

. probit y x

Probit regression                                       Number of obs =      5
                                                        LR chi2(1)    =   0.29
                                                        Prob > chi2   = 0.5893
Log likelihood = -3.2193281                             Pseudo R2     = 0.0433
------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
           x |      0.072      0.139     0.52   0.604       -0.200       0.344
       _cons |      0.292      0.589     0.50   0.620       -0.862       1.446
------------------------------------------------------------------------------

为了画出比较 Probit 和 Logit 的图表，我们可以使用以下命令:

. probit y x
. global p0 = e(b)[1,2]
. global p1 = e(b)[1,1]
. twoway (function y = exp(($b1 * x) + $b0) / (1 +exp(($b1 * x) + $b0) ), range(-60 60)) ///
>     (function y = normal(($p1 * x) + $p0), range(-60 60)), ytitle("Probability")       ///
>     legend(order(1 "Logit" 2 "Probit"))

然后我们就得到了下图：

在这里，我们可以看到这两个函数的累积密度函数 (CDF) 的差异。

3.4 范例 4：使用实际数据验证

下面让我们使用 Stata 的 union 数据集来验证一下。

. webuse union, clear
. logit union age grade south black
. probit union age grade south black

在 ML 形式中，我们可以使用上面定义的函数:

cap prog drop mllogit
program define mllogit
args lnf Xb
quietly replace `lnf' = ln(    invlogit(`Xb')) if $ML_y1==1
quietly replace `lnf' = ln(1 - invlogit(`Xb')) if $ML_y1==0
end

cap prog drop mlprobit
program mlprobit
args lnf xb
quietly replace `lnf' = ln(    normal(`xb')) if $ML_y1==1
quietly replace `lnf' = ln(1 - normal(`xb')) if $ML_y1==0
end

. ml model lf mllogit (union = age grade south black)
. ml maximize
. ml model lf mlprobit (union = age grade south black)
. ml maximize

如前所述，Stata probit 和 logit 命令利用了 ml 使用的 Mata 函数 optimize。

3.5 范例 5：线性回归模型的 MLE 估计

作者这里用 MLE 处理标准回归，使用 auto 数据集，具体形式如下：

. sysuse auto, clear
. capture program drop myols
. program myols
  1.     args lnf Xb lnsigma
  2.     local y "$ML_y1"
  3.     quietly replace `lnf' = ln(normalden(`y', `Xb', exp(`lnsigma')))
  4. end
. ml model lf myols (xb: mpg = weight foreign) (lnsigma:)
. ml maximize

                                                        Number of obs =     74
                                                        Wald chi2(2)  = 145.39
Log likelihood = -194.18306                             Prob > chi2   = 0.0000
------------------------------------------------------------------------------
         mpg | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
xb           |
      weight |     -0.007      0.001   -10.56   0.000       -0.008      -0.005
     foreign |     -1.650      1.054    -1.57   0.117       -3.716       0.416
       _cons |     41.680      2.121    19.65   0.000       37.522      45.837
-------------+----------------------------------------------------------------
lnsigma      |
       _cons |      1.205      0.082    14.66   0.000        1.044       1.366
------------------------------------------------------------------------------

. display exp([lnsigma]_cons)
3.3372821

其中 lnsigma 是回归的均方误差 (MSE)。现在，如果我们将其与标准 OLS 命令进行比较：

. regress mpg weight foreign

      Source |       SS           df       MS      Number of obs   =        74
-------------+----------------------------------   F(2, 71)        =     69.75
       Model |   1619.2877         2  809.643849   Prob > F        =    0.0000
    Residual |  824.171761        71   11.608053   R-squared       =    0.6627
-------------+----------------------------------   Adj R-squared   =    0.6532
       Total |  2443.45946        73  33.4720474   Root MSE        =    3.4071
------------------------------------------------------------------------------
         mpg | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
      weight |     -0.007      0.001   -10.34   0.000       -0.008      -0.005
     foreign |     -1.650      1.076    -1.53   0.130       -3.796       0.495
       _cons |     41.680      2.166    19.25   0.000       37.362      45.998
------------------------------------------------------------------------------

注意这里系数是完全相同的，但是标准误差和 MSE 项是不同的。对于 OLS，估计路径是非常不同的。regress 命令使用标准矩阵代数，而上面的 ml 选项是利用估计值的渐近特性。尽管差异非常小，但通常认为 ml 结果应该被认为比 OLS 估计值更好。

4. 参考文献

Mata 指南 -Link-
this Stata blog entry -Link-
《Maximum Likelihood Estimation with Stata, Fourth Edition》 -Link-

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh ml mate, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

Stata：最大似然估计(MLE)简易教程

1. MLE 简介

2.最大似然估计的简单例子

3.MLE 的 Stata 实现

3.1 范例 1：抛硬币例子的 MLE 估计

3.2 范例 2：Logit 模型的 MLE 估计

3.3 范例 3：Probit 模型的 MLE 估计

3.4 范例 4：使用实际数据验证

3.5 范例 5：线性回归模型的 MLE 估计

4. 参考文献

5. 相关推文

相关课程

免费公开课

最新课程-直播课

课程主页

课程主页

关于我们

热门资讯