Stata：数值求解极大值及 MLE 示例

发布时间：2020-12-01 阅读 4829

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

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

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

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

作者： 徐云娇 (厦门大学)，连玉君 (中山大学)
邮箱： jilyo@stu.xmu.edu.cn，arlionn@163.com

1. 简介

在求解极大值或极小值问题时，我们经常用到数值求解方法，事实上，Stata 中的 MLE 估计 (help ml) 和 NLS 估计 (help nl) 都是采用数值方法求取极值的。最简单明了的方式就是「网格搜索法」。在下图中，假设我们从 $b_{1}$ 出发来搜索极小值，希望找到最小值 $b^{*}$ 。基本步骤如下：

设定初始值 $b = b_{1}$ ，对应的函数值为 $Q (b_{1})$ ；
设定一个步长 $δ$ ，以及 $b_{2} = b_{1} + δ$ ，对应的函数值为 $Q (b_{2})$ 。显然， $Δ Q = Q (b_{2}) - Q (b_{1}) < 0$ ，表示我们的搜索方向是对的，反之，则需要将 $δ$ 设定为 $- δ$ ，即「掉头」。
重复第二步 $k$ 次，直到 $Δ Q = | Q (b_{k}) - Q (b_{k - 1}) | < ϵ$ 。其中， $ϵ$ 称为「收敛判据」，是一个预设的常数，比如 $0.001$ 。

通过上述描述可以看出，在数值求解过程中，最关键的几个参数是：初始值 $b_{1}$ ，步长 $δ$ 和收敛判据 $ϵ$ 。诸位可以想想三者有何关系？他们的参数设置或取值对于找到极值有何影响？

2. while 语句基本介绍

极值求解背后的思想是当两次搜索得到的数值之差大于收敛判据时，则继续搜索；当差值小于收敛判据时，则结束循环，输出最终搜索结果。

以上思想的实现过程需要借助 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，应注意对二者的区分。

3. 数值求解极大值

在理解了 while 语句的使用规则之后，我们便可以利用 Stata 采用数值方法求取函数极大值：

   *-示例: 采用数值法求取函数的极大值

   twoway function y = -0.2*exp(x) + ln(x^2) + 7, ///
            range(0 4) lw(*2)

假如我们想求解以上图示函数的极大值，那么在数值搜索的过程中涉及几个要点，分别是：初始值，步长，收敛判据。

求解的基本步骤是：
(1) 确定一个初始值 $x_{0}$ ，计算得到 $f (x_{0})$ ；
(2) 当 $x_{1} = x_{0} + d e l t a$ ，其中 $d e l t a$ 是步长，计算 $f (x_{1})$ ；
(3) 计算两次得到的函数值之差的绝对值， $| f (x_{1}) - f (x_{0}) |$ ，当大于收敛判据时，则继续搜索，否则就结束循环，输出最终搜索结果。

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。

需要指出的是，以上例子中只存在一个极值，但现实中还可能存在有多个极值点的情况，这个时候初始值、步长和收敛判据的选择就更加重要了。

4. 进一步讨论：全局最优解

上一节的分析对象相对比较简单，只有一个极大值。事实上，我们在进行数值求解时，是不知道目标函数的「长相」的。例如，下图中，若求取极小值，会有 B，D 和 F 几个「局部最小值」，而全局最小值在 D 处。此时，从哪里出发？步长多少？收敛判据取值多少？都会影响到我们是否能找到全局最小值。

例如，假设我们从 A 点出发，选择一个很小的步长，很可能在「好不容易」找到 B 点后便停下来「庆功会了」。然而，若我们选择一个较大的步长，则有可能直接跨过 C 点进行搜索，并很快找到 D 点。

当然，步长的选择有时也不是问题。比如，若我们从 C 点出发，即使选择一个很小的步长，也不需要花太多时间就可以到达 D 点。

情况还可以更为复杂，比如，我们要找到全局最大值，就可能面临边角解问题，如 G 点。

有些读者可能已经在感慨人生了，你的出生背景、教育环境多么像这里的「初始值」和「步长」呀，而我们又有多少人站在 E 点上开香槟，庆祝自己人生的巅峰时刻！

5. MLE估计方法的实现

回到现实中，代码还要一行一行地敲，^=^。本节介绍一下 Stata 中的 ml 命令，以便自己编写代码来实现 MLE 估计 —— 一种非常典型的数值求取极大值的方法。

5.1 MLE 简介

用一句话来概括极大似然估计，那就是：利用已知的样本结果，反推最有可能 (最大概率) 导致这样结果的参数值。

记已知的独立同分布样本集为 $D$ ：

D = {x_{1}, x_{2}, \dots, x_{N}}

与其对应的似然函数 (likelihood function) 为：

\begin{matrix} (1) & L (θ) = f (x_{1}, x_{2}, \dots, x_{N} ∣ θ) = \prod_{i = 1}^{N} f (x_{i} ∣ θ) \end{matrix}

极大似然估计 (Maximum Likelihood Estimation, MLE) 就是在求解一个最大化问题：

\begin{matrix} (2) & \hat{θ} = \arg max_{θ} L (θ) = \arg max_{θ} \prod_{i = 1}^{N} f (x_{i} ∣ θ) \end{matrix}

为了更便于计算，通常会对似然函数进行对数转换：

\begin{matrix} (3) & \hat{θ} = \arg max_{θ} \ln L (θ) = \arg max_{θ} \sum_{i = 1}^{N} \ln f (x_{i} ∣ θ) \end{matrix}

5.2 Stata 实例

利用 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
------------------------------------------------------------------------------

5.3 ml 命令中有关极大化设定的选项

A. 改变最优化算法

利用 technique(algorithm_spec) 选项可以指定似然函数最大化的方法，Stata 提供了四种不同的算法，分别是：

technique (nr): Stata's modified Newton-Raphson (NR) algorithm (默认算法)
technique (bhhh): Berndt-Hall-Hall-Hausman (BHHH) algorithm
technique (dfp): Davidon-Fletcher-Powell (DFP) algorithm
technique (bfgs): Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm

sysuse 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 算法则会出现无限循环的情况，最终无法收敛。可见：

(1) 采用不同的算法，迭代次数和结果会存在很大的差异；
(2) 采用不同的数值算法便于我们找到全局最大值(尤其是对于较难收敛的问题)；
(3) 若同时使用多种算法，Stata 默认情况下会将每种算法迭代 5 次，以便找出最佳的迭代方式。

B. 改变初始值

ml maximize：将初始值设定为 $θ = (0, 0, . . . . . ., 0)$ ；
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

不管我们如何改变初始值的设定方法，得到的估计结果都是一致的，说明此时估计方法的稳定性较高，利用数值搜索方法得到的极大值点很可能就是全局的极大值点，模型估计出的系数是可信的。

C. 改变步长、收敛判据和初始值

输入 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 默认的初始值设定 $θ = (0, 0, . . . . . ., 0)$ 。

6. 结语

除了本文介绍的方法外，在机器学习领域，还经常使用 模拟退火算法，蚁群算法，以便更为有效地找到全局最大值，日后的推文中会陆续介绍。

此外，有些读者会提出疑问：我们只需要求取目标函数的一阶偏导并令其为零不就可以得到极值了吗？对于有明确的目标函数，且能够求取偏导数的情形，这的想法是对的。但更多情况下，我们无法写出目标函数的一阶偏导，只能通过数值算法来求解。

参考资料

「极大似然估计详解，写的太好了！」
连享会推文，郭李鹏，「极大似然估计 (MLE) 及 Stata 实现」

专题	嘉宾	直播/回看视频
⭐ 最新专题		因果推断, 空间计量，寒暑假班等
⭕ 数据清洗系列	游万海	直播, 88 元，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]