Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
Stata 系列推文:
编译:王萃芳 (东北财经大学)
邮箱:wcfsophia@26.com
编者按:本文主要整理自 Maximum Likelihood Estimation (MLE),特此致谢!
目录
本文主要是介绍最大似然估计 (MLE) 的一些基础知识,并演示如何在 Stata 软件中进行最大似然估计。
MLE 的本质是参数估计,最合理的参数估计量应该使得从模型中抽取
让我们从抛硬币的经典例子开始。抛硬币有二种结果:正面 (
假设得到正面 (
如何利用 MLE 计算
假设概率是通过我们观察到的 5 次回合中正面的联合概率分布得出。从统计学上我们知道联合概率分布等于结果 1 的概率乘以结果 2 的概率乘以所有其他结果的概率。因此,我们的问题的似然函数可以写成:
或者是上述情况可以简化为:
其中,3 是
这是一个有未知数
接下来,对等式两边同时取对数,将似然函数转化为对数似然函数 (LL)。
对
实际上,从我们的示例中,由于我们知道
我们先生成一个变量表示投掷硬币的情况。由于我们知道对数似然 (LL) 函数等于
. 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}
是我们需要最大化的值,该值应始终写在大括号中。也可以有多个参数,每个参数都放在自己的大括号中。
. 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
。
MLE 的另一个常见示例是 Logit 或 Probit 回归。这是处理二分变量的标准工具。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
其中
由于对
这样通用 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
。我们分别指定因变量 $ML_y1
表示解释变量 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
------------------------------------------------------------------------------
除了 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) 的差异。
下面让我们使用 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
。
作者这里用 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 估计值更好。
Note:产生如下推文列表的 Stata 命令为:
lianxh ml mate, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
课程主页
课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh