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 系列推文:

PDF下载 - 推文合集

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

编者按:本文主要整理自 Maximum Likelihood Estimation (MLE),特此致谢!


目录


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

1. MLE 简介

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

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

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

假设得到正面 (H) 的概率等于某个未知值 p 或者 Pr(H)=p。那么,得到反面的概率,即没有得到正面的概率 (H 或者 1H) 为 1p,公式为 Pr(T)=Pr(1H)=1p

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

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

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

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

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

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

对 p 求导,并求得 p

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

3.MLE 的 Stata 实现

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

我们先生成一个变量表示投掷硬币的情况。由于我们知道对数似然 (LL) 函数等于 LL(p)=3×log(p)+2×log(1p),因此可以绘制下图:

. 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 的可能性来自逻辑分布,它采用以下函数形式:

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

或者 y=0 的概率:

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

在这两个似然分布中,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。对于上述观察结果,似然函数可以推导为:

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

这样通用 Logit LL 可以写成:

为了在 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 probitlogit 命令利用了 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小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh