温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
作者: 徐云娇 (厦门大学),连玉君 (中山大学)
邮箱: jilyo@stu.xmu.edu.cn,arlionn@163.com
目录
在求解极大值或极小值问题时,我们经常用到数值求解方法,事实上,Stata 中的 MLE 估计 (help ml
) 和 NLS 估计 (help nl
) 都是采用数值方法求取极值的。最简单明了的方式就是「网格搜索法」。在下图中,假设我们从
通过上述描述可以看出,在数值求解过程中,最关键的几个参数是:初始值
极值求解背后的思想是当两次搜索得到的数值之差大于收敛判据时,则继续搜索;当差值小于收敛判据时,则结束循环,输出最终搜索结果。
以上思想的实现过程需要借助 Stata 中的 while 循环语句,while 后面会跟一个进行逻辑判断的条件,该条件为真时,该命令会一直执行括号里面的内容,直到条件不被满足时,该循环终止。
以下这个简单的 while 循环可以帮助大家理解其工作原理:
local j = 0 //定义暂元 j,初值为 0
while `j' < 5{ //当 j 小于 5 时,执行以下命令
dis `j' //在结果窗口列示 j
local j = `j'+1 //令新的 j 等于原始值加上 1
}
*-或
local j = 0
while `j' < 5{
dis `j++' //效果等同于 dis `j' + local j = `j'+1
}
在 Stata 中运行以上代码,结果窗口便会列示 5 行,分别显示 0、1、2、3、4,可见一旦 j = 5
时,循环便会停止。
值得注意的是,除了程序中的 j++
外,还存在 ++j
这种语法,后者表示的是先把 j
加上 1 再列示 j
的当前值,所以当把以上程序中的 j++
修改为 ++j
后,结果窗口将列示 1、2、3、4、5,应注意对二者的区分。
在理解了 while 语句的使用规则之后,我们便可以利用 Stata 采用数值方法求取函数极大值:
*-示例: 采用数值法求取函数的极大值
twoway function y = -0.2*exp(x) + ln(x^2) + 7, ///
range(0 4) lw(*2)
假如我们想求解以上图示函数的极大值,那么在数值搜索的过程中涉及几个要点,分别是:初始值,步长,收敛判据。
求解的基本步骤是:
(1) 确定一个初始值,计算得到 ;
(2) 当,其中 是步长,计算 ;
(3) 计算两次得到的函数值之差的绝对值,,当大于收敛判据时,则继续搜索,否则就结束循环,输出最终搜索结果。
Stata 代码如下:
local delta = 0.05 // 步长
local x = 1 // x 的初始值
local j = 0 // 计数器:记录迭代次数
local e = 1 // y1-y0
local e0 = 0.01 // 收敛判据
while `e' > `e0'{
`trace'
local y0 = -0.2*exp(`x') + ln(`x'^2) + 7
local x = `x' + `delta'
local y1 = -0.2*exp(`x') + ln(`x'^2) + 7
local e = abs(`y1' - `y0')
dis in g "*" _c
local j = `j' + 1
}
dis "e = " `e'
dis "x = " `x' // x 的解
dis "y = " `y1' // y 的极大值
dis "j = " `j' // 迭代次数
*-图示:
twoway function y = -0.2*exp(x) + ln(x^2) + 7, ///
range(0 4) lw(thick) xline(`x') yline(`y1') ///
text(`=`y1'-0.5' `=`x'+0.8' "(`x',`y1')")
执行以上程序,经过 14 次搜索,函数差值的绝对值低至 0.0063,此时我们便近似找到了函数的极大值点 1.7,以及其对应的函数极大值 6.966467,用图形表示如下:
在整个分析的过程中,有几个关键要素,分别是:
初始值:越靠近真实极值点的初始值无疑会减少搜索的次数从而让我们更快得到极值。
步长:假如设定一个大步长,那搜索的速度就会更快;但大步长也有风险,在某些情境下可能会直接跨过极值点导致求解出不够精确的极值。
收敛判据:越严格 (即越小) 的收敛判据会使收敛条件更难被满足,可能会出现循环无法停止从而找不到极值的情况。
更换上述程序中有关这三个关键变量的初始设定,进一步测试程序的稳定性:
Q1: (更改初始值) 设定 x=2 为初始值,能否收敛?
执行程序之后未能收敛。观察函数图形易知,由于程序只能沿着坐标轴正向搜索,所以当把初始值设定为 x=2 时,两次搜索之间函数差值的绝对值只会越来越大,从而达不到收敛判据,最终未能收敛。
Q2: (更改步长) 设定 delta=0.1,结果有何变化?
重新运行程序,此时经过 8 次搜索,函数差值的绝对值低至 0.00082,函数的极大值点为 1.8,以及其对应的函数极大值 y = 6.9656438,略小于步长为 0.05 时得到的极大值。
Q3: (更改收敛判据) 设定 e0=0.0001,能否收敛?
执行程序之后未能收敛。在过小的收敛判据设定下,搜索过程中的收敛条件得不到满足,循环无法结束,所以极大值数值求解失败。
接下去进一步完善数值求解的过程,增加反向搜索功能:
*- 程序修改如下:
local h = 0.05 // 步长
local x = 1 // x 的初始值
local j = 0 // 计数器:记录迭代次数
local e = 1 // y1-y0
local e0 = `h'/10 // 收敛判据 (修改为动态数值)
while abs(`e')>`e0'{ // 修改:abs(`e')
local y0 = -0.2*exp(`x') + ln(`x'^2) + 7
local x = `x' + `h'
local y1 = -0.2*exp(`x') + ln(`x'^2) + 7
local e = `y1' - `y0' // 此前 e = abs(`y1'-`y0')
if (`e' < 0){
local h = -`h' // 新增:反向搜索
}
dis in g "*" _c
local j = `j' + 1
}
dis "e = " `e'
dis "x = " `x' // x 的解
dis "y = " `y1' // y 的极大值
dis "j = " `j' // 迭代次数
*-图示:
local x: dis %4.3f `x' // 新增:显示的更美观
local y: dis %4.3f `y1'
twoway function y = -0.2*exp(x) + ln(x^2) + 7, ///
range(0 4) lw(thick) xline(`x') yline(`y1') ///
text(`=`y'+0.5' `=`x'+0.4' "(`x', `y')")
执行以上程序,经过 15 次搜索,两次搜索的函数差值低至 0.00184,此时我们便近似找到了函数的极大值点 1.750,以及其对应的函数极大值 6.968,这比我们之前找到的极大值都略大一些,用图形表示如下:
修改后的程序不但帮助我们找到了更大的极大值,而且在其他方面优势也有很多,可以再次改变一些初始设定从而测试程序的稳定性:
(1) (更改初始值) 设定 x=3 为初始值,能否收敛?
由于新增了反向搜索的功能,所以此时执行程序之后可以收敛。经过 27 次搜索,再次锁定了极大值点 1.75,得到了与之前相同的极大值 6.968。
(2) (更改步长、收敛判据) 设定 delta=0.001,e0=0.0001,能否收敛?
可以收敛,执行以上程序,经过 691 次搜索,两次搜索的函数差值低至 0.0000986,小于收敛判据,此时我们确定了函数的极大值点 1.691,以及其对应的函数极大值 6.966。
需要指出的是,以上例子中只存在一个极值,但现实中还可能存在有多个极值点的情况,这个时候初始值、步长和收敛判据的选择就更加重要了。
上一节的分析对象相对比较简单,只有一个极大值。事实上,我们在进行数值求解时,是不知道目标函数的「长相」的。例如,下图中,若求取极小值,会有 B,D 和 F 几个「局部最小值」,而全局最小值在 D 处。此时,从哪里出发?步长多少?收敛判据取值多少?都会影响到我们是否能找到全局最小值。
例如,假设我们从 A 点出发,选择一个很小的步长,很可能在「好不容易」找到 B 点后便停下来「庆功会了」。然而,若我们选择一个较大的步长,则有可能直接跨过 C 点进行搜索,并很快找到 D 点。
当然,步长的选择有时也不是问题。比如,若我们从 C 点出发,即使选择一个很小的步长,也不需要花太多时间就可以到达 D 点。
情况还可以更为复杂,比如,我们要找到全局最大值,就可能面临边角解问题,如 G 点。
有些读者可能已经在感慨人生了,你的出生背景、教育环境多么像这里的「初始值」和「步长」呀,而我们又有多少人站在 E 点上开香槟,庆祝自己人生的巅峰时刻!
回到现实中,代码还要一行一行地敲,^=^。本节介绍一下 Stata 中的 ml
命令,以便自己编写代码来实现 MLE 估计 —— 一种非常典型的数值求取极大值的方法。
用一句话来概括极大似然估计,那就是:利用已知的样本结果,反推最有可能 (最大概率) 导致这样结果的参数值。
记已知的独立同分布样本集为
与其对应的似然函数 (likelihood function) 为:
极大似然估计 (Maximum Likelihood Estimation, MLE) 就是在求解一个最大化问题:
为了更便于计算,通常会对似然函数进行对数转换:
利用 Stata 实现 MLE 估计更加详细的介绍参见「极大似然估计 (MLE) 及 Stata 实现」,这里对 help ml
中有关最优化求解的选项进行介绍:
以简单的正态分布函数为例,并结合汽车价格数据,一个标准的 MLE 求解过程有两步:
第一步:将似然函数的基本设定编写为 ado 文件,并在 dofile 中选中相应代码,右键选择 "Tools" → "Execute selection quietly (run)";
cap program drop mymean_lf
program define mymean_lf
version 9.1 //设定 Stata 9.0版本,以便更高版本能够兼容
args lnf mu sigma //输入项:似然函数,参数1,参数2,......
quietly replace `lnf' = ln(normalden($ML_y1, `mu', `sigma'))
end
/*$ML_y1 默认写法,表示被解释变量*/
第二步:调入汽车价格的数据,利用 ml
命令进行极大似然估计。
sysuse auto, clear
ml model lf mymean_lf (price = mpg wei len) (sigma: )
/* lf 是最常用的架构,通常对样本满足独立同分布的条件时均可使用;
假定汽车价格 (price) 服从均值为 mu,方差为 sigma 的正态分布,
并将 mu 设定为汽车每公里耗油量 (mpg)、重量 (wei)、长度 (len) 的线性函数*/
ml search
ml maximize
执行以上代码,便可得到估计结果:
Number of obs = 74
Wald chi2(3) = 41.15
Log likelihood = -679.35161 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
price | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
eq1 |
mpg | -86.789 81.643 -1.06 0.288 -246.807 73.228
weight | 4.365 1.135 3.84 0.000 2.139 6.590
length | -104.868 38.633 -2.71 0.007 -180.587 -29.149
_cons | 14542.430 5729.204 2.54 0.011 3313.396 25771.464
-------------+----------------------------------------------------------------
sigma |
_cons | 2348.394 193.036 12.17 0.000 1970.050 2726.738
------------------------------------------------------------------------------
利用 technique(algorithm_spec)
选项可以指定似然函数最大化的方法,Stata 提供了四种不同的算法,分别是:
technique (nr)
: Stata's modified Newton-Raphson (NR) algorithm (默认算法)technique (bhhh)
: Berndt-Hall-Hall-Hausman (BHHH) algorithmtechnique (dfp)
: Davidon-Fletcher-Powell (DFP) algorithmtechnique (bfgs)
: Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithmsysuse auto, clear
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml maximize
est store nr
ml model lf mymean_lf (price = mpg wei len) (sigma:), technique(bfgs)
ml maximize
est store bfgs
ml model lf mymean_lf (price = mpg wei len) (sigma:), technique(nr bfgs)
ml maximize
est store nr_bfgs
local model nr bfgs nr_bfgs
esttab `model', mtitle(`model')
上面三种算法的估计结果完全一致,故不再展示。
但是在以上实例中,采用 bhhh
算法会在运行过程中提示出错,从而得不到结果;采用 dfp
算法则会出现无限循环的情况,最终无法收敛。可见:
ml maximize
:将初始值设定为 ml search
:Stata 自身提供了多种快速搜索初始值的方法;ml init
:可以根据理论或前期估计结果设定初始值。sysuse auto, clear
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml max //maximize 简写为 max
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml search
ml max
ml model lf mymean_lf (price = mpg wei len) (sigma:)
ml init /sigma = 3 //设定初始值
ml max
不管我们如何改变初始值的设定方法,得到的估计结果都是一致的,说明此时估计方法的稳定性较高,利用数值搜索方法得到的极大值点很可能就是全局的极大值点,模型估计出的系数是可信的。
输入 help ml##noninteractive_maxopts
或 help maximize
可以查看这部分选项的设定:
maximize_options: difficult, technique(algorithm_spec), iterate(#), [no]log, trace, gradient, showstep, hessian, showtolerance, tolerance(#), ltolerance(#), nrtolerance(#), nonrtolerance, and from(init_specs);
see [R] maximize. These options are seldom used.
虽然在日常分析中很少设定这些选项,但在面临一些复杂模型时,通过这些选项的设定可以让我们「暂时」获得一组参数估计值,以便判断程序设定是否正确,并酌情做出调整。这里,仅对最为关键的几个选项进行说明。
iterate(#)
:用于设定迭代次数,默认值为 16000。调试程序时,可以设定一个较小的数值,以便确认程序是否能够正常运行,同时,可以配合 timer
命令初步估算模型运行的耗时,这在配合使用 Bootstrap 获取标准误或置信区间时非常有用。
ltolerance(#)
:参数向量的收敛判据,当两次估计得到的参数向量的差值小于收敛判据时,便视为收敛。若采用 MLE 进行估计,即 ml
命令,#
的默认值为 1e-4,即 0.0001;对于其他极值问题,如采用 nl
执行 NLS 估计或采用 gmm
执行 GMM 估计,#
的默认值为 1e-6。
ltolerance(#)
:对数似然函数的收敛判据,默认值设定标准同上。显然,当模型长时间无法收敛时,我们可以设定一个大一点的 #
值,代价是估计出的参数不够准确,或无法找到全局最优解。
from(init_specs)
:人为指定参数的初始值。比如,我们可以先采用 OLS 或一些有偏但易于执行的估计量得到参数的初始值,这通常要优于 Stata 默认的初始值设定
除了本文介绍的方法外,在机器学习领域,还经常使用 模拟退火算法,蚁群算法,以便更为有效地找到全局最大值,日后的推文中会陆续介绍。
此外,有些读者会提出疑问:我们只需要求取目标函数的一阶偏导并令其为零不就可以得到极值了吗?对于有明确的目标函数,且能够求取偏导数的情形,这的想法是对的。但更多情况下,我们无法写出目标函数的一阶偏导,只能通过数值算法来求解。
连享会推文,郭李鹏,「极大似然估计 (MLE) 及 Stata 实现」
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟,课程主页 Stata 33 讲 - 连玉君, 每讲 15 分钟. Stata 小白的取经之路 - 龙志能,时长:2 小时,课程主页 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 因果推断, 空间计量,寒暑假班等 | |
⭕ 数据清洗系列 | 游万海 | 直播, 88 元,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等lianxh
命令发布了: 在 Stata 命令窗口中输入 ssc install lianxh
即可安装,随时搜索连享会推文、Stata 资源,详情:help lianxh
。
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD