🔥 热搜 🔥
上海习近平新疆鄂州父女瓜乌鲁木齐疫情H工口小学生赛高习明泽芊川一笑图包印尼排华
分类
社会娱乐国际人权科技经济其它
🔥 热搜 🔥
上海习近平新疆鄂州父女瓜乌鲁木齐疫情H工口小学生赛高习明泽芊川一笑图包印尼排华
分类
社会娱乐国际人权科技经济其它
查看原文
其他

Stata：最大似然估计入门教程(MLE)-ml

连享会 连享会 2023-10-24
收录于合集 #Stata命令 335个

 

👇 连享会 · 推文导航 | www.lianxh.cn

 

连享会 · 2022 空间计量专题

 

 

作者：曹昊煜 (兰州大学)
邮箱：caohy19@lzu.edu.cn

 

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

 

目录

  • 1. 最大似然估计与条件最大似然估计

    • 1.1 M 估计量

    • 1.2 最大似然估计

    • 1.3 条件最大似然估计

  • 2. 最大似然估计两种方法

    • 2.1 最大似然估计：极简示例

    • 2.2 最大似然估计：网格搜索

  • 3. Stata 范例

    • 3.1 范例一：估计正态分布的参数

    • 3.2 范例二：线性模型

    • 3.3 范例三：定性反应模型

    • 3.4 范例四：断尾回归

    • 3.5 范例五：截取回归

  • 4. 结束语

  • 5. 相关推文


 

 

1. 最大似然估计与条件最大似然估计

1.1 M 估计量

计量经济学中一类重要的估计量称为 “极值估计量”，该估计量通过最大化某个目标函数获得，即：

极值估计量又可以分为 M 估计量和 GMM 估计量。最大似然估计属于 M 估计量的一种，其目标函数的一般形式为样本均值：

1.2 最大似然估计

最大似然估计主要是将目标函数中的 替换为密度函数的对数形式。令 为独立同分布的序列，且概率密度为 ，则数据 的联合密度为：

当使用参数的假象值替代真实参数时，联合密度称为似然函数。真实参数 的 ML 估计是选择相应的 ，使得似然函数最大化。由于对数变换是单调的，最大化似然函数等价于对数似然函数最大化：

因此 的 ML 估计量为 M 估计量：

在一定的正则假定下，ML 估计量为一致渐进正态估计量。

1.3 条件最大似然估计

在一般的研究问题中，数据 可以区分为解释变量和被解释变量，而参数估计的主要目标是考察解释变量 对被解释变量 条件分布的影响。根据条件概率密度公式，数据的联合密度可以表示为：

对等式两边取对数，并使用假想的参数替代真实参数得到对数似然函数：

若条件概率密度和解释变量的概率密度之间不存在函数相关，忽略等式右侧的第二项得到简化的对数自然函数。对应的 M 估计量的形式为：

因此该方法称为条件极大似然估计。在一般的应用中，条件极大似然估计更为普遍，其原因在于我们对联合分布中的参数建模后，被解释变量的分布会依赖于解释变量的变化。

 

2. 最大似然估计两种方法

2.1 最大似然估计：极简示例

首先使用一个简单的例子来说明 MLE，该示例摘自 B 站 Up 主「小崔说数」的讲解内容。假设城市居民只有两种出行方式，分别是步行 () 或开车 ()。二者发生的概率分别为 和 。随机询问五个人，结果为 ，从而可以得到：


10011






在该模型中，相当于估计参数 。由于样本是随机抽取的，因此联合概率密度等于所有样本概率之积，定义似然函数：

对等式两边取对数，则有：

求取一阶条件：

估计结果 与 的样本均值相同，因为在二值观测中，期望与取值为 1 的概率相同。通过上面这个简单的例子，可以看出 MLE 的基本步骤：

  • 第一步：确定 的概率分布；
  • 第二步：写出样本中所有观测值的联合分布 ；
  • 第三步：极大化 从而得到 的估计。

2.2 最大似然估计：网格搜索

求解最大似然估计 (极大似然估计) 的结果可以使用一阶条件，或称为似然方程，但在某些情形下，并不存在参数的显式解，因此需要使用数值方法求解参数。一种简单的，便于考察 MLE 整个过程的方式是网格搜索，当真实参数范围已知时，可以通过搜索的方式估计参数。其过程如下：

  • 以某个特定的步长在参数的取值范围内选择不同的参数值；
  • 计算特定参数取值下的似然函数值；
  • 选取最大化似然函数取值对应的参数作为真实参数的估计。

最大似然估计的一个简单例子是估计正态分布的均值和标准差。由正态分布的概率密度函数，平均对数似然函数为：

假定数据生成过程服从均值为 1，标准差为 0.8 的正态分布，则网格搜索的基本思路是同时选择不同的 和 ，并计算对数似然值并比较其大小。过程如下：

. clear
. set obs 10000
. set seed 1010110 // 设定种子值
. gen y = rnormal(1,0.8) in 1/1000 // 生成随机样本 y ~ N(1,0.8)
. gen mu = .
. gen sigma = .
. gen logL = . // 生成用于保存均值、标准差和对数似然值的
. local n = 1

. * 开始网格搜索
. forvalue mu = 0.01(0.01)1{
2. forvalues sigma = 0.01(0.01)1{
3. gen term1 = (y-`mu')^2/(2*`sigma'^2)
4. egen term2 = sum(term1)
5. qui sum term2
6. local last = r(max)
7. local N = 1000
8. replace logL = -0.5*`N'*log(2*_pi) - 0.5*`N'*log(`sigma'^2) - `last' in `n' // 保存对数似然函数值
9. replace mu = `mu' in `n'
10. replace sigma = `sigma' in `n'
11. drop term1 term2
12. local n = `n'+1
13. }
14. }

. * 搜索结果
. gsort - logL
. list mu sigma logL in 1/10

+-------------------------+
| mu sigma logL |
|-------------------------|
1. | .99 .81 -1208.978 |
2. | .99 .82 -1209.109 |
3. | .99 .8 -1209.152 |
4. | .99 .83 -1209.527 |
5. | .99 .79 -1209.652 |
|-------------------------|
6. | .98 .81 -1209.714 |
7. | .98 .82 -1209.827 |
8. | .98 .8 -1209.907 |
9. | .99 .84 -1210.215 |
10. | .98 .83 -1210.228 |
+-------------------------+

可以看出，均值和标准差的搜索结果与真实的参数非常接近。当然，在具体的 MLE 中，我们并不知道真实参数的取值范围，因此需要使用牛顿-拉夫森等其他数值计算方法，在以下的例子中 Stata 会自动实现数值计算。

 

3. Stata 范例

在 Stata 中，可以使用 ml 命令灵活地实现各类最大似然估计。关于 ml 命令的使用和设定方式可以参考帮助文档 help ml 或连享会推文「Stata：数值求解极大值及 MLE 示例」。本文将通过五个简单的范例，介绍如何在 Stata 中实现基本的 MLE 方法。

3.1 范例一：估计正态分布的参数

第一个例子仍然是对正态分布中均值和标准差的估计。此时我们不再使用网格搜索，而是使用其他数值方法，概率密度函数仍然服从正态分布：

在 Stata 中使用 ml 命令，只需要定义特定观测值下的概率密度，而不是整体的对数似然函数，例如观测值 的概率密度是均值为 ，标准差为 的正态分布概率密度函数。

. clear
. set obs 1000
. set seed 1010110
. gen y = rnormal(1,0.8) // 生成随机样本 y ~ N(1,0.8)
. cap program drop mynormal_lf
. program define mynormal_lf // 定义对数似然函数
1. version 16
2. args lnf mu sigma
3. qui replace `lnf' = ln(normalden($ML_y1,`mu',`sigma'))
4. end
. ml model lf mynormal_lf (y = ) (sigma: ) // 模型形式
. ml maximize // 参数求解

Iteration 0: log likelihood = -1247.1045
Iteration 1: log likelihood = -1207.7554
Iteration 2: log likelihood = -1207.5499
Iteration 3: log likelihood = -1207.5496
Iteration 4: log likelihood = -1207.5496
Number of obs = 1,000
Wald chi2(0) = .
Log likelihood = -1207.5496 Prob > chi2 = .
------------------------------------------------------------------------------
y | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
_cons | 1.033 0.026 40.37 0.000 0.983 1.083
-------------+----------------------------------------------------------------
sigma |
_cons | 0.809 0.018 44.72 0.000 0.774 0.845
------------------------------------------------------------------------------

关于似然函数的定义需要注意以下几点：

  • cap program drop：用于删去可能存在的重复程序；
  • progname_lf：是程序命名的一般形式，即自己确定的程序名称加上下划线和 lf 后缀；
  • version 16：声明 Stata 版本；
  • args：后面紧跟对数似然函数值 lnf 和参数名称；
  • $ML_y1：代表被解释变量。

估计结果中显示了 MLE 的迭代过程，最大化的对数似然值为 -1207.5496，估计出的均值为 1.033，标准差为 0.809，与设定的数据生成过程基本相同。

3.2 范例二：线性模型

对于均值和标准差的估计属于无条件最大似然估计，因为我们并没有对分布中的参数进行建模，也相当于对截距项回归。接下来考察普通的线性模型：

假设随机干扰项 。从而 ，即对正态分布的均值进行了线性建模。相应的条件对数似然函数为：

使用 Stata 进行线性模型的 ML 估计，并与 regress 命令做对比：

. sysuse auto.dta, clear
. cap program drop mylinear_lf
. program define mylinear_lf
1. version 16
2. args lnf xb sigma
3. qui replace `lnf' = ln(normalden($ML_y1,`xb',`sigma'))
4. end
. ml model lf mylinear_lf (price = foreign trunk weight length) (sigma: )
. ml maximize, nolog

Number of obs = 74
Wald chi2(4) = 90.07
Log likelihood = -666.25274 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
price | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
foreign | 3580.051 624.022 5.74 0.000 2356.990 4803.111
trunk | -10.282 78.824 -0.13 0.896 -164.774 144.210
weight | 5.769 0.934 6.18 0.000 3.938 7.600
length | -89.653 34.534 -2.60 0.009 -157.339 -21.968
_cons | 4672.931 3852.845 1.21 0.225 -2878.506 12224.368
-------------+----------------------------------------------------------------
sigma |
_cons | 1967.417 161.721 12.17 0.000 1650.451 2284.384
------------------------------------------------------------------------------

. reg price foreign trunk weight length

Source | SS df MS Number of obs = 74
-------------+---------------------------------- F(4, 69) = 21.00
Model | 348631330 4 87157832.5 Prob > F = 0.0000
Residual | 286434066 69 4151218.35 R-squared = 0.5490
-------------+---------------------------------- Adj R-squared = 0.5228
Total | 635065396 73 8699525.97 Root MSE = 2037.5
------------------------------------------------------------------------------
price | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
foreign | 3580.051 646.236 5.54 0.000 2290.845 4869.256
trunk | -10.282 81.630 -0.13 0.900 -173.129 152.565
weight | 5.769 0.967 5.96 0.000 3.839 7.699
length | -89.653 35.763 -2.51 0.015 -160.999 -18.307
_cons | 4672.931 3989.999 1.17 0.246 -3286.899 12632.762
------------------------------------------------------------------------------

可见 MLE 与 OLS 的估计结果完全相同，仅仅在标准误上存在细微的差异，如果考察二者的一阶条件，则可以发现二者面临的优化问题是相同的。

3.3 范例三：定性反应模型

接下来考察几个非线性模型的例子。当被解释变量取值为 0-1 时，需要构建定性反应模型，其基本的概率分布为：

对参数 进行建模，由于其代表的是 取 1 的概率，因此其函数形式可以定义为正态分布或 Logistic 分布。当假定 服从正态分布时，该模型称为 Probit 模型：

当假定 服从 Logistic 分布时，该模型称为 Logit 模型：

其中，

在 Stata 中，normal 和 invlogit 分别表示正态分布累积分布函数和拟 Logit 变换，probit 和 logit 命令分别用于估计上述模型。

. * probit 模型
. sysuse auto.dta, clear
. cap program drop myprobit_lf
. program define myprobit_lf
1. version 16
2. args lnf xb
3. qui replace `lnf' = ln(normal(`xb')) if $ML_y1 == 1
4. qui replace `lnf' = ln(1 - normal(`xb')) if $ML_y1 == 0
5. end
. ml model lf myprobit_lf (foreign = price trunk weight length)
. ml maximize, nolog

Number of obs = 74
Wald chi2(4) = 15.11
Log likelihood = -17.72839 Prob > chi2 = 0.0045
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.32 0.001 0.000 0.001
trunk | 0.010 0.083 0.12 0.901 -0.153 0.174
weight | -0.004 0.002 -2.71 0.007 -0.007 -0.001
length | 0.033 0.045 0.74 0.460 -0.055 0.121
_cons | 1.201 5.070 0.24 0.813 -8.735 11.138
------------------------------------------------------------------------------

. probit foreign price trunk weight length, nolog

Probit regression Number of obs = 74
LR chi2(4) = 54.61
Prob > chi2 = 0.0000
Log likelihood = -17.72839 Pseudo R2 = 0.6063
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.32 0.001 0.000 0.001
trunk | 0.010 0.083 0.12 0.901 -0.153 0.174
weight | -0.004 0.002 -2.71 0.007 -0.007 -0.001
length | 0.033 0.045 0.74 0.460 -0.055 0.121
_cons | 1.201 5.070 0.24 0.813 -8.735 11.138
------------------------------------------------------------------------------. * logit 模型
. cap program drop mylogit_lf
. program define mylogit_lf
1. version 16
2. args lnf xb
3. qui replace `lnf' = ln(invlogit(`xb')) if $ML_y1 == 1
4. qui replace `lnf' = ln(1 - invlogit(`xb')) if $ML_y1 == 0
5.
. * qui replace `lnf' = ln(exp(`xb')/(1+exp(`xb'))) if $ML_y1 == 1
. * qui replace `lnf' = ln(1 - exp(`xb')/(1+exp(`xb'))) if $ML_y1 == 0
. end
. ml model lf mylogit_lf (foreign = price trunk weight length)
. ml maximize, nolog

Number of obs = 74
Wald chi2(4) = 12.61
Log likelihood = -17.786494 Prob > chi2 = 0.0134
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.14 0.002 0.000 0.002
trunk | 0.009 0.142 0.06 0.949 -0.269 0.287
weight | -0.007 0.003 -2.55 0.011 -0.013 -0.002
length | 0.050 0.081 0.61 0.539 -0.109 0.209
_cons | 3.388 9.308 0.36 0.716 -14.855 21.631
------------------------------------------------------------------------------

. logit foreign price trunk weight length, nolog

Logistic regression Number of obs = 74
LR chi2(4) = 54.49
Prob > chi2 = 0.0000
Log likelihood = -17.786494 Pseudo R2 = 0.6050
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.14 0.002 0.000 0.002
trunk | 0.009 0.142 0.06 0.949 -0.269 0.287
weight | -0.007 0.003 -2.55 0.011 -0.013 -0.002
length | 0.050 0.081 0.61 0.539 -0.109 0.209
_cons | 3.388 9.308 0.36 0.716 -14.855 21.631
------------------------------------------------------------------------------

从估计结果中，无论是 Probit 还是 Logit 模型，直接编写的最大似然估计结果与官方命令提供的结果基本没有差异。

3.4 范例四：断尾回归

有时候，我们获取的样本并非从总体中随机抽样得到，而是从总体的子样本中获取的，例如规模以上的企业调查。假设随机变量 的分布为 ，则 的断尾分布密度函数为：

假设无断尾 的分布为 ，则在 处断尾的分布为：

从而对数似然函数为：

在 Stata 中，truncreg 命令可以用于估计断尾回归，因此可以使用该命令验证自己编写的断尾回归代码。

. lxhuse laborsub.dta, clear
. keep if whrs > 0 // 假设在 0 处断尾
. cap program drop mytrunc_lf
. program define mytrunc_lf
1. version 16
2. args lnf xb sigma
3. qui replace `lnf' = ln(((1/`sigma')*normalden(($ML_y1-`xb')/`sigma'))/(1-normal((0-`xb')/`sigma')))
4. end
. ml model lf mytrunc_lf (whrs = kl6 k618 wa we) (sigma: )
. ml maximize, nolog

Number of obs = 150
Wald chi2(4) = 10.05
Log likelihood = -1200.9157 Prob > chi2 = 0.0395
------------------------------------------------------------------------------
whrs | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
kl6 | -803.004 321.361 -2.50 0.012 -1432.860 -173.148
k618 | -172.875 88.729 -1.95 0.051 -346.780 1.031
wa | -8.821 14.368 -0.61 0.539 -36.983 19.341
we | 16.529 46.504 0.36 0.722 -74.617 107.674
_cons | 1586.260 912.354 1.74 0.082 -201.922 3374.441
-------------+----------------------------------------------------------------
sigma |
_cons | 983.726 94.443 10.42 0.000 798.621 1168.830
------------------------------------------------------------------------------

. truncreg whrs kl6 k618 wa we, ll(0) nolog

Truncated regression
Limit: Lower = 0 Number of obs = 150
Upper = +inf Wald chi2(4) = 10.05
Log likelihood = -1200.9157 Prob > chi2 = 0.0395
------------------------------------------------------------------------------
whrs | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
kl6 | -803.004 321.361 -2.50 0.012 -1432.861 -173.147
k618 | -172.875 88.729 -1.95 0.051 -346.781 1.031
wa | -8.821 14.368 -0.61 0.539 -36.983 19.341
we | 16.529 46.504 0.36 0.722 -74.617 107.674
_cons | 1586.260 912.355 1.74 0.082 -201.923 3374.442
-------------+----------------------------------------------------------------
/sigma | 983.726 94.443 10.42 0.000 798.621 1168.831
------------------------------------------------------------------------------

3.5 范例五：截取回归

截取由 Tobin (1958) 首次引入经济学，与断尾分布不同，此时样本仍然是从总体中随机抽取的，而当被解释变量的取值小于某个特定的值时，全部被归并到该取值上。截取模型又称为归并模型或 Tobit 模型。其形式为潜变量模型：

可以将 的分布视为两个部分， 的部分是正常的观测值，其分布与潜变量 相同，而 的概率与 的概率相同。因此 Tobit 模型的概率密度由两个部分组成：

相应地，条件对数似然函数为：

使用 tobit 命令来验证自己编写的 MLE 程序：

. lxhuse womenwk.dta, clear
. cap program drop mytobit_lf
. program define mytobit_lf
1. version 16
2. args lnf xb sigma
3. qui replace `lnf' = ln((1/`sigma')*normalden(($ML_y1-`xb')/`sigma')) if $ML_y1 != 0
4. qui replace `lnf' = ln(normal((0-`xb')/`sigma')) if $ML_y1 == 0
5. end
. ml model lf mytobit_lf (lwf = age married children education) (sigma: )
. ml maximize, nolog

Number of obs = 2,000
Wald chi2(4) = 472.59
Log likelihood = -3349.9685 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
lwf | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
age | 0.052 0.006 9.08 0.000 0.041 0.063
married | 0.484 0.104 4.68 0.000 0.281 0.687
children | 0.486 0.032 15.33 0.000 0.424 0.548
education | 0.115 0.015 7.62 0.000 0.085 0.145
_cons | -2.808 0.263 -10.67 0.000 -3.324 -2.292
-------------+----------------------------------------------------------------
sigma |
_cons | 1.873 0.040 46.80 0.000 1.794 1.951
------------------------------------------------------------------------------

. tobit lwf age married children education, ll(0) nolog

Tobit regression Number of obs = 2,000
Uncensored = 1,343
Limits: Lower = 0 Left-censored = 657
Upper = +inf Right-censored = 0
LR chi2(4) = 461.85
Prob > chi2 = 0.0000
Log likelihood = -3349.9685 Pseudo R2 = 0.0645
------------------------------------------------------------------------------
lwf | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
age | 0.052 0.006 9.08 0.000 0.041 0.063
married | 0.484 0.104 4.68 0.000 0.281 0.687
children | 0.486 0.032 15.33 0.000 0.424 0.548
education | 0.115 0.015 7.62 0.000 0.085 0.145
_cons | -2.808 0.263 -10.67 0.000 -3.324 -2.291
-------------+----------------------------------------------------------------
var(e.lwf)| 3.507 0.150 3.225 3.814
------------------------------------------------------------------------------

MLE 的估计结果与 Tobit 模型的估计结果完全相同。

 

4. 结束语

本文主要介绍了最大似然估计和条件最大似然估计的理论过程，并使用网格搜索对 MLE 中的数值计算进行了简单的示例，最后使用五个范例说明了在 Stata 中如何使用 ml 命令进行最大似然估计。

本文的范例主要涉及线性和非线性的单方程模型。实际上多方程模型同样可以使用 MLE，具体理论包括似无相关回归、完全信息极大似然估计和有限信息极大似然估计等。此外，ml 命令还可以进行更加复杂的设定，我们将在之后的推文中介绍。

 

5. 相关推文

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

  • 专题：Stata命令
    • Stata新命令：面板-LogitFE-ProbitFE
  • 专题：Stata程序
    • Stata：数值求解极大值及 MLE 示例
    • 极大似然估计 (MLE) 及 Stata 实现
  • 专题：回归分析
    • 多层级Tobit模型及Stata应用
  • 专题：内生性-因果推断
    • Stata: 一文读懂 Tobit 模型
  • 专题：交乘项-调节
    • Logit-Probit中的交乘项及边际效应图示
  • 专题：Probit-Logit
    • Logit-Probit：非线性模型中交互项的边际效应解读
    • 详解 Logit/Probit 模型中的 completely determined 问题
    • 二元选择模型：Probit 还是 Logit？
    • Stata：二元Probit模型
    • 动态 Probit 模型及 Stata 实现

 

课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp

 

New！ Stata 搜索神器：lianxh 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

 

🍏 关于我们

  • 连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。
  • 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。


您可能也对以下帖子感兴趣

文章有问题？点此查看未经处理的缓存