本文主要介绍XGBoost和其他传统GBDT的比较的优劣
XGBoost又叫(Newton Boosting Tree)

• GBDT推导: ML——GBDT-梯度提升树-推导过程
• XGBoost推导: ML——XGBoost-推导过程

XGBoost的优点

参考博客: https://www.cnblogs.com/massquantity/p/9794480.html

• XGBoost损失函数是二阶泰勒展开(与牛顿法对应),GBDT是一阶泰勒展开(与梯度下降法对应)

  • 传统 GBDT 在优化时只用到一阶导数信息, 所以传统GBDT也叫 (Gradient Boosting)
  • XGBoost 则对目标函数进行了二阶泰勒展开，同时用到了一阶和二阶导数, 所以XGBoost又叫(Newton Boosting Tree)

• XGBoost加了正则项 ，普通的GBDT没有,所以XGBoost学出来的模型更简单,能防止过拟合,提高模型的泛化性能

• XGBoost Shrinkage(缩减)

  • 每次进行完一次迭代后,将叶子节点的权重乘以该系数(一般叫做eta \(\eta\))
  • 可以理解为这里Shrinkage是将学习速率调小,从而需要的迭代次数增多
  • 减小学习率实际上是减弱了每棵树对整体的影响,从而让后面的树有更多的学习空间
  • 下面的表述还有待确定:
    • 进一步惩罚决策树叶节点的值(惩罚的意思是叶节点越大,惩罚越多,损失函数越大)
    • 对叶节点的惩罚本身可以理解为一个正则化

• 结点分裂的增益计算公式不同

  • 传统 GBDT 一般采用的是最小二乘法作为内部分裂的增益计算指标(因为用的都是回归树)
    • 注意: 这里本文中描述的最小绝对偏差回归(LAD)是第一步损失函数的定义,不是这一步中的信息增益计算
    • 查看源码: sklearn.ensemble.GradientBoostingClassifier在分裂结点时可以选择三种方式:
      • friedman_mse(默认), mean squared error with improvement score by Friedman
      • mse: mean squared error
      • mae: mean absolute error
  • 而 XGBoost 使用的是经过优化推导后的式子
    $$
    \begin{align}
    Gain = \frac{G_L^2}{H_L+ \lambda} + \frac{G_R^2}{H_R+ \lambda} - \frac{(G_L + G_R)^2}{H_L+ H_R + \lambda} - \gamma
    \end{align}
    $$
    • 注意: XGBoost中的信息增益计算形式固定为上面的计算方式,但是具体的值与损失函数的定义相关(因为 \(g_i\) 和 \(h_i\) 的是损失函数的一阶和二阶梯度)

• XGBoost支持自定义的损失函数 ，支持一阶和二阶可导就行

  • 注意,这里的损失函数指的是 \(l(y_i,\hat{y}_i)\)，单个样本预测值与目标值的差异, 也就是单个样本的损失函数
  • 从ML——XGBoost-推导过程中可知:
    • \(g_i = l’(y_i,\hat{y}_i^{t-1})\) 为 \(l(y_i,\hat{y}_i)\) 对 \(\hat{y}_i\) 的一阶导数在 \(\hat{y}_i = \hat{y}_i^{t-1}\) 处的值
    • \(h_i = l’’(y_i,\hat{y}_i^{t-1})\) 是 \(l(y_i,\hat{y}_i)\) 对 \(\hat{y}_i\) 的二阶导数在 \(\hat{y}_i = \hat{y}_i^{t-1}\) 处的值
  • XGBoost中只要损失函数二次可微分即可得到 \(g_i\) 和 \(h_i\)
    • \(g_i\) 和 \(h_i\) 本身与损失函数的定义形式无关, 只要求损失函数二阶可微分即可
  • 只要有了 \(g_i\) 和 \(h_i\) 我们即可
    • 根据预先推导的叶子节点分数表达式 \(w_j^{\star} = -\frac{G_j}{H_j+\lambda}\) 求得叶子结点的分数
    • 根据预先推导的信息增益公式 \(Gain = \frac{G_L^2}{H_L+ \lambda} + \frac{G_R^2}{H_R+ \lambda} - \frac{(G_L + G_R)^2}{H_L+ H_R + \lambda} - \gamma\) 确定分裂特征和分裂点
  • GBDT 损失函数关系一般选择最小二乘回归或者最小绝对偏差回归
    • 最小方差回归 : (Least-Squares Regression, LSR)
      $$\begin{align} L(y,F(x)) = \frac{1}{2}(y-F(x))^{2} \end{align}$$
    • 最小绝对偏差回归 : (Least Absolute Deviation Regression, LAD)
      $$\begin{align} L(y,F(x)) = |y-F(x)| \end{align}$$
    • 查看源码: sklearn.ensemble.GradientBoostingClassifier的损失函数是定义好的, 不能自己定义, 详细如下源码

      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

      LOSS_FUNCTIONS = {'ls': LeastSquaresError, 'lad': LeastAbsoluteError, 'huber': HuberLossFunction, 'quantile': QuantileLossFunction, 'deviance': None, # for both, multinomial and binomial 'exponential': ExponentialLoss, } _SUPPORTED_LOSS = ('deviance', 'exponential') .... if (self.loss not in self._SUPPORTED_LOSS or self.loss not in LOSS_FUNCTIONS): raise ValueError("Loss '{0:s}' not supported. ".format(self.loss))

• XGBoost 借鉴了随机森林的做法，支持列采样 ，不仅能降低过拟合，还能减少计算量，这也是 XGBoost 异于传统 GBDT 的一个特性

  • 列采样: 这借鉴于随机森林中的做法,每棵树不使用所有特征,而是部分特征参与训练
  • 可以减少计算量,同时还能降低过拟合,简直优秀

• XGBoost 有缺失值自动处理 ， 在计算分裂增益时不会考虑带有缺失值的样本，这样就减少了时间开销,在分裂点确定了之后，将带有缺失值的样本分别放在左子树和右子树，比较两者分裂增益，选择增益较大的那一边作为默认分裂方向

  • 进一步理解稀疏数据的支持: [待更新]

• 并行化处理 ：由于 Boosting 本身的特性，传统 GBDT 无法像随机森林那样树与树之间的并行化

  • XGBoost 的并行主要体现在特征粒度上，在对结点进行分裂时，由于已预先对特征排序并保存为block 结构，每个特征的增益计算就可以开多线程进行，极大提升了训练速度

• 剪枝策略不同

  • 传统 GBDT 在损失不再减少时会停止分裂，这是一种预剪枝的贪心策略，容易欠拟合
  • XGBoost采用的是后剪枝的策略，先分裂到指定的最大深度 (max_depth) 再进行剪枝
    • 而且和一般的后剪枝不同， XGBoost 的后剪枝是不需要验证集的[待更新:XGBoost剪枝的策略是怎样的?只依赖信息增益指标吗?]
    • 和博客作者指出的一样,我这里并不是”纯粹的”后剪枝,因为提前设定了最大深度

• 基分类器的选择不同:

  • 传统GBDT中原始论文使用树回归 ，本文见Firedman 1999,后来作者提出可以使用逻辑回归 ，本文见Friedman 2000
  • XGBoost后面的各种损失计算等都包含着树模型的复杂度,叶子节点分类等,所以是只能用CART,不能使用逻辑回归的
  • (从函数空间定义和后面的公式推导来看)XGBoost中基函数只使用CART回归树，不能使用逻辑回归
  • 但是事实上XGBoost的实现中是支持线性分类器作为基分类器的, 参数booster[default='gbtree'],可选为booster=gblinear
    • 使用线性分类器作为基分类器时, XGBoost相当于带有L1正则化和L2正则化的:
      • Logistic回归(分类问题)
      • 线性回归(回归问题)

• 分桶策略算法不同

  • 传统的GBDT分桶时每个样本的权重都是相同的
  • XGBoost中每个样本的权重为损失函数在该样本点的二阶导数(对不同的样本,计算得到的损失函数的二阶导数是不同的), 这里优点AdaBoost的思想,重点关注某些样本的感觉
  • 这里影响的是划分点的位置(我们划分划分点[桶]时都是均匀划分样本到桶里面,当不同样本的权重不同时,每个桶里面的样本数量可能会不同)
  • 下图是一个示例

XGBoost的缺点

注意这里是

• 空间消耗大
  • 因为要在训练之前先对每个特征进行预排序并将结果存储起来，所以空间消耗较大
  • GBDT无需预排序,但是每次重新排序很耗时间
• 速度慢
  • 虽然XGBoost速度比传统 GBDT 快了不少, 但是不如 LightGBM 快, 且 LightGBM 占用内存更低

XGBoost为什么能够并行?

而GBDT是不能并行的，原因是：https://www.136.la/shida/show-187480.html

• GBDT不能并行的原因是没有预排序（XGB的预排序结果会存储到block结构）等，在有了预排序结果后，同一个特征的切割方式可以并行尝试计算增益
• 决策树最耗时间(包括GBDT)的步骤是对特征值的排序
  • 用于确定最佳分割点
• XGBoost训练前,预先对数据进行了排序,称为预排序
  • 将预先排序的结果保存为block结构, 后面迭代的时候重复使用这个结构,从而实现一次排序,多次使用,大大减少计算量
• 这个结构减少计算量的同时还为并行化提供了可能([待更新]实际上不用预排序也能并行的吧?只是每次需要先使用一个单一线程排序,然后再多个线程并行?只是不够并行)
  • 进行结点的分裂时,需要计算每个特征的增益,然后选择增益最大的那个特征分裂
  • 这里我们可以同时使用多个线程计算不同特征的增益, 从而实现并行化
• 总结为三方面的并行, ([待更新]但是具体实现了哪些并行不确定)
  • 同一层级的结点间每个结点的分裂可以并行
  • 同一个结点内部不同特征增益的计算可以并行
  • 同一个结点同一个特征的不同分裂点的增益计算可以并行

GBDT为什么不能自定义损失函数?

GBDT为什么不能像XGBoost一样自定义损失函数?

• 查看sklearn.ensemble.GradientBoostingClassifier的源码发现, 确实不支持自定义的损失函数
• [待更新],因为涉及到后面的参数更新方式?

XGBoost如何使用自定义的损失函数?

模型直接调用fit函数无法传入自定义的损失函数, 需要在模型开始定义的时候传入或者使用xgb.train函数调用

• 使用方法1:

  1 2 3

  from xgboost import XGBClassifier clf = XGBClassifier(objective=MyLossFunction)

• 使用方法2:

  1 2 3 4 5

  import xgboost as xgb from xgboost import XGBClassifier clf = XGBClassifier() xgb.train(xgb_model=clf, obj=MyLossFunction)

各种树的介绍

树

• 一个根节点,每个结点可能有多个子节点

二叉树

• 一个根节点,或者为空
• 每个结点只有两个子节点

二叉搜索树

也叫二叉查找树

• 首先是一棵二叉树
• 左边孩子结点的值都小于当前结点
• 右边孩子结点的值都大于当前结点

缺点

• 极端情况下,树模型会退化成链表,查找变成了 O(n) 复杂度的

平衡二叉搜索树(AVL树)

也叫平衡二叉查找树

• 首先是一棵二叉搜索树
• 对每个结点而言, 左右孩子结点的深度差值不能超过1 ， 从而保证查找是 O(log n) 的
• 控制平衡方法: 参考链接AVL树详解
  • 左-左型: 右旋
  • 右-右型: 左旋
  • 左-右型: 左旋 + 右旋
    • 第一步左旋后面部分,变成 左-左 型
    • 第二步使用右旋修正 左-左 型
  • 右-左型: 右旋 + 左旋
    • 第一步右旋后面部分,变成 右-右 型
    • 第二步使用左旋修正 右-右 型

缺点

• 每棵树的左右子树高度最多差1这个要求太严了
• 几乎每次插入或者删除结点时都会造成规则破坏, 也就需要频繁的通过左旋或者右旋操作来修正
• 插入和删除太频繁的场景中不太适合使用AVL树, 性能会因为左右子树高度最多差1这个规则频繁被打破而降低

红黑树

• 首先是一棵二叉搜索树
• 每个结点不是黑色就是红色
• 根节点为黑色
• 每个叶子结点都为黑色的空结点(NULL)
  • 注意: 叶子节点不存数据
• 任何相邻结点不同时为红色
  • 注意,相邻结点可以为黑色,且可以某条路径上全是黑色
• 对每个结点而言,当前结点到叶子结点的所有路径包含相同的黑色结点数

优点

• 能保证查找时间复杂度是 O(log n) 的, 和AVL树差不多
  • 证明: [待更新]
• 插入和删除操作中不会频繁破坏红黑树的规则

红黑树的应用

• 容器集合 HashMap, TreeMap等
  • HashMap是 链表 + 红黑树的实现, 当冲突时就需要使用红黑树加快检索
  • HashMap中 链表长度太短时使用链表, 太长时使用红黑树, 这个阈值一般设置为8

二三树

• 红黑树是二三树的一个变形,一般用红黑树就够了
  [待更新]

B树

• B树在大量的数据库和文件系统场景中都有使用
  [待更新]

B+树

[待更新]

总结

• 可以说红黑树是一种不严格的平衡树

C++和Java等语言都有++和–操作,为什么以方便自居的Python却没有这种操作呢?

Python的数值对象

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

a = 1
b = 1
c = 1000123
d = 1000123

print id(a)
print id(b)
print id(c)
print id(d)

## output:
# 94181316498840
# 94181316498840
# 94181323965720
# 94181323965720

• Python数值对象都是不可变类型,与String一样,所以不能修改对象内部数据
• C++中的i++修改内存中对象本身,数值增加1,而Python不能修改对象本身
• 与C++中字符串可以修改,Python中不能修改是一个道理

本文介绍Python中优先队列的用法

• 注意: queue并不能算是Python标准库,所以在LeetCode等OJ环境中不能使用, 想要使用优先队列的话可以使用Python的标准库heapq
• heapq的使用请参考 Python——heapq模块-最大堆最小堆

导入包

1
2
3

from Queue import PriorityQueue
# or
from queue import PriorityQueue

• queue包名已经弃用,测试发现本地Python2.7环境可以用,但是LeetCode线上环境不能用
• 推荐使用Queue

使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

from Queue import PriorityQueue

pq = PriorityQueue()

pq.put((10, "start"))
pq.put((5, "b", 12, 123))
pq.put((5, "a", 6))
pq.put(1)
pq.put(4)
pq.put([0, "a"])
pq.put([8, "b"])
pq.put("avb")
pq.put(None)

# while pq.not_empty()
# while not pq.empty()
while pq.qsize():
    print pq.get()
## output:
# None
# 1
# 4
# [0, 'a']
# [8, 'b']
# avb
# (5, 'a', 6)
# (5, 'b', 12, 123)
# (10, 'start')

• 不能使用not pq这样的语句判断优先队列是否收敛,他不是普通的内嵌对象(list,str等是内嵌对象),除非pq == None否则,双端队列对象永远为not pq == False

• 下面用法最优:

  1	while pq.qsize():

• PriorityQueue中默认递增排序(这一点与Python中的sorted函数和sort()函数一样),每次get(),移除并返回最小的对象

• PriorityQueue中,可以同时添加不同类别的对象

• PriorityQueue会将对象首先按照类别排序,然后各个类别内部按照不同数值排序

• 若传入对象是可以直接比较大小的类型即可直接传入,包括tuple, list, str, int(long)等类型

  • Python中list,tuple,str等都是可以直接比较大小的,默认使用他们的第一个元素比较大小,如果第一个元素相等,则比较第二个元素,以此类推
  • 详细情况参考本文后面的说明

Python中的内嵌对象比较大小

• Python中类别间也可以比较大小,默认类别间大小为:tuple > str > list > int(long) > None, 但是记不清楚的话不建议使用Python的这个特性,容易造成错误
• Python中list,tuple,str等都是可以直接比较大小的,默认使用他们的第一个元素比较大小,如果第一个元素相等,则比较第二个元素,以此类推
• Python中类别间也可以比较大小,默认类别间大小为:tuple > str > list > int(long) > None, 但是记不清楚的话不建议使用Python的这个特性,容易造成错误

  1 2 3 4 5

  ls = [(1, 'b'), (1, 'a'), (2, 'a'), [1, 'a'], [1, 'b'], [2, 'a'], '1a', '1b', '2a', 1, 2, 3, None] ls.sort() print ls # # output # [None, 1, 2, 3, [1, 'a'], [1, 'b'], [2, 'a'], '1a', '1b', '2a', (1, 'a'), (1, 'b'), (2, 'a')]

Python自定义对象比较大小

• 在做算法题时,没必要的情况下不建议使用

• 在做工程时建议使用这种方式

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

  import Queue class Node(): def __init__(self, val): self.val = val def __lt__(self, other): return self.val < other.val pq = Queue.PriorityQueue() pq.put(Node(5)) pq.put(Node(1)) while pq.qsize(): print pq.get().val # # output # 1 # 5

• 注意__lt__函数中是小于号,说明递增排序,大于号,说明递减排序

• PriorityQueue对象不是普通Python内嵌对象,不能使用Python内嵌的len函数

Python的list有很多强大的功能,有些比较罕见的操作可能很有用,需要我们记住

list的常见操作

list子列表

• 注意使用子列表时是一个新对象,操作子列表与原始list无关
• 在快速排序和归并排序中不可将子列表传入,以期待可以从函数中修改原始列表的值

list反序子列表

• list 反序列示例：

  1 2 3 4 5 6 7

  l = [1, 2, 3] l1 = l[::-1] print l print l1 # # output: # [1, 2, 3] # [3, 2, 1]

list的罕见操作

remove(object)

• 移除列表中第一个与object相等的对象

  1 2 3 4 5 6

  l = [1, 2, 2, 3, 4] l.remove(2) print l # # output: # [1, 2, 3, 4]

pop(index)

• 从列表中移除一个元素,并返回该元素,index为索引
• 默认移除最后一个

  1 2 3 4 5 6 7 8

  l = [1, 2, 2, 3, 4, 5] l.pop(0) print l l.pop() print l # # output: # [2, 2, 3, 4, 5] # [2, 2, 3, 4]

本文介绍Python中双端队列(double-ended queue, 简称为deque)的用法

导入包

1
2
from collections import deque
import collections

使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import collections
dq = collections.deque()
dq.append(3)
dq.append(4)
dq.append(1)
dq.appendleft(9)
dq.appendleft(10)
print dq
print dq.pop()
print dq
print dq.popleft()
print dq

while len(dq):
    print dq.pop()


# # output:
# deque([10, 9, 3, 4, 1])
# 1
# deque([10, 9, 3, 4])
# 10
# deque([9, 3, 4])
# 4
# 3
# 9

• 注意,deque没有qsize()函数,但是可以像普通队列一样使用Python内嵌的len函数

整体说明

• 打包本地项目为的 pip 包时，可以使用 pip install . 命令
• 执行 pip install . 命令时，Python 包管理工具 pip 会根据当前目录中的 setup.py 文件来安装包
  • 这个过程涉及多个步骤，包括解析包的元数据、处理依赖关系、构建和安装包等
  • 安装的内容包括包的所有模块、数据文件、编译文件以及命令行工具
  • 默认情况下，只安装核心依赖，可选依赖需要用户显式指定

构建一个包的步骤

第一步：解析和构建

• 解析 setup.py (setup.py 的详细示例见附录)：
  • pip 首先会查找并解析当前目录中的 setup.py 文件
  • 读取包的元数据（如包名、版本、作者信息等）和安装选项（如 install_requires、extras_require 等）
• 构建包：
  • pip 会使用 setuptools 或 distutils 来构建包
  • 包的构建包括编译C扩展（如果存在）、处理包内的数据文件等
  • 生成一个源代码分发包（Source Distribution，简称 sdist）或一个二进制分发包（Binary Distribution，简称 wheel）

第二步：处理依赖

• 安装核心依赖
  • pip 会根据 setup.py 中的 install_requires 列表安装包的核心依赖
  • 这些依赖是包运行所必须的
• 可选依赖
  • 如果用户在命令中指定了可选依赖组（如 pip install .[dev]），pip 会安装对应的 extras_require 可选依赖
  • 否则，默认情况下不会安装 extras_require 中的可选依赖

第三步：安装包**

• 导入包和模块
  • pip 会将构建好的包安装到Python环境的 site-packages 目录中
  • 包中的所有模块和子包都会被导入
• 包数据文件
  • 如果 setup.py 中设置了 include_package_data=True 或指定了 package_data，这些数据文件也会被安装
• 命令行工具
  • 如果 setup.py 中定义了 entry_points，对应的命令行工具会被安装并注册

第四步：打包内容

• 包和模块
  • 所有通过 packages 或 find_packages() 指定的包和模块
• 数据文件
  • 通过 package_data 或 include_package_data 指定的额外数据文件
• 编译文件
  • 如果包中包含 C/C++ 扩展模块，这些模块会被编译并打包

附录：setup.py 文件的示例

• setup.py 文件是 Python 项目中用于定义包的元数据、依赖关系及其他安装相关信息的脚本
• setup.py 文件主要通过 setuptools 或 distutils 库来实现包的配置和分发

setup.py 文件示例

• setup.py 文件基本结构示例：

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

  from setuptools import setup, find_packages setup( name="package_name", # 包名 version="0.1.0", # 版本号 author="Author Name", # 作者 author_email="author@example.com", # 作者邮箱 description="A short description", # 简短描述 long_description=open("README.md").read(), # 长描述（通常从README文件导入） long_description_content_type="text/markdown", # 长描述类型（如Markdown） url="https://github.com/username/repo", # 项目URL packages=find_packages(), # 自动发现并包含所有包 classifiers=[ # 分类标签（如开发状态、受众、许可证） "Programming Language :: Python :: 3", "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", ], python_requires='>=3.6', # Python版本要求 install_requires=[ # 安装依赖 "requests>=2.25.0", "numpy", ], extras_require={ # 可选依赖（如开发、测试） "dev": ["pytest", "sphinx"], }, entry_points={ # 命令行工具入口 "console_scripts": [ "mycommand=mypackage.module:function", ], }, include_package_data=True, # 包含包内数据文件 package_data={ # 指定包数据文件 "mypackage": ["data/*.dat"], }, zip_safe=False, # 是否允许打包为zip文件 )

主要内容和功能说明

• 元数据：
  • name ： 包的名称
  • version ： 包的版本号，通常遵循语义版本规范（如 MAJOR.MINOR.PATCH）
  • author 和 author_email ： 作者的姓名和联系邮箱
  • description 和 long_description ： 包的简短和详细描述
  • url ： 项目的主页或代码仓库地址
• 包和模块：
  • packages ： 指定需要包含的包列表，通常使用 find_packages() 自动发现
  • package_data ： 指定需要包含的非代码文件（如数据文件、模板等）
• 依赖管理：
  • install_requires ： 安装包时需要满足的依赖列表
  • extras_require ： 可选依赖，按功能分组（如开发依赖、测试依赖）
• Python 版本：
  • python_requires ： 指定包支持的Python版本范围
• 分类标签：
  • classifiers ： 用于描述包的特性和分类，帮助用户和工具识别包的适用性
• 命令行工具：
  • entry_points ： 定义包提供的命令行工具及其入口函数
• 数据文件：
  • include_package_data ： 是否包含包内的额外数据文件
• 打包选项：
  • zip_safe ： 指定包是否可以安全地打包为zip文件

附录：setup.py 中的 extras_require 使用

• 在 setup.py 中定义的 extras_require 提供了一种机制，可以按需安装额外的依赖项

• extras_require 允许定义一些可选的依赖组，比如开发依赖、测试依赖等

• 举例来说，前面附录小节的示例中 extras_require 中定义了一个名为 “dev” 的组和对应的依赖

• 默认情况下 ，安装包时不会自动安装 extras_require 中定义的可选依赖项

  • 直接使用 pip install . 或 pip install package_name 安装包即可
  • 在这种情况下，只有 install_requires 中定义的核心依赖会被安装

• 安装可选依赖，以 安装 “dev” 组的可选依赖 为例：

  • 使用 pip 安装时，可以通过在包名后加上 [dev] 来指定安装这些可选依赖

    1	pip install .[dev]

  • 如果包已经发布到PyPI或其他包管理平台，可以使用以下命令：

    1	pip install package_name[dev]

本文主要区分隐马尔可夫模型(HMM),最大熵马尔可夫模型(MEMM),和条件随机场(CRF)

HMM

• 生成式模型
• 有向图模型,贝叶斯网络

模型描述

• 模型参数: \(\lambda = (A, B, \pi)\)
  • A为状态转移矩阵
  • B为观测概率矩阵
  • \(\pi\) 为初始状态概率向量
• HMM对

假设

• 观测序列之间独立
• 当前状态仅仅依赖于上一个状态

问题

• 概率计算问题: 给定模型 \(\lambda = (A, B, \pi)\) 和观测序列 \(O=(o_{1}, o_{2},,,o_{n})\)，计算在给定模型下观测序列出现的概率 \(P(O|\lambda)\)
• 学习问题: 已知观测序列 \(O=(o_{1}, o_{2},,,o_{n})\)，估计模型参数 \(\lambda = (A, B, \pi)\)，使得在该模型下观测序列出现的概率 \(P(O|\lambda)\) 最大.(极大似然法,EM算法)
• 预测问题(解码问题): 给定模型 \(\lambda = (A, B, \pi)\) 和观测序列 \(O=(o_{1}, o_{2},,,o_{n})\)，求状态序列 \(I=(i_{1}, i_{2},,,i_{n})\)，使得 \(P(I|O;\lambda)\) 最大,(维特比算法)

序列标注问题

• 序列标注问题是已知观测序列 \(O=(o_{1}, o_{2},,,o_{n})\)，求状态序列 \(I=(i_{1}, i_{2},,,i_{n})\)，使得 \(P(I|O)\) 最大
• 实际上序列标注问题包括学习问题和预测问题两个问题:
  • 学习问题: 根据观测序列确定模型参数 \(\lambda = (A, B, \pi)\)，极大似然法或EM算法(EM算法会同时估计得到最优状态序列(隐变量))
  • 预测问题: 根据模型参数和观测序列确定最优状态序列 \(I=(i_{1}, i_{2},,,i_{n})\)，维特比算法

优点

• 算法成熟
• 效率高, 模型简单, 容易训练

缺点

• 序列标注问题中,当前状态(标注)往往不仅仅和前一个状态相关,还可能和观察序列相关(这里指的是整个序列)
• 也就是说每个状态可能还与整个观察序列的(除了当前观察值以外的)其他观察值(观察序列上下文)相关

MEMM

• 判别式模型
• 有向图模型,贝叶斯网络

假设

• 当前状态仅依赖于上一状态和当前观测值(或所有观测值)
• (问题: 为什么有些书上画出的图是当前状态依赖上一状态和所有观测值?,这里应该是当前观测值和所有观测值两种情况都是MEMM,<<百面>>画的是所有观测值的情况)

问题

• MEMM似乎只用于序列标注,也就是在已知观测序列的情况下,寻找最优的状态序列
• 有其他应用的话再添加

序列标注问题

• 用于序列标注时,一般也包括两个问题: 学习问题和预测问题
  • 学习问题: 根据观测序列确定模型参数(每条边的(概率)值和初始状态?)
  • 预测问题: 根据模型和观测序列确定维特比算法

优点

• 解决了观测独立性问题(观测独立性是只当前观测序列只与当前状态相关)[问题: 在MEMM中并不关心观测序列由谁影响,而是关心观测序列如何影响了状态序列]

缺点

• 标签偏置(labeling bias)问题
  • MEMM中概率最大路径往往容易出现在转移少的状态中
  • MEMM归一化在加和函数 \(\sum\) 计算内部,而CRF的归一化在加和函数 \(\sum\) 的外部,这使得MEMM只会关注加和函数[原始建模问题概率值 \((y_{1…n}|x_{1…n})\) ]的局部特征,而不是的整体特征,所以MEMM存在偏置问题
• 比HMM复杂

CRF

• 判别式模型
• 无向图模型,马尔可夫网络

假设

• 当前状态仅依赖于上一状态和当前观测值(或所有观测值)
• (问题: 为什么有些书上画出的图是当前状态依赖上一状态和所有观测值?,这里应该是当前观测值和所有观测值两种情况都是线性CRFs,<<百面>>画的是所有观测值的情况)
• 与MEMM的区别就是无向图模型与有向图模型的区别

问题

序列标注问题

优点

• 模型复杂,能建模更多可能的特征
• 全局归一化(这里与MEMM的区别是,MEMM归一化在加和函数[原始建模问题概率值 \((y_{1…n}|x_{1…n})\) ] \(\sum\) 计算内部,而CRF的归一化在加和函数[原始建模问题概率值 \((y_{1…n}|x_{1…n})\) ] \(\sum\) 的外部)

缺点

• 模型复杂,速度慢

期望最大化(Exception Maximization Algorithm)EM算法

不同教材的不同形式

李航统计学习方法

算法步骤

• 输入: 观测变量数据Y

• 输出: 模型(参数 \(\theta\))

• E步: 计算 \(Q(\theta, \theta^{i})\)
  $$
  \begin{align}
  Q(\theta, \theta^{i}) &= \mathbb{E}_{Z}[logP(Y,Z|\theta)|Y,\theta^{i}] \\
  &= \mathbb{E}_{Z\sim P(Z|Y,\theta^{i})}[logP(Y,Z|\theta)] \\
  &= \sum_{Z} P(Z|Y,\theta^{i})logP(Y,Z|\theta) \\
  &= \sum_{Z} logP(Y,Z|\theta)P(Z|Y,\theta^{i})
  \end{align}
  $$

• M步: 求使得 \(Q(\theta, \theta^{i})\) 极大化的参数 \(\theta=\theta^{i+1}\)
  $$\theta^{i+1} = \mathop{\arg\max}_{\theta}Q(\theta, \theta^{i})$$

• 重复E步和M步,直到收敛

• 理解: \(Q(\theta, \theta^{i})\) 可以理解为 \(Q(\theta|\theta^{i})\)，表示在参数 \(\theta^{i}\) 已知的情况下,对数似然函数关于隐变量后验分布的期望函数,函数的参数为 \(\theta\)

隐变量的期望还是分布?

参考博客:https://www.jianshu.com/p/c3ff1ae5cb66

仅考虑隐变量的期望

• 应用场景为k-means聚类,但是k-means聚类E步求的是最可能的 \(Z\) 值(概率最大的 \(Z\) 值),而不是 \(Z\) 的期望

考虑隐变量的分布

• 应用场景为GMM模型聚类

推导

• 已知数据是观测数据Y,像极大似然法一样,使得似然函数最大化即可,这里为了方便计算使用对数似然

• 我们的终极目标与极大似然法一样,求一个使得似然函数最大(可能是极大)的参数$$\theta^{\star}=\mathop{\arg\max}_{\theta}L(\theta)$$
  其中
  $$
  \begin{align}
  L(\theta)&=logP(Y|\theta)\\
  &=log\sum_{Z}P(Y,Z|\theta)\\
  &=log\left(\sum_{Z}P(Y|Z,\theta)P(Z|\theta)\right)
  \end{align}
  $$

• 显然,上述似然函数比较难以求解,非凸,且涉及和(积分或加法)的对数等操作,难以展开(可能还有其他的原因)

• 所以我们使用EM算法迭代不断逼近原始似然函数的最优解 \(\theta^{\star}\) (注意:我们只能得到局部最优,但是一般来说局部最优也够用了)

• 可用Jensen不等式得到 \(L(\theta)\) 的下界来作为优化目标不断迭代
  $$
  \begin{align}
  L(\theta)&=log\left(\sum_{Z}P(Y|Z,\theta)P(Z|\theta)\right)\\
  &=log\left(\sum_{Z}P(Z|Y,\theta^{i})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{i})}\right)\\
  &\geq \sum_{Z}P(Z|Y,\theta^{i})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{i})}\\
  &=B(\theta,\theta^{i})
  \end{align}
  $$

• 由于此时 \(B(\theta, \theta^{i})\) 是 \(L(\theta)\) 的下界(\(B(\theta, \theta^{i})\) 是固定 \(\theta^{i}\) 时关于 \(\theta\) 的凸函数且容易求导,可以求极大值),所以使得前者增大的参数 \(\theta\) 也能使得后者增大,为了使得后者尽可能的增大,我们对前者取极大值

  • 下面的推导基于事实: 消去与 \(\theta\) 无关的项,极大值点(\(\theta^{i+1}\))不变

$$
\begin{align}
\theta^{i+1} &= \mathop{\arg\max}_{\theta}B(\theta,\theta^{i})\\
&=\mathop{\arg\max}_{\theta}\left( \sum_{Z}P(Z|Y,\theta^{i})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{i})}\right)\\
&=\mathop{\arg\max}_{\theta}\left( \sum_{Z}P(Z|Y,\theta^{i})logP(Y|Z,\theta)P(Z|\theta)-\sum_{Z}P(Z|Y,\theta^{i})logP(Z|Y,\theta^{i})\right)\\
&= \mathop{\arg\max}_{\theta}\left( \sum_{Z}P(Z|Y,\theta^{i})logP(Y|Z,\theta)P(Z|\theta)\right)\\
&=\mathop{\arg\max}_{\theta}Q(\theta,\theta^{i})
\end{align}
$$

• 由上式可知,我们的EM算法步骤中E步和M步是正确的
• 问题:推导第二步中为什么选择 \(P(Z|Y,\theta^{i})\) 而不是其他分布呢?
  • 解答: \(B(\theta, \theta^{i})\) 和 \(L(\theta)\) 什么时候相等呢?(前面推导中Jensen不等式什么时候能取等号呢?)
  • Jensen不等式取等号当且仅当Jensen不等式中函数的值为常数,此处函数的值为 \(log\frac{P(Y,Z|\theta)}{Q(Z)}\)
    $$
    \begin{align}
    L(\theta)&=log\left(\sum_{Z}P(Y,Z|\theta)\right)\\
    &=log\left(\sum_{Z}Q(Z)\frac{P(Y,Z|\theta)}{Q(Z)}\right)\\
    &\geq \sum_{Z}Q(Z)log\frac{P(Y,Z|\theta)}{Q(Z)}\\
    \end{align}
    $$
  • 不等式中当且仅当 \(log\frac{P(Y,Z|\theta)}{Q(Z)}\) 为常数,也就是 \(\frac{P(Y,Z|\theta)}{Q(Z)}=c\)，c为常数,时等号成立
  • 此时由于 \(Q(Z)\) 是一个分布(注意:正因为 \(Q(Z)\) 是一个分布才能用Jensen不等式),所以有
    $$
    \begin{align}
    Q(Z)=\frac{P(Y,Z|\theta)}{\sum_{Z}P(Y,Z|\theta)}=\frac{P(Y,Z|\theta)}{P(Y|\theta)}=P(Z|Y,\theta)
    \end{align}
    $$

吴恩达CS229

• E步: 计算 \(Q_{i}(Z)=P(Z|Y,\theta^{i})\)
• M步: 求使得原始似然函数下界极大化的参数 \(\theta=\theta^{i+1}\)
  $$
  \begin{align}
  \theta^{i+1} &= \mathop{\arg\max}_{\theta}\sum_{Z}P(Z|Y,\theta^{i})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{i})} \\
  &= \mathop{\arg\max}_{\theta}\sum_{Z}Q_{i}(Z)log\frac{P(Y|Z,\theta)P(Z|\theta)}{Q_{i}(Z)}
  \end{align}
  $$
  • 进一步消除与 \(\theta\) 无关的项可以得到
    $$
    \begin{align}
    \theta^{i+1} &= \mathop{\arg\max}_{\theta}\sum_{Z}Q_{i}(Z)logP(Y|Z,\theta)P(Z|\theta)
    \end{align}
    $$
• 推导步骤和李航统计学习方法一样,核心是运用Jensen不等式

总结

• 以上两个不同课程的E步不同,但完全等价,吴恩达CS229课程中E步计算 \(Q_{i}(Z)=P(Z|Y,\theta^{i})\) 就等价于计算出了李航统计学习方法中的 \(Q(\theta, \theta^{i})\)，二者关系如下:
  $$
  \begin{align}
  Q(\theta, \theta^{i})&=\sum_{Z}P(Z|Y,\theta^{i})logP(Y|Z,\theta)P(Z|\theta) \\
  &= \sum_{Z}Q_{i}(Z)logP(Y|Z,\theta)P(Z|\theta)
  \end{align}
  $$
• M步中,二者本质上完全相同,但是吴恩达CS229中没消去与 \(\theta\) 无关的项,所以看起来不太简洁

实例

实例一

• 三个硬币的抛硬币问题
  • 第一次:抛硬币A,决定第二次抛C还是B,选中B的概率为 \(\pi\)
  • 第二次:抛硬币B或C,正面为1,反面为0
• (第一个硬币A抛中B的概率为隐变量)

实例二

• 200个人不知道男女的身高拟合问题(性别为隐变量)
• 这里是高斯混合模型的代表

实例三

• K-Means聚类
• 参考<<百面机器学习>>P100-P101

进一步理解

我的理解

• 初始化参数 \(\theta^{0}\)

• 1.根据参数 \(\theta^{i}\) 计算当前隐变量的分布函数 \(Q_{i}(Z)=P(Z|Y,\theta^{i})\)

  • 这一步的本质是使得在参数 \(\theta = \theta^i\) 时, 求得一个隐变量 \(Z\) 的分布,使得原始式子中的不等式取等号

• 2.根据 \(Q_{i}(Z)=P(Z|Y,\theta^{i})\) 得到对数似然函数下界函数 \(B(\theta, \theta^{i})\) (求原始似然函数的下界B函数是因为直接对原始似然函数求极大值很难)或者 \(Q(\theta,\theta^{i})\)
  $$\mathop{\arg\max}_{\theta}B(\theta, \theta^{i})=\mathop{\arg\max}_{\theta}Q(\theta,\theta^{i})$$
  (\(Q(\theta,\theta^{i})\) 可以看作是 \(B(\theta, \theta^{i})\) 的简化版, \(B(\theta, \theta^{i})\) 才是原始似然函数的下界, \(Q(\theta,\theta^{i})\) 不是原始似然函数的下界)

  • 这里 \(B(\theta, \theta^{i})\) 就是原始似然函数的下界(也就是不等式取到等号)

• 3.求使得函数 \(B(\theta, \theta^{i})\) 极大化的参数 \(\theta=\theta^{i+1}\)

  • 这一步是在固定隐变量 \(Z\) 的分布时, 用极大似然求一个使得下界 \(B(\theta, \theta^{i})\) 最大的参数 \(\theta = \theta^{i+1}\) 使得 \(B(\theta^{i+1}, \theta^{i}) = \text{max} B(\theta, \theta^{i})\)

• 4.循环1,2,3直到收敛(相邻两次参数的变化或者是似然函数的变化足够小即可判断为收敛)

  • \(||\theta^{i+1}-\theta^{i}||<\epsilon_{1}\) 或者 \(||Q(\theta^{i+1},\theta^{i})-Q(\theta^{i},\theta^{i})||<\epsilon_{2}\)

• 总结:

  • 在吴恩达CS229课程中: E步包含1, M步包含2,3,其中第2步中求的是 \(B(\theta, \theta^{i})\)
  • 在李航<<统计学习方法>>中: E步包含1,2, M步包含3,其中第2步中求的是 \(Q(\theta,\theta^{i})\)
  • 两种表达等价

图示理解

• 迭代图示如下图(图来自博客:https://www.cnblogs.com/xieyue/p/4384915.html)

• 也可参考李航<<统计学习方法>>第160页的图和解释

收敛性

• 参考李航<<统计学习方法>>第160页-第162页推导过程
• 参考<<百面机器学习>>第P099页-第P100页
  • 核心思想原始函数单调有界
  • 原始函数为 \(L(\theta)\)，函数下界为 \(B(\theta,\theta^{i})\)
  • E步:
    • 找到使得在当前 \(\theta^{i}\) 确定时,原始函数的下界 \(B(\theta, \theta^{i})\)，在 \(\theta^{i}\) 处有
      \(\)
      $$
      \begin{align}
      L(\theta^{i}) = B(\theta,\theta^{i})
      \end{align}
      $$
  • M步:
    • 找到使得函数 \(B(\theta,\theta^{i})\) 取得极大值的 \(\theta^{i+1}\)
  • i = i + 1,然后重新开始E和M步
    $$
    \begin{align}
    L(\theta^{i+1}) >= L(\theta^{i+1})
    \end{align}
    $$
  • 所以函数是单调的
  • 由于 \(L(\theta)\) 有界(这里原始函数有界可以从似然函数的定义得到)
  • 函数单调有界=>函数收敛(数学分析中的定理)

优劣性

优势

• 简单性
• 普适性

劣势

• 不能保证收敛到最大值,只能保证收敛到极大值
• 对初始值敏感,不同初始值可能收敛到不同极值点
• 实际使用时通常采用多次选择不同的初始值来进行迭代,最终对估计值选最好的

EM算法的推广

• 引入F函数

GEM1

F函数极大-极大法

• 初始化参数 \(\theta^{0}\)
• E步: 求使得 \(F(\tilde{P},\theta^{i})\) 极大化的 \(\tilde{P}^{i+1}\)
• M步: 求使得 \(F(\tilde{P}^{i+1},\theta)\) 极大化的 \(\theta^{i+1}\)
• 重复E,M,直到收敛

GEM2

• 初始化参数 \(\theta^{0}\)
• E步: 计算 \(Q(\theta,\theta^{i})\)
• M步: 求 \(\theta^{i+1}\) 使得 \(Q(\theta^{i+1},\theta^{i}) > Q(\theta^{i},\theta^{i})\)
• 重复E,M,直到收敛
• 总结: \(Q(\theta,\theta^{i})\) 的极大化难求时,这种方法可以简化计算

GEM3

• 初始化参数 \(\theta^{0}\)
• E步: 计算 \(Q(\theta,\theta^{i})\)
• M步: 对参数 \(\theta^{i}\) 的每一个维度k,固定参数的其他维度,求使得 \(Q(\theta,\theta^{i})\) 极大化的 \(\theta_{k}^{i+1}\)，最终得到 \(\theta^{i+1}\)
  • 使得 \(Q(\theta^{i+1},\theta^{i}) > Q(\theta^{i},\theta^{i})\)
• 重复E,M,直到收敛
• 总结: 一种特殊的GEM算法,将M步分解为参数 \(\theta\) 的维度次来条件极大化

本文介绍MCMC采样法和他的两个常用方法:MH(Metropolis-Hastings)采样法和Gibbs采样法
马尔可夫蒙特卡罗(Markov Chain Monte Carlo, MCMC)采样法

• 对于高维空间中的随机向量,拒绝采样和重要性采样经常难以找到合适的参考分布,容易导致采样效率低下(样本的接受概率太小或者重要性权重太低),此时可以考虑马尔可夫蒙特卡罗采样法(MCMC)
• MCMC中常见的有两种,MH(Metropolis-Hastings)采样法和Gibbs采样法

MCMC概述

• 马尔可夫蒙特卡罗(Markov Chain Monte Carlo, MCMC)采样法可分为两个部分(两个MC)描述

蒙特卡罗法

• 蒙特卡罗法(Monte Carlo)是指: 基于采样的数值型近似求解方法

马尔可夫链

又称离散时间马尔可夫链(discrete-time Markov chain，缩写为DTMC),或者马氏链

• 马尔可夫链(Markov Chain)是指: 状态空间中经过从一个状态到另一个状态的转换的随机过程, 该随机过程满足马尔可夫性
• 马尔可夫性(Markov property)是指: 当一个随机过程在给定现在状态及所有过去状态情况下，其未来状态的条件概率分布仅依赖于当前状态；换句话说，在给定现在状态时，它与过去状态（即该过程的历史路径）是条件独立的，那么此随机过程即具有马尔可夫性质。具有马尔可夫性质的过程通常称之为马尔可夫过程
  • 马尔可夫性的简单理解: “无记忆性”: 下一状态的概率分布只能由当前状态决定，在时间序列中它前面的事件均与之无关

MCMC基本思想

• 针对采样的目标分布,构造一个马尔可夫链 ，使得该马尔可夫链的平稳分布就是目标分布
• 从任何一个初始状态出发,沿着马尔可夫链进行状态转移,直到马尔可夫链收敛(到达平稳状态)
• 继续在马尔可夫链上进行状态转移,收敛后继续采样得到的样本就是原始目标分布的样本
  • burn-in处理 : 现实应用中,我们需要丢弃收敛前的样本,只保留收敛后的样本
  • burn-in 原意为”老化,定型”之意,在这里表示我们只取后面马氏链定型(收敛)后采样得到的样本
  • 假设采样到收敛用了n次采样,那么服从原始分布的k个样本为 \((x^{n+1}, x^{n+2},,,, x^{n+k})\) 有时候为了得到近似独立的样本,可以间隔每r次再取出其中一个样本 \((x^{n+r}, x^{n+2r},,,, x^{n+kr})\)
  • 真正独立同分布的k个样本可用多条k条不同的收敛后的马氏链得到,不同马氏链采样得到的样本是独立同分布的
• 核心: 马尔可夫链的构造 ，也就是确定马尔可夫链的状态转移概率

常见的MCMC方法

MH(Metropolis-Hastings)采样法

• 对于原始目标分布 \(p(x)\)
• 选择一个参考条件分布 \(q(x^{\star}|x)\)，定义接受概率 \(A(x,x^{\star})\)：
  • (注意这里是参考条件分布,因为是马尔可夫链,所以每个状态都由上一个状态转移而来,需要定义的参考分布应该是条件分布,不是一般拒绝采样中的普通参考分布)
    $$
    A(x,x^{\star}) = min\left ( 1, \frac{p(x^{\star})q(x|x^{\star})}{p(x)q(x^{\star}|x)} \right )
    $$
• MH采样法构建满足平稳分布就是目标分布 \(p(x)\) 的秘诀就是让每次采样时,当前状态以一定概率停留在上一状态
  • 与拒绝采样对应: 接受意味着跳转到下一状态,拒绝意味着停留在当前状态

采样过程

• 随机选取初始样本x^{0}
• for t = 1,2,3,…:
  • 参考条件分布采样 \(x^{\star} \sim q(x^{star}|x^{t-1})\)
  • 均匀分布采样 \(u \sim U(0,1)\)
  • 判断是否接受: 如果 \(u < A(x^{t-1}, x^{\star})\)，则接受,令: \(x^{t} = x^{\star}\)，否则拒绝,令: \(x^{t}=x^{t-1}\)
• burn-in处理 : 丢弃采样到平稳分布前的样本, 只保留平稳分布后的样本即为服从原始分布 \(p(x)\) 的样本
  • 假设采样到收敛用了n次采样,那么服从原始分布的k个样本为 \((x^{n+1}, x^{n+2},,,, x^{n+k})\) 有时候为了得到近似独立的样本,可以间隔每r次再取出其中一个样本 \((x^{n+r}, x^{n+2r},,,, x^{n+kr})\)
  • 真正独立同分布的k个样本可用多条k条不同的收敛后的马氏链得到,不同马氏链采样得到的样本是独立同分布的
• 采样次数一般来说是凭经验选择一个足够大的值,现实是现实可以使用一些参数变化量这类的指标来判断采样是否收敛,参考NLP——LLDA的Gibbs采样实现

与拒绝采样的区别

• MH采样基于拒绝采样来逼近平稳分布
• 拒绝采样中: 如果样本某一步被拒绝,那么该步不会产生新的样本,需要重新对当前步进行采样
• MH中: 每一步都会产生一个样本,被拒绝后,就令当前样本和上一个样本相同即可
  • 因为这里是为了使得每个状态的转出概率等于转入概率,所以拒绝就意味着当前采样步骤状态不跳转
  • MH采样法最核心的思想就是一定概率停留在上一个状态来实现对马尔可夫链的构建的

MH采样法正确性证明

• MH采样法构造的马尔可夫链(状态转移概率矩阵)是正确的吗?

• 细致平稳条件 ， 如果非周期的马氏链的状态转移矩阵P和分布 \(\pi(x)\) 满足下面的式子对任意的 \(i,j\) 都成立:
  $$\pi(x^{i})P_{ij} = \pi(x^{j})P_{ji}$$

  • 上式为细致平稳分布条件(detailed balance condition)
  • 其中 \(\pi(x)\) 为马氏链的平稳分布,在这里等于我们的原始分布 \(p(x)\)

• 证明 \(\pi(x)\) 为马氏链的平稳分布:
  $$
  \begin{align}
  \sum_{i=1}^{\infty}\pi(x^{i})P_{ij} = \sum_{i=1}^{\infty}\pi(x^{j})P_{ji} = \pi(x^{j})\sum_{i=1}^{\infty}P_{ji} = \pi(x^{j}) \\
  => \pi(x) P = \pi(x)
  \end{align}
  $$

  • 由于 \(\pi(x)\) 为方程 \(\pi(x) P = \pi(x)\) 的解,所以 \(\pi(x)\) 是状态转移矩阵P对应的马氏链的平稳分布

• 在MH采样法中

  • 参考条件分布函数本对应状态转移矩阵的一个元素, \(q(x^{i}|x^{j}) = P_{ij}\) (注意: 实际上一般不相等)
  • 但是很难构造这样方便采样的函数,于是我们使用一个接受率来修正 \(q(x^{i}|x^{j})\alpha(x^{j}, x^{i}) = P_{ij}\)
    • \(\alpha(x^{j}, x^{i})\) 表示从 \(x^{j}\) 跳转到 \(x^{i}\) 的接受率, 其值可如下求得:
      $$
      \begin{align}
      \pi(x^{i}) P_{ij} &= \pi(x^{j})P_{ji} \\
      \pi(x^{i}) q(x^{j}|x^{i})\alpha(x^{i}, x^{j}) &= \pi(x^{j})q(x^{i}|x^{j})\alpha(x^{j}, x^{i})
      \end{align}
      $$
  • 显然,直接取:
    $$
    \begin{align}
    \alpha(x^{i}, x^{j}) &= \pi(x^{j})q(x^{i}|x^{j})\\
    \alpha(x^{j}, x^{i}) &= \pi(x^{i})q(x^{j}|x^{i})\\
    \end{align}
    $$
  • 即可

• 但是由于 \(\alpha(x^{j}, x^{i})\) 一般来说太小,所以我们考虑将 \(\alpha(x^{j}, x^{i})\) 和 \(\alpha(x^{i}, x^{j})\) 同时扩大M倍,使得其中大的那个为1,即可得到最大的接受率

  • 使用原始接受率的方法称为一般MCMC方法
  • 使用扩大M被接受率的方法称为MCMC的改进方法: MH方法

• 改进后的接受率为从 \(x^{i}\) 跳转到 \(x^{j}\) 的接受率:
  $$
  \begin{align}
  A(x^{i}, x^{j}) = min\left ( 1, \frac{p(x^{j})q(x^{i}|x^{j})}{p(x^{i})q(x^{j}|x^{i})} \right )
  \end{align}
  $$

  • 理解:
    $$
    \begin{align}
    p(x^{j})q(x^{i}|x^{j}) &> p(x^{i})q(x^{j}|x^{i})时: A(x^{i}, x^{j}) = 1 \\
    p(x^{j})q(x^{i}|x^{j}) &< p(x^{i})q(x^{j}|x^{i})时: A(x^{i}, x^{j}) = \frac{p(x^{j})q(x^{i}|x^{j})}{p(x^{i})q(x^{j}|x^{i})}
    \end{align}
    $$

• 在MH中表现为从 \(x\) 跳转到 \(x^{\star}\) 的接受率:
  $$
  A(x,x^{\star}) = min\left ( 1, \frac{p(x^{\star})q(x|x^{\star})}{p(x)q(x^{\star}|x)} \right )
  $$

Gibbs采样法

Gibbs采样是MH采样法的一个特例,每次只更新样本的一个维度

• 针对维度很高的多维向量,同时采样多个维度难度较高,且接受率很小
• 使用Gibbs采样每次采样一个维度可以解决这个问题

准备

• 求得已知其他维度下,每一维度的条件概率 \(p(x_{i}|x_{1},,,x_{i-1}, x_{i+1},,,x_{d})\)

采样过程

• 随机选择初始状态 \(x^{0} = (x_{1}^{0}, x_{2}^{0}, x_{3}^{0},,,, x_{d}^{0})\)

• for t = 1,2,3,…:

  • 基于前一次产生的样本: \(x^{t-1} = (x_{1}^{t-1}, x_{2}^{t-1}, x_{3}^{t-1},,,, x_{d}^{t-1})\)，依次采样和更新每个维度的值,即依次按照:
    • \(x_{1}^{t} \sim p(x_{1}|x_{2}^{t-1},,,x_{d}^{t-1})\)
    • \(x_{2}^{t} \sim p(x_{2}|x_{1}^{t}, x_{3}^{t-1},,,x_{d}^{t-1})\)
    • \(x_{i}^{t} \sim p(x_{i}|x_{1}^{t},,,x_{i-1}^{t}, x_{i+1}^{t-1},,,x_{d}^{t-1})\)
    • \(x_{d}^{t} \sim p(x_{d}|x_{1}^{t}, x_{2}^{t},,,x_{d-1}^{t})\)
  • 得到新的一个样本: \(x^{t} = (x_{1}^{t}, x_{2}^{t}, x_{3}^{t},,,, x_{d}^{t})\)

• burn-in处理 : 丢弃采样到平稳分布前的样本, 只保留平稳分布后的样本即为服从原始分布 \(p(x)\) 的样本

  • 假设采样到收敛用了n次采样,那么服从原始分布的k个样本为 \((x^{n+1}, x^{n+2},,,, x^{n+k})\)，有时候为了得到近似独立的样本,可以间隔每r次再取出其中一个样本 \((x^{n+r}, x^{n+2r},,,, x^{n+kr})\)
  • 真正独立同分布的k个样本可用多条k条不同的收敛后的马氏链得到,不同马氏链采样得到的样本是独立同分布的
  • 注意: 采样完成部分维度生成的中间样本如 \((x_{1}^{t},,,x_{i-1}^{t}, x_{i}^{t}, x_{i+1}^{t-1},,,x_{d}^{t-1})\) 是不能作为最终样本的,因为他们之间(同一轮次的多个中间样本)相互依赖性太强,不具有独立同分布的(虽然完全采样完成的也不能视为具有独立同分布,但是可近似的认为是独立同分布的)

• 采样次数一般来说是凭经验选择一个足够大的值,现实是现实可以使用一些参数变化量这类的指标来判断采样是否收敛,参考NLP——LLDA的Gibbs采样实现