最大似然估计（Maximum Likelihood Estimation; MLE）

class: center, middle, inverse, title-slide

.title[
# 最大似然估计（Maximum Likelihood Estimation; MLE）
]
.author[
### 文旷宇
]
.institute[
### 华中科技大学
]
.date[
### 2023
]

---

class: left, middle
layout: true

<div>
<style type="text/css">.xaringan-extra-logo {
width: 55px;
height: 64px;
z-index: 0;
background-image: url(../images/hust-logo.png);
background-size: contain;
background-repeat: no-repeat;
position: absolute;
top:1em;right:1em;
}
</style>
<script>(function () {
  let tries = 0
  function addLogo () {
    if (typeof slideshow === 'undefined') {
      tries += 1
      if (tries < 10) {
        setTimeout(addLogo, 100)
      }
    } else {
      document.querySelectorAll('.remark-slide-content:not(.title-slide):not(.inverse):not(.hide_logo)')
        .forEach(function (slide) {
          const logo = document.createElement('div')
          logo.classList = 'xaringan-extra-logo'
          logo.href = null
          slide.appendChild(logo)
        })
    }
  }
  document.addEventListener('DOMContentLoaded', addLogo)
})()</script>
</div><script>document.addEventListener('DOMContentLoaded',function(){new xeBanner(JSON.parse('{"left":"","right":"","exclude":["title-slide"],"position":"top"}'))})</script>
<script>document.addEventListener('DOMContentLoaded',function(){new xeBanner(JSON.parse('{"left":"文旷宇/华中科技大学","exclude":["title-slide"],"position":"bottom"}'))})</script>

---
### 主要内容

- 简单线型回归的最大似然估计

- 最大似然估计的一般框架

- 最大似然估计量的计算

- 最大似然估计量的理论性质

- 基于最大似然估计的统计推断：似然比检验

- 拟最大似然估计的理论性质

- 模型选择：交叉验证与赤池信息准则

---
### 以简单线性回归为例

- 独立同分布样本 `$(Y_i,X_i): i=1,\cdots,n$`，考虑简单回归 
$$
Y_i = \beta_0 +\beta_1 X_i + e_i，
$$
其中残差项 `$e_i$` 服从正态分布，即 `$e_i \sim N(0,\sigma^2)$`。

- 未知参数 `$(\beta_0,\beta_1,\sigma^2)$`；如何进行最大似然估计？

- 关键步骤：推导出 `$(Y_i,X_i)$` 的联合概率密度函数，进而写出似然函数。

- 一些记号：

- `$(Y_i,X_i)$` 的联合概率密度函数： `$f_{XY}$`；
  - `$X_i$` 的边缘概率密度函数： `$f_X$`；
  - `$Y_i|X_i$` 的条件概率密度函数： `$f_{Y|X}$`。

---

### 似然函数

- 条件概率密度公式： `$f_{XY}(X_i,Y_i) = f_{Y|X}(Y_i|X_i) f_X(X_i)$`。

- 根据假设， `$Y_i|X_i$` 服从正态分布，即 `$Y_i|X_i \sim N(\beta_0+\beta_1 X_i,\sigma^2)$`。那么，
$$
f_{Y|X}(Y_i|X_i) = \frac{1}{\sqrt{2\pi}\sigma}\exp \left(-\frac{(Y_i-\beta_0-\beta_1 X_i)^2}{2\sigma^2} \right)。
$$

- 由于样本独立同分布，对数似然函数可表示为 `$$\begin{aligned} & l(\beta_0,\beta_1,\sigma^2) = \sum_{i=1}^{n}\log f_{XY}(X_i,Y_i) \\ &= \sum_{i=1}^{n}\log f_{Y|X}(Y_i|X_i) + \sum_{i=1}^{n}\log f_X(X_i) \\ &= -\frac{n}{2}\log(2\pi)-\frac{n}{2}\log(\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(Y_i-\beta_0-\beta_1 X_i)^2  + \sum_{i=1}^{n}\log f_X(X_i)。 \end{aligned}$$`
注意：最后一项，即边缘密度 `$f_X$` 与未知参数无关。

---

### 最大似然估计量
- 通过最大化对数似然函数 `$l(\beta_0,\beta_1,\sigma^2)$` 得出最大似然估计量，记作 `$(\hat \beta_{0,MLE},\hat \beta_{1,MLE},\hat \sigma_{MLE}^2)$`。
- 我们有 `$$(\hat \beta_{0,MLE},\hat \beta_{1,MLE},\hat \sigma_{MLE}^2) =  \arg\max_{\beta_0,\beta_1 \in \mathbb{R},\sigma^2 \in \mathbb{R}_{+}} l(\beta_0,\beta_1,\sigma^2)。$$`

- 一阶偏导数 `$$\begin{aligned} \frac{\partial l}{\partial \beta_0} &= \frac{1}{\sigma^2}\sum_{i=1}^{n}(Y_i-\beta_0-\beta_1 X_i)  \\ \frac{\partial l}{\partial \beta_1} &= \frac{1}{\sigma^2}\sum_{i=1}^{n}X_i(Y_i-\beta_0-\beta_1 X_i) \\ \frac{\partial l}{\partial \sigma^2} &= -\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2}\sum_{i=1}^{n}(Y_i-\beta_0-\beta_1 X_i)^2 \end{aligned}$$`

---

### 最大似然估计量

- 一阶条件：令 `$\frac{\partial l}{\partial \beta_0}=0$` 以及 `$\frac{\partial l}{\partial \beta_1}=0$` 得出 `$$\begin{aligned} & \bar{Y}=\hat \beta_{0,MLE} + \hat \beta_{1,MLE} \bar{X} \\ & \hat \beta_{1,MLE} = \frac{\sum_{i=1}^n (Y_i-\bar Y)(X_i-\bar X)}{\sum_{i=1}^n(X_i-\bar X)^2}， \end{aligned}$$` 其中 `$\bar X = n^{-1}\sum_{i=1}^n X_i$` 和 `$\bar Y = n^{-1}\sum_{i=1}^n Y_i$`。

- 一阶条件：令 `$\frac{\partial l}{\partial \sigma^2}=0$`，得出 `$$\hat \sigma_{MLE}^2 = \frac{1}{n} \sum_{i=1}^{n}(Y_i-\hat \beta_{0,MLE}-\hat \beta_{1,MLE} X_i)^2。$$`

---

### OLS还是MLE？

- 普通最小二乘估计量（OLS）无需假设残差项的分布，而最大似然估计量（MLE）需要假设残差项的分布。

- 线性回归中，回归系数的OLS估计量与假设残差项为正态分布的MLE估计量相同，此时OLS估计量可以视作MLE估计框架下的一个特例。

- 当残差项不是正态分布时：
  - 如果残差项的分布类型有充分的理论支持，那么MLE吧，接下来你会发现MLE相对于OLS的优势；
  - 如果你无法确信残差项的分布的类型，那么OLS吧，做“最保守”的选择，并且OLS的计算比MLE简单太多。

---

### 从简单线性回归例子过渡到一般性最大似然框架

|         | 简单线性回归  | 一般性最大似然框架|
| - |:-:| :-:|
| 样本      | `$(Y_i,X_i) \in \mathbb{R}^2,i=1,\cdots,n$`     |     `$Z_i \in \mathbb{R}^d,i=1,\cdots,n$`       |
| 模型      | `$Y_i \vert X_i \sim N(\beta_0+\beta_1 X_i,\sigma^2)$`       |   `$Z_i\sim f(\cdot;\theta)$` |
| 未知参数 |  `$(\beta_0,\beta_1,\sigma^2)$`      |    `$\theta$` |
| 参数空间 | `$\beta_0,\beta_1\in \mathbb{R},\sigma^2 \in \mathbb{R}_+$`   | `$\theta \in \Theta \subset \mathbb{R}^p$` |

---

### 最大似然估计量

- 独立同分布样本 `$Z_i,i=1,\cdots,n$`；随机向量 `$Z_i \in \mathcal{Z} \subset \mathbb{R}^d$`。

- 模型设定正确：
  - 记 `$Z_i$` 的概率密度函数为 `$g$`； 
  - 假设 `$g$` 属于特定参数分布，即 `$g=f(\cdot;\theta_0)$`，其中 `$\theta_0 \in \Theta \subset \mathbb{R}^p$`。
  - 未知参数 `$\theta_0$` 需要根据样本估计出。

- 似然函数和对数似然函数 
`$$L(\theta) = \prod_{i=1}^n f(Z_i;\theta)$$`
`$$l(\theta) = \log L(\theta)= \sum_{i=1}^n \log f(Z_i;\theta)。$$`

- MLE估计量 `$\hat \theta = \arg \max_{\theta \in \Theta} l(\theta)$`。

---

### MLE估计量的计算

- 与前面线性回归例子不同，一般来说，MLE估计量 `$\hat \theta$` 没有显示解析解，其计算依赖于对 `$l(\theta)$` 进行数值优化寻找最大值。根据定义， `$\hat \theta$` 满足 $$ \frac{\partial l(\hat \theta)}{\partial \theta}=0。$$

- 牛顿迭代算法（Newton-Raphson algorithm）
  
  给定 `$\theta^\dagger$`，将上式在 `$\theta^\dagger$` 处进行一阶泰勒展开可得
  $$
  0 = \frac{\partial l(\hat \theta)}{\partial \theta} \approx \frac{\partial l(\theta^\dagger)}{\partial \theta}+\frac{\partial^2 l(\theta^\dagger)}{\partial \theta \partial \theta^\top}(\hat \theta - \theta^\dagger)，
  $$
  经整理后，我们有
  $$
  \hat \theta \approx \theta^\dagger - \left[\frac{\partial^2 l(\theta^\dagger)}{\partial \theta \partial \theta^\top}\right]^{-1} \cdot \frac{\partial l(\theta^\dagger)}{\partial \theta} = \theta^\dagger + J(\theta^\dagger)^{-1} U(\theta^\dagger)。
  $$
  
  
---

### MLE估计量的计算

- 得分向量（score vector） 
`$$U(\theta) = \frac{\partial l}{\partial \theta}(\theta)，$$` 
即对数似然函数关于未知参数的一阶偏导数。对任意给定的 `$\theta$`，得分向量 `$U(\theta)$` 是一个 `$p$` 维向量。

- 可观测信息矩阵（observed information matrix） 
`$$J(\theta) = -\frac{\partial^2 l}{\partial \theta \partial \theta^\top}(\theta)，$$` 
即对数似然函数关于未知参数的二阶偏导数的相反数。对任意给定的 `$\theta$`，可观测信息矩阵 `$J(\theta)$` 是一个 `$p\times p$` 维的矩阵。

---

### MLE估计量的计算

- 牛顿迭代算法的具体步骤
  - 确定一个初始值，记作 `$\theta^{(0)}$`；
  - 对 `$k\ge 0$`，给定第 `$k$` 步的值 `$\theta^{(k)}$`，通过 `$\theta^{(k+1)} = \theta^{(k)}+J(\theta^{(k)})^{-1}U(\theta^{(k)})$` 计算第 `$k+1$` 步的值；
  - 如果 `$\left|l(\theta^{(k+1)})-l(\theta^{(k)})\right|<\epsilon$`则停止迭代，其中 `$\epsilon$` 是预先设定的一个极小的数。

- 牛顿迭代法只能保证找到局部最优（最大或最小）解。

- 虽然很多情况下对数似然函数 `$l(\theta)$` 可被证明是严格凹函数（意味着牛顿迭代法找到的解就是全局最大解），但是在一般情况下需要注意：
  - 给定一个初始值并通过牛顿迭代找出一个局部最优解后，计算可观测信息矩阵，若该矩阵是一个正定矩阵，则说明该最优解是一个局部最大解；
  - 多次变更初始值，通过比较不同初始值所得出的局部最大解确定全局最大解。

---

### MLE估计量的计算

- 一般地，当 `$\Theta$` 存在约束条件时，可考虑令 `$\theta = T(\gamma)$`，其中 `$T(\cdot)$` 是一对一单调增函数且 `$\gamma \in \mathbb{R}^p$`。通过参数变换，可以在 `$\gamma$` 尺度上最大化似然函数 `$l(T(\gamma))$`，这样可使带约束条件的数值优化转变成不受约束的数值优化，防止牛顿迭代过程中参数取值落到参数空间之外。

- 常用参数变换：考虑一个未知参数的情形，即 `$\theta \in \Theta \subset \mathbb{R}$` 
  - `$\Theta = [0,+\infty)$`，可令 `$\theta = \exp(\gamma)$`，其中 `$\gamma \in \mathbb{R}$`。
  - `$\Theta=[0,1]$`，可令 `$\theta = \frac{\exp(\gamma)}{1+\exp(\gamma)}$`，其中 `$\gamma \in \mathbb{R}$`。

- 步长调整：牛顿迭代过程中，可将步长 `$J(\theta^{(k)})^{-1}U(\theta^{(k)})$` 替换为 `$cJ(\theta^{(k)})^{-1}U(\theta^{(k)})$`，在每一次迭代里选择最优的 `$c$` 以最大化对数似然函数，这样可以缩短迭代步数。

---

### MLE估计量的一致性

- 目标：证明 `$\hat \theta \to^p \theta_0$`。

- 假设条件，包括但不限于：
  - **（真实参数）** `$\theta_0$` 位于参数空间 `$\Theta$` 内部，即没有处于 `$\Theta$` 的边界；
  
  - **（模型）** 1. 不同参数对应不同模型，即 `$\theta \ne \theta' \Rightarrow f(\cdot;\theta) \ne f(\cdot;\theta')$`；2. 对于所有参数 `$\theta \in \Theta$`，概率密度函数 `$f(\cdot;\theta)$` 具有相同支撑集，概率密度函数 `$f(\cdot;\theta)$` 的支撑集是集合 `$\{z \in \mathcal{Z} :f(z;\theta)>0 \}$`；
  
  - **（估计量）** MLE估计量 `$\hat \theta$` 是最大化对数似然函数 `$l(\theta)$` 的唯一解。

- 简要证明思路。

对于任意参数 `$\theta$`，考虑 `$n^{-1}l(\theta)$`。根据大数定律，我们有
  $$
  \frac{1}{n}l(\theta) = \frac{1}{n}\sum_{i=1}^n \log f(Z_i;\theta) \to^p E\left[ \log f(Z_i;\theta)\right]。
  $$
---

### MLE估计量的一致性

- 简要证明思路（续）

根据MLE估计量的定义， `$\hat \theta$` 最大化 `$n^{-1}l(\theta)$`，如果我们能够证明 `$\theta_0$` 最大化 `$E\left[ \log f(Z_i;\theta)\right]$`，那么就能在一定假设条件下推出 `$\hat \theta \to^p \theta_0$`。事实上，
`$$\begin{aligned} 
  & E\left[ \log f(Z_i;\theta)\right] - E\left[ \log f(Z_i;\theta_0)\right]\\
  =& E\left[\log \frac{f(Z_i;\theta)}{f(Z_i;\theta_0)} \right] \le \log E\left[ \frac{f(Z_i;\theta)}{f(Z_i;\theta_0)} \right] \\ =& \log \int \frac{f(z;\theta)}{f(z;\theta_0)}f(z;\theta_0)dz = \log \int f(z;\theta)dz\\ =& 0，
\end{aligned}$$`
其中，第二行的小于等于关系是根据Jensen不等式得出的。所以， `$E\left[ \log f(Z_i;\theta)\right]$` 在 `$\theta_0$`处取得最大值，从而完成证明。

*Jensen不等式*：对于随机变量 `$W$` 和凹函数 `$h(\cdot)$`，有 `$E[h(W)]\le h(E[W])$`；对于凸函数 `$h(\cdot)$`，有 `$E[h(W)]\ge h(E[W])$`。

---

### MLE估计量的渐近正态分布

- 目标：证明 `$I(\theta_0)^{1/2}(\hat \theta-\theta_0) \to^d N(0,I_p)$`，其中 `$I(\theta)$` 是费雪信息矩阵， `$N(0,I_p)$` 代指 `$p$` 维标准正态分布。换句话说，当样本容量较大时， `$\hat \theta$` 近似服从 `$N(\theta_0,I(\theta_0)^{-1})$` 分布。

- 费雪信息矩阵（Fisher information matrix）
$$
I(\theta) = E\left[J(\theta) \right]=-E\left[\frac{\partial^2 l}{\partial \theta \partial \theta^\top}(\theta)\right]
$$
是可观测信息矩阵的期望。

- 额外假设条件：
  - **（模型光滑度）** 存在 `$\theta_0$` 的一个邻域 `$\mathcal{B}$`，在该邻域内对数似然函数 `$l(\theta)$` 关于 `$\theta$` 的前三阶导数几乎必然存在，且对于 `$\theta \in \mathcal{B}$` 和 `$r,s,t=1,\cdots,p$`， 
`$n^{-1}E\left\{|\partial^3l(\theta)/\partial \theta_r \partial \theta_s \partial \theta_t| \right\}$` 一致有界。
  - **（费雪信息正定且有限）** 在邻域 `$\mathcal{B}$` 内，费雪信息矩阵有限且是正定的。

---

### MLE估计量的渐近正态分布

- 简要证明思路
  
  根据MLE估计量的定义可知 `$U(\hat \theta)=0$`。将 `$U(\hat \theta)$` 在 `$\theta_0$` 附近一阶泰勒展开：
$$
0=U(\hat \theta) \approx U(\theta_0) + \frac{\partial^2 l(\theta_0)}{\partial \theta \partial \theta^\top}(\hat \theta-\theta_0),
$$
经整理后有
$$
\hat \theta-\theta_0 \approx J(\theta_0)^{-1}U(\theta_0)。
$$
那么
`$$\begin{aligned}
I(\theta_0)^{1/2}(\hat \theta-\theta_0) &= I(\theta_0)^{1/2}J(\theta_0)^{-1}U(\theta_0) \\
&= I(\theta_0)^{1/2}J(\theta_0)^{-1}I(\theta_0)^{1/2}I(\theta_0)^{-1/2}U(\theta_0)。
\end{aligned}$$`
只需要说明 `$I(\theta_0)^{1/2}J(\theta_0)^{-1}I(\theta_0)^{1/2} \to^p I_p$` 并且 `$I(\theta_0)^{-1/2}U(\theta_0) \to^d N(0,I_p)$` 即可。

---

### MLE估计量的渐近正态分布

- 简要证明思路（续）
  
  考虑可观测信息矩阵
`$$\begin{aligned}
J(\theta_0) &=n\cdot \frac{1}{n} \left[-\frac{\partial^2 l}{\partial \theta \partial \theta^\top}(\theta_0)\right]= n\cdot \frac{1}{n}\sum_{i=1}^n -\frac{\partial^2}{\partial \theta \partial \theta^\top}\log f(Z_i;\theta_0)\\
&\to^p n E\left[-\frac{\partial^2 }{\partial \theta \partial \theta^\top}\log f(Z_i;\theta_0) \right]=I(\theta_0)，
\end{aligned}$$`
因此， `$I(\theta_0)^{1/2}J(\theta_0)^{-1}I(\theta_0)^{1/2} \to^p I_p$`。
接下来考虑 `$n^{-1/2}U(\theta_0)=n^{-1/2} \frac{\partial l(\theta_0)}{\partial \theta}=n^{-1/2}\sum_{i=1}^n \frac{\partial}{\partial \theta} \log f(Z_i;\theta_0)$`。由于
`$$\begin{aligned}
E\left[\frac{\partial}{\partial \theta} \log f(Z_i;\theta_0) \right]&=\int f(z;\theta_0) \frac{\partial}{\partial \theta} \log f(z;\theta_0)dz \\ &=\int \frac{\partial}{\partial \theta}f(z;\theta_0)dz = \frac{\partial}{\partial \theta}\int f(z;\theta_0)dz=0，
\end{aligned}$$`

---
### MLE估计量的渐近正态分布

- 简要证明思路（续）

根据中心极限定理，我们有
`$$n^{-1/2}U(\theta_0) \to^d N\left(0,E\left[\frac{\partial}{\partial \theta} \log f(Z_i;\theta_0) \frac{\partial}{\partial \theta^\top} \log f(Z_i;\theta_0)\right] \right)。$$` 
如果上式中的渐近方差等于 `$n^{-1}I(\theta_0)$`，即 `$n^{-1/2}U(\theta_0) \to^d N\left(0,n^{-1}I(\theta_0) \right)$`，那么
`$$\begin{aligned}
I(\theta_0)^{-1/2}U(\theta_0) &= n^{1/2}I(\theta_0)^{-1/2}n^{-1/2}U(\theta_0) \\
&\to^d n^{1/2}I(\theta_0)^{-1/2}N\left(0,n^{-1}I(\theta_0) \right)\\
&= N(0,I_p)。
\end{aligned}$$`
合并整理上述结果完成证明。

---
### MLE估计量的渐近正态分布

- 简要证明思路（续）

现在我们验证
`$$\begin{aligned}
& n^{-1} I(\theta_0)=n^{-1}\sum_{i=1}^n E\left[-\frac{\partial^2}{\partial \theta \partial \theta^\top}\log f(Z_i;\theta_0) \right]\\
=&n^{-1}\sum_{i=1}^n E\left[\frac{1}{f(Z_i;\theta_0)^2}\frac{\partial f(Z_i;\theta_0)}{\partial \theta}\frac{\partial f(Z_i;\theta_0)}{\partial \theta^\top}-\frac{1}{f(Z_i;\theta_0)}\frac{\partial^2 f(Z_i;\theta_0)}{\partial \theta \partial \theta^\top} \right]\\
=&E\left[\frac{\partial}{\partial \theta} \log f(Z_i;\theta_0) \frac{\partial}{\partial \theta^\top} \log f(Z_i;\theta_0)\right]-n^{-1}\sum_{i=1}^n \int \frac{\partial^2 f(z;\theta_0)}{\partial \theta \partial \theta^\top}dz。
\end{aligned}$$`
又根据 `$\int f(z;\theta_0)dz=1$`，等式两边对 `$\theta$` 同时求二阶偏导数可得
`$$\frac{\partial^2}{\partial \theta \partial \theta^\top}\int f(z;\theta_0)dz= \int \frac{\partial^2 f(z;\theta_0)}{\partial \theta \partial \theta^\top}dz=0，$$`
从而完成验证。

---

### MLE估计量渐近方差的直觉解释

为什么费雪信息矩阵能够刻画MLE估计量的渐近方差？为什么渐近方差与费雪信息矩阵是求逆的关系？

- 由前述证明可知 `$I(\theta_0)=E[U(\theta_0)U(\theta_0)^\top]$`。得分向量 `$U(\theta_0)$` 衡量了当参数 `$\theta$` 微小偏离 `$\theta_0$`后对数似然函数的变化，那么 `$I(\theta_0)$` 可以理解为对数似然函数的平均（对样本 `$Z_i$`）平方变化。

- 又由定义可知 `$I(\theta_0)=E[J(\theta_0)]$`，即对数似然函数在 `$\theta_0$` 附近的平均（对样本 `$Z_i$`）曲率，曲率越大，变化越陡峭。

- 无论从上述两个维度去理解费雪信息，我们都能得出：费雪信息越“大”，对数似然函数在真实参数 `$\theta_0$` 附近的区分度越高，越能识别出真实参数，因此估计得越精准，估计量的方差越小。

---

### Cramer-Rao下界：MLE估计量的渐近有效性

- ***（Cramer-Rao下界）*** 给定独立同分布样本 `$Z_i, i=1,\cdots,n$`，假设 `$Z_i$` 的概率密度函数 `$g$` 可以由参数分布 `$f(\cdot;\theta_0)$`表示，其中 `$\theta_0 \in \Theta$`。记 `$\tilde \theta$` 为真实参数 `$\theta_0$` 的一个***无偏估计量***，在一些正则假设条件下，我们有
$$
Var(\tilde \theta) \ge I(\theta_0)^{-1}，
$$
其中 `$Var(\tilde \theta)$` 是估计量 `$\tilde \theta$` 的方差矩阵，此处“ `$\ge$` "应理解为 `$Var(\tilde \theta)-I(\theta_0)^{-1}$` 是半正定矩阵。

- Cramer-Rao下界告诉了我们一个无偏估计量能做到的“最好”程度。

- MLE估计量是渐近无偏的，且渐近方差正好是 `$I(\theta_0)^{-1}$`，因此我们称MLE是渐近有效的。

---

### Cramer-Rao下界：MLE估计量的渐近有效性

- ***（Cramer-Rao下界的拓展）*** 记 `$\tilde \theta$` 为真实参数 `$\theta_0$` 的一个估计量（可能无偏也可能有偏），那么在一些正则假设条件下，我们有
$$
MSE(\tilde \theta) \ge I(\theta_0)^{-1}，
$$
其中 `$MSE(\tilde \theta)$`是是估计量 `$\tilde \theta$` 的均方误差。

- 均方误差的定义及其分解
`$$\begin{aligned}
MSE(\tilde \theta) &= E\left[(\tilde \theta - \theta_0)(\tilde \theta - \theta_0)^\top \right]\\
&= (E \tilde \theta - \theta_0)(E\tilde \theta - \theta_0)^\top + E\left[(\tilde \theta -E \tilde \theta)(\tilde \theta -E\tilde \theta)^\top \right]\\
&= Bias(\tilde \theta)Bias(\tilde \theta)^\top +Var(\tilde \theta)，
\end{aligned}$$`
其中， `$Bias(\tilde \theta)$` 代指偏差。

---

### OLS的BLUE性质与MLE的渐近有效

- 回到简单线性回归
$$
Y_i = \beta_0 +\beta_1 X_i +e_i, \quad\quad i=1,\cdots,n
$$
假设残差服从 `$t(\sigma,\nu)$` 分布，其中 `$\sigma$` 是尺度参数， `$\nu$` 是自由度参数。

- 该模型适用于 `$Y|X$` 的分布是厚尾的情形，即因变量 `$Y$` 出现较多异常值的情形。

- BLUE性质告诉我们OLS估计量是最有效的，而MLE又是渐近有效的，到底哪个估计量最有效？两个结论矛盾吗？

- 不矛盾。对于该回归模型，假设残差 `$t$` 分布的MLE估计量更加有效，但这并不代表BLUE性质不对，因为MLE估计量不是线性估计量。

---

### 基于MLE的统计推断：似然比统计量

- 针对未知总体参数 `$\theta_0$`，考虑简单假设检验：
`$$H_0: \theta_0 = \theta_\star \quad \quad \quad
H_1: \theta_0 \ne \theta_\star。$$`

- 一个自然的比较标准是似然比：
`$$LR(\theta_\star) = \frac{L(\theta_\star)}{\sup_{\theta \in \Theta}L(\theta)}=\frac{L(\theta_\star)}{L(\hat \theta)}。$$`
似然比 `$LR(\theta_\star)\in [0,1]$`，当原假设为真时， `$LR(\theta_\star)$` 应接近1，当原假设为假时， `$LR(\theta_\star)$` 应远离1。那么如何度量是“接近”还是“远离”呢？这需要建立关于似然比的渐近分布。

- （Wilks定理）假设建立MLE估计量渐近正态分布所需要的假设条件成立，当原假设为真时，随着样本容量 `$n \to \infty$`，
`$$W(\theta_\star) = -2\log LR(\theta_\star) = 2l(\hat \theta)-2l(\theta_\star) \to^d \chi_p^2，$$`
其中自由度 `$p$` 是未知参数的个数。

---

### 基于MLE的统计推断：似然比统计量

- `$W(\theta_\star)\to^d \chi_p^2$` 的证明思路
  
  我们将对数似然函数 `$l(\theta_\star)$` 在 `$\hat \theta$` 处二阶泰勒展开得到
`$$l(\theta_\star) = l(\hat\theta)+(\theta_\star-\hat \theta)^\top U(\hat \theta)-\frac{1}{2}(\theta_\star-\hat \theta)^\top J(\tilde \theta)(\theta_\star-\hat \theta)，$$`
其中， `$\tilde \theta$` 位于 `$\theta_\star$` 和 `$\hat \theta$` 之间。当原假设为真，即 `$\theta_\star=\theta_0$` 时， `$\hat \theta \to^p \theta_\star$` 意味着 `$\tilde \theta \to^p \theta_\star$`。又由于 `$U(\hat \theta)=0$`，有
`$$\begin{aligned}
W(\theta_\star)&=(\hat \theta-\theta_\star)^\top J(\tilde \theta)(\hat \theta-\theta_\star)\\
&=\left[I(\theta_\star)^{1/2}(\hat \theta-\theta_\star)\right]^\top \left[I(\theta_\star)^{-1/2} J(\tilde \theta)I(\theta_\star)^{-1/2}\right]\left[I(\theta_\star)^{1/2}(\hat \theta-\theta_\star)\right]。
\end{aligned}$$`
根据MLE的渐近正态分布可知， `$I(\theta_\star)^{1/2}(\hat \theta-\theta_\star) \to^d N(0,I_p)$`。由于 `$\tilde \theta \to^p \theta_\star$` 且 `$J(\theta_\star)\to^p I(\theta_\star)$`，那么 `$J(\tilde \theta)\to^p I(\theta_\star)$`，进而 `$I(\theta_\star)^{-1/2} J(\tilde \theta)I(\theta_\star)^{-1/2}\to^p I_p$`。汇总以上结果得出 `$W(\theta_\star)\to^d \chi_p^2$`。

---

### 基于MLE的统计推断：似然比统计量

- 当原假设为假时， `$LR(\theta_\star)$` 远离1，使得 `$W(\theta_\star)$` 沿着实数轴向右变动。

- 给定显著性水平 `$\alpha$`，临界值为 `$\chi_p^2$` 分布的 `$1-\alpha$` 分位数，记作 `$\chi^2_{p,1-\alpha}$`。当 `$W(\theta_\star)>\chi^2_{p,1-\alpha}$`，拒绝原假设；否则不能拒绝原假设。

---

### 模型设定错误：拟最大似然估计（Quasi-MLE）

- 独立同分布样本 `$Z_i,i=1,\cdots,n$`，其中 `$Z_i$` 的概率密度函数为 `$g$`。

- 假设参数模型为 `$f_\theta\equiv f(\cdot;\theta)$`，其中参数 `$\theta \in \Theta \subset \mathbb{R}^p$`。但是真实密度函数 `$g$` 不能由该类参数模型表示，即对所有 `$\theta \in \Theta$` 都有 `$g\ne f_\theta$`。

- 拟最大似然估计量（QMLE）：在“错误”的参数模型基础上最大化对数似然得出的估计量，即
`$$\check \theta = \arg \max_{\theta \in \Theta} l(\theta)，$$` 其中
`$$l(\theta)=\sum_{i=1}^n \log f(Z_i;\theta)。$$`

- 由于真实概率密度函数不能够由所假设的参数模型表示，那么就不存在“真实参数”，那么QMLE到底在估计什么？

---

### QMLE在干什么？——矮子里头拔将军

- Kullback-Leibler（KL）散度

$$
KL(g \parallel f_\theta)=\int g(z)\log \frac{g(z)}{f(z;\theta)}dz。
$$
- 一些性质

- `$KL(g\parallel f_\theta)\ge 0$`，当且仅当 `$g=f_\theta$`取等号。
  - KL散度不是距离：非对称 `$KL(g\parallel f_\theta)\ne KL(f_\theta \parallel g)$`且不满足三角不等式。
  - KL散度度量了使用 `$f_\theta$` 近似 `$g$` 的信息损失。

- 既然模型假设错误了，那么就将错误降至最低，于是有了 
`$$\theta_g=\arg\min_{\theta \in \Theta}KL(g\parallel f_\theta)。$$`
注意到 `$KL(g\parallel f_\theta)=E\left[\log g(Z_i)\right]-E\left[\log f(Z_i;\theta)\right]$`，换句话说 `$\theta_g=\arg\max_{\theta \in \Theta} E\left[\log f(Z_i;\theta)\right]$`。

---

### QMLE的理论性质

- QMLE在估计 `$\theta_g$`，即 `$\check \theta \to^p \theta_g$`

根据大数定律，
`$$\frac{1}{n}l(\theta) = \frac{1}{n}\sum_{i=1}^n \log f(Z_i;\theta) \to^p E\left[\log f(Z_i;\theta)\right]，$$`
由于 `$\check \theta$` 最大化 `$\frac{1}{n}l(\theta)$`，而 `$\theta_g$`最大化 `$E\left[\log f(Z_i;\theta)\right]$`，那么在一些正则条件下就可以得出 `$\check \theta \to^p \theta_g$`。

---

### QMLE的理论性质

- QMLE在大样本下依然是近似正态分布，但是渐近方差改变了。

根据最大化对数似然函数的一阶条件，有 `$\sum_{i=1}^n \frac{\partial \log f(Z_i;\check \theta)}{\partial \theta}=0$`，将等式左边在 `$\theta_g$` 处泰勒展开，有
`$$0=\sum_{i=1}^n \frac{\partial \log f(Z_i;\check \theta)}{\partial \theta} \approx \sum_{i=1}^n \frac{\partial \log f(Z_i;\theta_g)}{\partial \theta}+\sum_{i=1}^n \frac{\partial^2 \log f(Z_i;\theta_g)}{\partial \theta \partial \theta^\top}(\check \theta -\theta_g)，$$`
经整理后得到
`$$\sqrt{n}(\check \theta - \theta_g) \approx \left[-\frac{1}{n}\sum_{i=1}^n \frac{\partial^2 \log f(Z_i;\theta_g)}{\partial \theta \partial \theta^\top}\right]^{-1}\left[\frac{1}{\sqrt{n}}\sum_{i=1}^n \frac{\partial \log f(Z_i;\theta_g)}{\partial \theta}\right]。$$`

---

### QMLE的理论性质
- QMLE的渐近正态分布（续）
  
  我们类似地定义得分向量 `$U(\theta_g)$` 和可观测信息矩阵 `$J(\theta_g)$`：
`$$\begin{aligned}
U(\theta_g) = \sum_{i=1}^n \frac{\partial}{\partial \theta}\log f(Z_i;\theta_g)\\
J(\theta_g) = -\sum_{i=1}^n \frac{\partial^2}{\partial \theta\partial \theta^\top}\log f(Z_i;\theta_g)。
\end{aligned}$$`
进一步地，我们定义
`$$\begin{aligned}
I_g(\theta_g) &= E\left[J(\theta_g)\right]= n \int -\frac{\partial^2 \log f(z;\theta_g)}{\partial \theta \partial \theta^\top}g(z)dz\\
K_g(\theta_g) &=E\left[U(\theta_g)U(\theta_g)^\top\right]= n \int \frac{\partial \log f(z;\theta_g)}{\partial \theta}\frac{\partial \log f(z;\theta_g)}{\partial \theta^\top}g(z)dz。
\end{aligned}$$`

---
### QMLE的理论性质
- QMLE的渐近正态分布（续）

由于 `$\theta_g$`最大化 `$E\left[\log f(Z_i;\theta)\right]=\int \log f(z;\theta)g(z)dz$`，其一阶条件为
`$$\frac{\partial}{\partial \theta}\int \log f(z;\theta_g)g(z)dz = \int \frac{\partial}{\partial \theta}\log f(z;\theta_g)g(z)dz\\=E\left[\frac{\partial \log f(Z_i;\theta_g)}{\partial \theta}\right]=E\left[U(\theta_g)\right]=0，$$`
那么 `$K_g(\theta_g)$` 就是得分向量 `$U(\theta_g)$` 的方差矩阵。

根据大数定律
`$$-\frac{1}{n}\sum_{i=1}^n \frac{\partial^2 \log f(Z_i;\theta_g)}{\partial \theta \partial \theta^\top} \to^p \frac{1}{n}I_g(\theta_g)。$$`

---
### QMLE的理论性质
- QMLE的渐近正态分布（续）

根据中心极限定理，
`$$\frac{1}{\sqrt{n}}\sum_{i=1}^n \frac{\partial \log f(Z_i;\theta_g)}{\partial \theta} \to^d N(0,n^{-1}K_g(\theta_g))。$$`
合并上述结果有
`$$\begin{aligned}
\sqrt{n}(\check \theta - \theta_g) & \to^d n I_g(\theta_g)^{-1}\cdot N(0,n^{-1}K_g(\theta_g)) \\
&= N(0,n I_g(\theta_g)^{-1}K_g(\theta_g)I_g(\theta_g)^{-1})，
\end{aligned}$$`
也就是说大样本下（ `$n\to \infty$` ），QMLE估计量的渐近分布为
$$ \check \theta \sim N(\theta_g,I_g(\theta_g)^{-1}K_g(\theta_g)I_g(\theta_g)^{-1})。$$

- QMLE的渐近方差是典型的“三明治”形式。当模型设定正确时， `$\theta_g = \theta_0$` 并且 `$K_g(\theta_g)=I_g(\theta_g)$`，渐近方差形式退化到了模型设定正确的情形。

---

### MLE/QMLE的GMM解释

- 无论是MLE还是QMLE，待估计的未知参数都满足矩条件：
  - MLE： `$E\left[\frac{\partial}{\partial \theta} \log f(Z_i;\theta_0) \right]=0$`；
  - QMLE： `$E\left[\frac{\partial}{\partial \theta} \log f(Z_i;\theta_g) \right]=0$`，
  
  所对应的样本矩是 `$n^{-1}U(\theta)$`。

- 基于上述矩条件构造GMM估计量
`$$\hat \theta_{GMM} = \arg \min_{\theta \in \Theta} U(\theta)^\top W U(\theta)。$$`
由于有 `$p$` 个未知参数 `$p$` 个矩条件，属于恰好识别类型，该GMM本质上就是矩估计。又因为二次型 `$U(\theta)^\top W U(\theta)\ge 0$`，因此 `$\hat \theta_{GMM}$` 满足 `$U(\hat \theta_{GMM})=0$`，此时权重矩阵的选择无关。MLE/QMLE估计量同时可以理解为GMM估计量。

- *在GMM框架下，估计量的大样本近似分布是什么？与MLE/QMLE的大样本近似分布是否一致？*

---

### 模型选择

- 在实际应用中，通常我们无法得知生成数据的真实模型 `$g$`，如何从一系列候选模型中选择一个**最优**的模型呢？

- 如何衡量“最优”？——KL散度。记某一候选模型为 `$f_\theta$`，其中 `$\theta \in \Theta$`，我们将该模型视作一个近似模型，可能包括也有可能不包括真实模型 `$g$`，最优参数记为 `$\theta_g$`，相应的QMLE估计量记为 `$\check \theta$`。定义
`$$KL(g \parallel f_{\check \theta}) = \int g(z)\log \frac{g(z)}{f(z;\check \theta)}dz。$$`
在一系列模型中，选择 `$KL(g \parallel f_{\check \theta})$` 最小的模型作为最优模型。

- 奥卡姆剃刀原理（Ockham's razor）
  > 如无必要，勿增实体。
  > It is vain to do with more what can be done with fewer.
  
  在模型选择中，如果几个模型拟合数据相同地好，那么选择那个最简单的模型。

---

### 模型选择

- 将KL散度进一步整理，有 
`$$KL(g \parallel f_{\check \theta})=\int g(z)\log g(z)dz - \int g(z)\log f(z;\check \theta)dz,$$`
其中第一项与候选模型无关可以忽略，那么最小化KL散度等同于最大化
`$$K(f) \equiv \int g(z)\log f(z;\check \theta)dz。$$`
- 可以将 `$K(f)$` 写为期望的形式
`$$K(f) = E_Z\left[\log f(Z;\check \theta) \right]，$$`
其中 `$Z$` 是与样本 `$Z_i, i=1,\cdots,n$` 独立同分布的一个随机向量，因此 `$Z$` 与 `$\check \theta$` 是独立的，上式是对 `$Z$` 求期望。

---
### 交叉验证法（cross validation）

- 如何合理地估计 `$K(f)$` 是构建模型选择方法的重点。一个直接的估计量是 `$$\bar K(f) = n^{-1}l(\check \theta)=n^{-1}\sum_{i=1}^n \log f(Z_i;\check \theta)。$$`

估计量 `$\bar K(f)$` 存在严重偏差，直观的原因是样本 `$Z_i, i=1,\cdots,n$` 被使用了两次，一次是计算QMLE估计量，一次是利用样本平均估计期望，那么 `$Z_i$` 与 `$\check \theta$` 不是独立的，不符合 `$K(f)$`的定义。
  
- 为了解决 `$\bar K(f)$` 中 `$Z_i$` 与 `$\check \theta$` 不独立的问题，交叉验证法将样本分为两个部分：一部分是 `$Z_i$`，另外一部分是其余所有样本点 `$Z_j,j\ne i$`，用于QMLE估计，估计量记作 `$\check \theta_{-i}$`。交叉验证法写作
`$$CV(f)=n^{-1}\sum_{i=1}^n \log f(Z_i;\check \theta_{-i})，$$`
显然 `$Z_i$` 与 `$\check \theta_{-i}$` 是独立的。选取 `$CV(f)$` 最大的模型。

---
### 赤池信息准则（AIC，Akaike Information Criterion）

- 交叉验证法的一个缺点是，对于每个候选模型都需要计算QMLE估计量 `$n$` 次，运算量大。另外一个思路是对 `$\bar K(f)$` 做偏差修正，这样对于每个模型只需要计算一次QMLE，这引出了AIC。

- `$\bar K(f)$` 相对 `$K(f)$` 的偏差 
  
  首先，将 `$\log f(z;\check \theta)$` 在 `$\theta_g$` 处泰勒展开，有
`$$\begin{aligned}
K(f) &= \int g(z)\log f(z;\check \theta)dz \\
&\approx \int g(z)\left[\log f(z;\theta_g)+(\check \theta - \theta_g)^\top \frac{\partial \log f(z;\theta_g)}{\partial \theta} \right. \\ & \left.+\frac{1}{2}(\check \theta - \theta_g)^\top \frac{\partial^2 \log f(z;\theta_g)}{\partial \theta\partial \theta^\top}(\check \theta - \theta_g)\right]dz\\
&= \int g(z)\log f(z;\theta_g)dz -\frac{1}{2n}(\check \theta - \theta_g)^\top I(\theta_g)(\check \theta - \theta_g)。
\end{aligned}$$`

---
### 赤池信息准则（AIC，Akaike Information Criterion）

- `$\bar K(f)$` 相对 `$K(f)$` 的偏差（续）

其次，将 `$l(\theta_g)$` 在 `$\check \theta$` 处泰勒展开，有
`$$\begin{aligned}
\bar K(f) &= n^{-1}\sum_{i=1}^n \log f(Z_i;\check \theta) = n^{-1} l(\check \theta)\\
&\approx n^{-1}\left[l(\theta_g) - (\theta_g-\check \theta)^\top U(\check \theta)+\frac{1}{2}(\theta_g-\check \theta)^\top J(\check \theta)(\theta_g-\check \theta)\right]\\
&= \frac{1}{n}l(\theta_g)+\frac{1}{2n}(\theta_g-\check \theta)^\top J(\check \theta)(\theta_g-\check \theta)\\
&\approx \frac{1}{n}l(\theta_g)+\frac{1}{2n}(\check \theta-\theta_g)^\top I(\theta_g)(\check \theta-\theta_g)，
\end{aligned}$$`
其中最后一步利用了 `$J(\check \theta) \approx I(\theta_g)$`。

---

### 赤池信息准则（AIC，Akaike Information Criterion）

- `$\bar K(f)$` 相对 `$K(f)$` 的偏差（续）

`$$\begin{aligned}
& E \left[\bar K(f) - K(f)\right]\\ =& E\left[\frac{1}{n}l(\theta_g)-\int g(z)\log f(z;\theta_g)dz+\frac{1}{n}(\check \theta-\theta_g)^\top I(\theta_g)(\check \theta-\theta_g)\right]\\
=& \frac{1}{n}E\left[(\check \theta-\theta_g)^\top I(\theta_g)(\check \theta-\theta_g)\right]= \frac{1}{n}E\left[\text{trace}\left[(\check \theta-\theta_g)^\top I(\theta_g)(\check \theta-\theta_g)\right]\right]\\
=& \frac{1}{n}E\left[\text{trace}\left[ I(\theta_g)(\check \theta-\theta_g)(\check \theta-\theta_g)^\top\right]\right]= \frac{1}{n}\text{trace}\left[E\left[ I(\theta_g)(\check \theta-\theta_g)(\check \theta-\theta_g)^\top\right]\right]\\
=& \frac{1}{n}\text{trace}\left[I(\theta_g)E\left[ (\check \theta-\theta_g)(\check \theta-\theta_g)^\top\right]\right]=\frac{1}{n}\text{trace}\left[I(\theta_g)Var( \check \theta)\right]\\
=& \frac{1}{n}\text{trace}\left[K_g(\theta_g)I(\theta_g)^{-1}\right]。
\end{aligned}$$`
- 经偏差校正后， `$K(f)$` 的一个估计量是
`$$\hat K(f) = n^{-1} l(\check \theta) - n^{-1}\text{trace}\left[K_g(\theta_g)I(\theta_g)^{-1}\right]。$$`

---

### 赤池信息准则（AIC，Akaike Information Criterion）
- 选择模型 `$f$` 最小化 `$KL(g \parallel f_{\check \theta})$` 等同于最大化 `$\hat K(f)$`，也就是最小化 `$-2n\hat K(f)$`，即$$-2 l(\check \theta)+2\text{trace}\left[K_g(\theta_g)I(\theta_g)^{-1}\right]。$$

- 当候选模型 `$f_\theta$` 包含真实模型 `$g$` 时，我们有 `$K_g(\theta_g) = I(\theta_g)$`，因此 `$\text{trace}\left[K_g(\theta_g)I(\theta_g)^{-1}\right]=\text{trace}[I_p]=p$`。当候选模型不包含真实模型时，我们也可以使用近似 `$\text{trace}\left[K_g(\theta_g)I(\theta_g)^{-1}\right]\approx p$`。

- 在一系列候选模型中最小化AIC
`$$-2 l(\check \theta)+2p，$$`
其中，第一项可以看作样本内的拟合优度，第二项是针对模型复杂度的惩罚项。

- 与交叉验证相比，AIC计算量更小，但是AIC的推导过程使用了大量的近似替换，在可靠性上可能差于交叉验证。

---

### NIC，Network Information Criterion

- NIC可以看作AIC的拓展，它将未知的 `$K_g(\theta_g)$` 和 `$I(\theta_g)$` 替换为相应的估计量，记作
`$$\widehat{K_g(\theta_g)}= \sum_{i=1}^n \frac{\partial \log f(Z_i;\check \theta)}{\partial \theta}\frac{\partial \log f(Z_i;\check \theta)}{\partial \theta^\top}\\
\widehat{I_g(\theta_g)}= - \sum_{i=1}^n \frac{\partial^2 \log f(Z_i;\check \theta)}{\partial \theta \partial \theta^\top}。$$`

- NIC最小化$$-2 l(\check \theta)+2\text{trace}\left[\widehat{K_g(\theta_g)}\widehat{I_g(\theta_g)}^{-1}\right]。$$