Machine Learning-贝叶斯规则分类器(上)

贝叶斯规则分类器

一、简介

在医学、证券、电信、保险等领域，对预测模型的精确度要求很高的同时，也要求模型具有很高的可解释性，目前主流的机器学习模型虽然在稳定性和精确度上表现较好，但是存在一个比较严重的问题，就是可解释性很差，本文介绍一种准确率很高且可以高度解释的模型————贝叶斯规则分类器模型（BRL）。此预测模型在医学的中风预测中，表现完胜现行的CHADS模型和其他机器学习模型，且具有高度的可解释性。

对于决策树、决策列表、关联分类器模型，使用的是贪婪算法减少计算量，但是贪婪算法严重影响模型的准确性和可解释性，如果不使用贪婪算法，在全样本空间进行搜索，计算量又是一个问题。贝叶斯规则分类器模型找到了两者之间的平衡,因为他不是以一种贪婪的、涉及分裂和剪枝的方式来构建的。贝叶斯规则分类器是使用了预先定义的规则，将模型空间减少到规则的排列，而不是所有的分割集。

此模型涉及到几大知识点，关联分析、Beta分布的共轭分布是多项式分布、MCMC采样，下面先介绍这三块的一些概念和知识。

二、模型相关基础知识

2.1关联分析

关联分析是一种简单、实用的分析技术，就是发现存在于大量数据集中的关联性或相关性，从而描述了一个事物中某些属性同时出现的规律和模式。其是从大量数据中发现项集之间有趣的关联和相关联系。关联分析的一个典型例子是购物篮分析。该过程通过发现顾客放入其购物篮中的不同商品之间的联系，分析顾客的购买习惯。通过了解哪些商品频繁地被顾客同时购买，这种关联的发现可以帮助零售商制定营销策略。其他的应用还包括价目表设计、商品促销、商品的排放和基于购买模式的顾客划分。

为了讲清关联分析，我们引入以下概念。

事物：

每一个交易记录被称为一个事物，在购物篮案例中，每一个购物的交易记录就是一个事物。

事物库：

一些事物的集合被称为事物件库。

项：

事物中的每一个元素称为一个项，在购物篮案例中，每一个商品就是一个项。

项集：

项组成的集合称为项集。

项集的支持度计数：

一个项集在事物库中出现的次数，叫做这个项集的支持度计数。

项集的支持度：

一个项集的支持度计数除以事物库中事物的总数，就是这个项集的支持度。

频繁项集：

如果我们对项目集的支持度设定一个最小阈值，那么所有支持度大于这个阈值的项集就是频繁项集。

关联规则：

项集和项集之间可以定义一个规则，如项集X，Y，可以定义关联规则：X->Y.

关联规则的支持度：

对于关联规则X->Y ,X的项和Y的项构成的新的项集的支持度为关联规则的支持度。公式如下：
$s(X->Y)=\frac{\sigma(X \bigcup Y )}{N}$
其中N为事物库中事物的总数， $\sigma(X \bigcup Y )$ 为X的项和Y的项构成的项集。

关联规则的置信度：

项集Y在包含项集X的事物中出现的频繁程度定义为关联规则的置信度。表达式如下：
$c(X->Y)=\frac{\sigma(X \bigcup Y )}{\sigma(X )}$

支持度很低的规则只能偶然出现，支持度通常用来删除那些无意义的规则。还具有一种期望的性质，可以用于关联规则的发现。
置信度度量通过规则进行推理具有可靠性。对于给定的规则，置信度越高，Y在包含X的事务中出现的可能性越大。置信度也可以估计Y在给定X的条件下概率。

对于关联规则定义这两个度量很有意义的。首先，通过对规则支持度支持度的限定滤去没有意义的规则。我们从商家的角度出发，数据挖掘意义是通过挖掘做出相应的战略决策产生价值。如果一个规则支持度很低，说明顾客同时购买这些商品的次数很少，商家针对这个规则做决策几乎没有意义。其次，置信度越大说明这个规则越可靠。

至此，我们可以通过定义上面的概念，来做关联分析，一个简单粗暴的算法是，遍历事物库中所有的关联规则，选择较大支持度和较大置信度的关联规则做关联分析。但是这个做法虽然简单，但是计算量极其大，所以引入以下两个算法：Apriori算法、FP-growth算法。

Apriori算法

这一算法思想很简单，就是为了减少暴力搜索的计算量，我们不在全部项集规则的排列组合中遍历，使用关联规则的支持度和置信度进行剪枝，使搜索范围收缩到比较小的空间，具体操作分两步。

1、使用项集的支持度寻找频繁项集

一个项集如果使频繁项集，那么他的子集一定是频繁项集。如果一个项集是非频繁项集，那么他的超级也一定是非频繁项集。

根据以上结论，我们使用多次搜索，第k次搜索时，搜索长度为k得项集（项集得长度就是项集中含有项得数量），如：第一次搜索时，搜索长度为1得项集。在每次搜索时，我们根据事先给定的项集支持度的阈值，删除阈值以下的项集，找到第k次搜索的长度为k的频繁项集。

因为非频繁项集的超级一定是非频繁项集，这样，我们在第k+1次搜索时，就没必要搜索包含非频繁项集的项集。这样我们在寻找频繁项集时，大大减少了搜索空间。

2、使用关联规则的置信度寻找关联规则

第一步中我们得到了事物库中所有的频繁项集，我们使用频繁项集参数关联规则，我们对这些规则进行遍历，找到符合置信度要求的关联规则，在遍历时，因为不满足置信度的关联规则的子关联规则一定不满足置信度要求。所以我们在遍历每一步时，记录剔除不满足置信度的关联规则及其子关联规则。

有了上面两步，我们就可以找到强的关联规则。上面两步搜索的流程如下两图：

[图片上传中...(3.png-884e62-1625712621167-0)]

3.png

FP-growth算法

对于Apriori算法已经大大减少了计算量，FP-growth算法在Apriori算法的基础上，进一步减少了计算量，使事物库仅需要遍历两边。

FP-growth算法之所以大大减少计算量，是因为FP-growth算法使用的数据是一种数据压缩的存储形式FP树。在这种数据结构下面，我们可以很便利的找到频繁模式。从而找到强的关联规则。与此同事，FP-growth算法使用了分治策略，将一个大问题分解成几个不相交的小问题，从而使我们在处理问题时，可以只关注和我们感兴趣的子问题。

FP 树

FP树是频繁模式树，其是一种压缩存储数据的数据结构。其构造方式为增量构造。这种树的优点是，可以很容易的通过其找到频繁模式。具体构造算法如下：

1、遍历一次数据库，导出频繁项（1项集）的集合和支持度计数（频率）。并且以降序排列得到项头表L。

2、创建根节点NULL.

3、将事务库中的事物按照项头表L中项目顺序重新排了事物，注意这一步需要提出非频繁的1项集。

4、将事物库中事物一条条添加到根节点，添加方法每一个项向前面的项链接，第一个项链接根节点。当第k条事物中前面的项在前k-1条事物中出现过时，直接链接到前面的节点上，节点上记录项的同事，记录项在事务中出现的次数。

5、当遍历完所有事物时，将树上相同项用虚线相连，就形成了FP树。

FP树如下图：

微信图片_20210520102130.png

有了FP树，我们可以使用FP树找到频繁项集，但是还需要两个概念，一是条件模式基，二是条件树。

条件模式基

条件模式基就是包含FP-Tree中与后缀模式一起出现的前缀路径的集合。

例如，上图中包含项a的子节点对应的条件模式基为：{f:1,c:1,a:1,m:1}、{f:1,b:1,m:1}、{c:1,b:1,m:1}。一般的我们不写最后一个项m:1,所以项m的条件模式基是：{f:1,c:1,a:1}、{f:1,b:1}、{c:1,b:1}

条件树

有了条件模式基，我们就可以用这些基重新构造一颗FP树，接着上面的例子，项m的条件树为（设支持度的阈值为大于等于2）：

微信图片_20210520152230.png

寻找频繁项集

由上图可以得到，包含m的全部频繁模式为：{m:3}、{f:2,m:2}、{c:2,m:2}、{b:2,m:2}。同理使用递归方法，可以找到所有的频繁项集。

至此我们使用FP算找到了频繁项集，接下来，我们可以对频繁项集建立关键规则。计算每个关联规则的置信度，找到符合置信度阈值的关联规则。

至此关联分析的一些基础知识介绍完毕，我们本文中主要使用关联分析的一些概念还有FP算法，用FP算法可以高效的找到一些频繁模式和强关联规则。

接下来我们介绍另一个基础知识，共轭分布。

2.2Beta分布和多项式分布（共轭分布）

Beta分布

与贝叶斯公式紧密联系在一起的一个分布是Beta分布，这一分布有几个优良的性质，这些性质决定了Beta分布经常作为贝叶斯公式里面的先验分布。我们从Beta分布的公式推导出发，来介绍Beta分布。

Beta分布的推导

像理解二项分布一样，我们先构建一个随机事件，再求得对应随机变量的分布即为二项分布。二项分布对应的随机事件是，抛n次硬币，有k次正面向上的事件。而Beta分布对应的随机事件是，一些列随机变量的顺序统计量中，第k大的数是多少。

我们构造以下问题，来推导Beta分布。

设 $X_1,X_2,X_3...X_n$ 是区间[0,1]里面的n个随机数，服从均匀分布， $X_{（1）},X_{（2）},X_{（3）}...X_{（n）}$ 是其顺序统计量。求 $X_{（k）}$ 的分布是什么？

实际上求 $X_{（k）}$ 分布就是Beta分布，以下是推导过程。

设 $X_{1}$ 所在的区间为 $[x,x+\Delta(x)]$

当区间 $(x,x+\Delta(x))$ 里面只有 $X_{1}$ 时，

$X_{1}$ 为第k大的概率为：
$P（x+\Delta(x) >X_{1}>x ）= x^{k-1}(n-x+ \Delta(x) )^{n-k} \Delta(x) = x^{k-1}(1-x )^{n-k}\Delta(x)+o(\Delta(x) )$

其中： $o(\Delta(x))$ 为 $\Delta(x)$ 的高阶无穷小量。

当区间 $(x,x+\Delta(x))$ 里面除了 $X_{1}$ ，还有另外一个值时，则 $X_{1}$ 为第k大的概率为：
$P（x+\Delta(x) >X_{1}>x ）=\frac{1}{2} x^{k-2}(1-x+ \Delta(x) )^{n-k} \Delta(x) ^2= o(\Delta(x) )$

所以只要落在区间 $(x,x+\Delta(x))$ 里面的点大于等于2，求得对应的概率就是 $\Delta(x)$ 的高阶无穷小量。

所以，最终 $X_{1}$ 为第k大的概率为:
$P（x+\Delta(x) >X_{1}>x ） = x^{k-1}(1-x)^{n-k}\Delta(x)$

同理 $X_{i}$ 为第k大的概率也是：
$P（x+\Delta(x) >X_{i}>x ） = x^{k-1}(1-x)^{n-k} \Delta(x)$

所以有n种情况的概率总和是：
$P（x+\Delta(x) >X_{i}>x ） = n x^{k-1}(1-x)^{n-k} \Delta(x)$

又，为了保证 $X_{i}$ 为第k大，所以需要在区间[0,x]上必须有k-1个数，在剩下的n-1个数中选取k-1个数，有 $C_{n-1}^{k-1}$ 种可能。

所以最终的 $X_{(k)}$ 的概率为：
$P（x+\Delta(x) >X_{(k)}>x ）= n C_{n-1}^{k-1} x^{k-1}(1-x)^{n-k} \Delta(x)$

扩展到连续的情况, $X_{(k)}$ 的密度函数为：
$f(x)=\lim_{\Delta(x)->0} \frac{n C_{n-1}^{k-1} x^{k-1}(1-x)^{n-k} \Delta(x)}{\Delta(x)}\\ =n C_{n-1}^{k-1} x^{k-1}(1-x)^{n-k}\\ =n*\frac{(n-1)!}{(k-1)!(n-k)!}x^{k-1}(1-x)^{n-k}\\ = \frac{(n)!}{(k-1)!(n-k)!}x^{k-1}(1-x)^{n-k}\\ =\frac{\Gamma(n+1) }{ \Gamma(k) \Gamma(n-k+1) }x^{k-1}(1-x)^{n-k}$

令 $\alpha=k,\beta=n-k+1$ ,有:
$f(x) =\frac{\Gamma(\alpha+\beta) }{ \Gamma(\alpha) \Gamma(\beta) }x^{\alpha-1}(1-x)^{\beta-1}$

上面公式就是Beta分布的密度函数。至此，我们就推导除了Beta分布的密度函数。

Beta分布的均值和方差

Beta分布的均值为：
$E(x)=\frac{\alpha}{\alpha+\beta}$
Beta分布的方差为：
$Var(x)=\frac{\alpha \beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$
Beta分布的众数为：
$\frac{\alpha -1}{(\alpha+\beta-2)}$

Beta分布的性质

Beta分布有两个比较常用的性质，也正是因为这两种性质，导致在贝叶斯推断里面，Beta分布是最常见的分布。

1、Beta分布根据 $\alpha 、\beta$ 的取值不同，可以模拟的各种各样的分布。

2、Beta分布是二项分布的共轭先验分布。

对于这两个性质，我们下面展开讨论。

我们对于第一个性质，用Python画出不同 $\alpha 、\beta$ 取值的概率密度曲线图形。

import numpy as np
from scipy.stats import beta
from matplotlib import pyplot as plt

alpha_values = [1/3,2/3,1,1,2,2,4,10,20]
beta_values = [1,2/3,3,1,1,6,4,30,20]
colors =  ['blue', 'orange', 'green', 'red', 'purple', 
           'brown', 'pink', 'gray', 'olive']
x = np.linspace(0, 1, 1002)[1:-1]

fig, ax = plt.subplots(figsize=(14,9))

for a, b, c in zip(alpha_values, beta_values, colors):
    dist = beta(a, b)
    plt.plot(x, dist.pdf(x), c=c,label=r'$\alpha=%.1f,\ \beta=%.1f$' % (a, b))

plt.xlim(0, 1)
plt.ylim(0, 6)

plt.xlabel('$x$')
plt.ylabel(r'$p(x|\alpha,\beta)$')
plt.title('Beta Distribution')

ax.annotate('Beta(1/3,1)', xy=(0.014, 5), xytext=(0.04, 5.2),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(10,30)', xy=(0.276, 5), xytext=(0.3, 5.4),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(20,20)', xy=(0.5, 5), xytext=(0.52, 5.4),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(1,3)', xy=(0.06, 2.6), xytext=(0.07, 3.1),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(2,6)', xy=(0.256, 2.41), xytext=(0.2, 3.1),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(4,4)', xy=(0.53, 2.15), xytext=(0.45, 2.6),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(1,1)', xy=(0.8, 1), xytext=(0.7, 2),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(2,1)', xy=(0.9, 1.8), xytext=(0.75, 2.6),
            arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(2/3,2/3)', xy=(0.99, 2.4), xytext=(0.86, 2.8),
            arrowprops=dict(facecolor='black', arrowstyle='-'))

plt.legend(loc=0)
plt.show()

output_79_0.png

由上面beta分布的图像，我们可以看出，当 $\alpha、\beta$ 取不同值的时候，可以得到完全不一样的分布，当 $\alpha=1、\beta=1$ 时，刚好得到均匀分布。

下面我们证明Beta分布是二项分布的共轭先验，在证明之前引入几个概念。

似然（likelihood）与概率（probability）

在使用贝叶斯推断之前，必选理解这两个概念。

似然（Likelihood）：已知数据的概率值，我们反向估计所属的概率分布（参数）。

概率（Probability）：已知数据的概率分布我们去求数据的集体概率值。

贝叶斯推断

根据贝叶斯公式，我们有以下表达式：
$P(\theta|x )= \frac{ P( x| \theta) P( \theta)}{P(x)}$

其中： $P(\theta|x)$ 被称为后验概率（posterior）, $P( x|\theta)$ 被称为似然（likelihood）， $P(\theta)$ 被称为先验概率（prior）， $P(x)$ 被称为边缘似然（ marginal likelihood）可以理解为是用来归一化的。

共轭分布与共轭先验

在贝叶斯统计中，如果后验分布与先验分布属于同类（分布形式相同），则先验分布与后验分布被称为共轭分布，而先验分布被称为似然函数的共轭先验。

Beta分布是二项分布的共轭先验

设 $\theta$ 服从 $Beta(\alpha,\beta )$ 分布，x服从参数为 $\theta$ 的二项分布，即 $Bin( n,\theta)$ .

有贝叶斯公式可知， $\theta$ 的后验分布为：
$P(\theta|x ) \propto P(x| \theta)P( \theta)$

由于 $P(x| \theta)P( \theta) = C_{n}^{x} \theta^x(1-\theta)^{n-x} \frac{\Gamma(\alpha+\beta) }{ \Gamma(\alpha) \Gamma(\beta) } \theta^{\alpha-1}(1- \theta)^{\beta-1} \\ = C_{n}^{x} \frac{\Gamma(\alpha+\beta) }{ \Gamma(\alpha) \Gamma(\beta) } \theta^{x+\alpha-1}(1- \theta)^{n-x+\beta-1}\\ \propto \theta^{x+\alpha-1}(1- \theta)^{n-x+\beta-1}$

所以有：
$P(\theta|x ) \propto \theta^{x+\alpha-1}(1- \theta)^{n-x+\beta-1}$

归一化得到：
$P(\theta|x )= \frac{\Gamma(n+\alpha+\beta) }{ \Gamma(x+\alpha) \Gamma(n-x+\beta) } \theta^{x+\alpha-1}(1- \theta)^{n-x+\beta-1}$

所以 $P(\theta|x )$ 服从 $Beta(\alpha+x,\beta+n-x )$ 分布，所以Beta分布是二项分布的共轭先验。

Beta分布的两个性质，上面就介绍完了，有了这两个性质，我们可以对其进行应用，以棒球击球问题为例子，具体看一下，Beta分布怎么应用在实际中。

棒球击球率问题

一个关于击中棒球的例子。在棒球运动中，有个叫平均击球率的概念。就是用一个运动员击中棒球的次数除以他总的击球数量。一般情况下，棒球运动员的击球概率在0.266左右。高于这个值就是不错的运动员了。假设我们要预测一个运动员在某个赛季的击球率，我们可以使用已有的数据计算。但是在赛季刚开始的时候，他击球次数少，因此无法准确预测。比如他只打了一次球，那击球率就是1或者0，这个显然是不对的，我们也不会这么预测。因为我们都有一个先验期望。即根据历史情况，我们认为一个运动员大概的击球率应当是在0.215到0.360之间。因此，当一个运动员在赛季开始就被三振出局，那么我们可以预期这个运动员的击球率可能会略低于平均值，但他不可能是0。那么，在这个运动员的例子中，关于在赛季开始的击球情况，可以使用二项式分布表示，也就是一系列击球成功和失败的实验（假设之间相互独立）。同时，我们也会给这个数据一个先验期望（即统计中的先验知识），这个先验的分布一般就是Beta分布。这里的Beta分布就是用来修正我们观测到的运动员的击球率的（简单来说就是即便开始这个运动员被三振出局了，我们也只会预测他的击球率可能低于平均水平，但不会是0）。

对于以上问题，我们可以假设运动员赛季开始时击球成功与否服从二项分布Bin(n,p),运动员的击球率的先验概率服从 $Beta(\alpha,\beta)$

那么我们根据上面的贝叶斯推断有,运动员的击球率的后验概率为：
$P=Beta( \alpha+np ,\beta +n(1-p) )$

2.3Dirichlet分布和多项式分布（共轭分布）

二项分布在多维上的扩展为多项式分布，Beta分布在多维上的扩展为Dirichlet分布，所以有上面一样的结论，Dirichlet分布是多项式分布共轭先验。

Dirichlet分布

与推导Beta分布的过程一样，只不过我们在顺序统计量中，关心的不再是第k大的变量 $X_{(k)}$ 的分布，而是第k大的变量 $X_{(k)}$ 和第m大的变量 $X_{(m)}$ 的联合分布（或者更多顺序统计量的联合分布）。

与上面推导过程一样，我们可以得到这两个顺序统计量（可以更多顺序统计量）的联合分布的密度函数如下：
$f(x_1,x_2,x_3)=\frac{\Gamma(\alpha_1+\alpha_2+\alpha_3) }{\Gamma(\alpha_1)\Gamma(\alpha_2)\Gamma(\alpha_3)} x_1^{\alpha_1-1} x_2^{\alpha_2-1} x_3^{\alpha_3-1}$

其中： $x_1+x_2+x_3=1$

Dirichlet分布其实就是Beta分布在多维上的扩展，所以Dirichlet分布随着参数 $\alpha_1，\alpha_2，\alpha_3$ 的不同而形状各异。

Dirichlet分布的均值

Dirichlet分布的均值为：
$E(\vec x)=(\frac{\alpha_1}{\sum \alpha_i }、\frac{\alpha_2}{\sum \alpha_i }、...\frac{\alpha_n}{\sum \alpha_i } )$

Dirichlet分布是多项式分布的共轭先验(以二维为例)

设 $\theta_1、\theta_2$ 服从Dirichlet分布 $Dir(\alpha_1、\alpha_2、\alpha_3 )$ .x服从参数为 $\theta_1、\theta_2$ 的多项式分布 $Mult(N:\theta_1、\theta_2)$ .

有贝叶斯公式可知， $\theta$ 的后验分布为：
$P(\theta_1 ,\theta_2 |x ) \propto P(x| \theta_1 ,\theta_2)P( \theta_1 ,\theta_2)$

由于
$P(x| \theta_1 ,\theta_2)P( \theta_1 ,\theta_2) = \frac{\Gamma(\alpha_1+\alpha_2+\alpha_3) }{\Gamma(\alpha_1)\Gamma(\alpha_2)\Gamma(\alpha_3)} \theta_1^{\alpha_1-1} \theta_2^{\alpha_2-1} (1-x_2-x_1)^{\alpha_3-1} \frac{N!}{(x_1)!(x_2)!} \theta_1^{x_1}\theta_2^{x_2}$

和上面化简一样，我们可以得到后验概率为：
$P( \theta_1 ,\theta_2|x )= \frac{\Gamma(\alpha_1+\alpha_2+\alpha_3+N) }{\Gamma(\alpha_1+x_1)\Gamma(\alpha_2+x_2)\Gamma(\alpha_3+(N-x_1-x_2))}\theta_1^{\alpha_1+x_1-1} \theta_2^{\alpha_2+x_2-1} （N-\theta_1-\theta_2）^{\alpha_3+1-x_2-x_1-1}$

所以后验分布也是Dirichlet分布 $Dir(\alpha_1+x_1、\alpha_2+x_2、\alpha_3+N-x_1-x_2 )$

所以Dirichlet分布是多项式分布的共轭先验。

有了上面两个数学模型的基础知识，我们下面介绍贝叶斯规则分类器。

2.4马尔科夫链蒙特卡洛方法(MCMC)

马尔科夫链蒙特卡洛方法包含两部分，即蒙特卡洛方法和马尔科夫链。下面一一展开。

蒙特卡洛方法（随机模拟）

蒙特卡洛方法就是随机模拟，给定一个分布p(x),我们直接用这个p(x)生成一些随机样本来解决问题的方法，叫做蒙特卡洛。

怎么根据给定的分布p(x)生成随机样本。对于均匀分布，我们可以根据线性同余发生器生成伪随机数，这些伪随机数可以认为是服从均匀分布。

如何生成服从其他分布的随机数呢？我们有以下定理：

（Box-Muller定理）如果随机变量 $U_1,U_2$ 独立，且 $U_1,U_2$ 服从Unifrom[0,1].定义：
$Z_0=\sqrt{-2linU_1}cos(2\pi U_2)$
$Z_1=\sqrt{-2linU_1}sin(2\pi U_2)$
则， $Z_0,Z_1$ 独立且服从标准正态分布。

常见的连续分布，如指数分布、Gamma分布、t分布、F分布、Beta分布、Dirichlet分布都可以通过上面类似的变换得到。

蒙特卡洛方法简单使用，但是在一些情况下，蒙特卡洛方法的使用会受到限制，如：
1、分布p(x)不能显示表达时。
2、P(x)是高维分布时。

对于上面两个问题，当p(x)分布不能显示表达或者计算机不容易产生随机数时，我们一般采用拒绝接受采样的方法，对于高维分布，经常使用的方法是马尔科夫链蒙特卡洛方法(MCMC)和Gibbs Sampling

拒绝接受采样

当一个分布p(x)不能显示表达，或者计算机不容易产生随机数时，我们可以采用拒绝接受采样的方法。

拒绝接受采样的思想很简单，对于一个计算机不容易产生随机数的分布p(x)，其是先构建一个计算机容易产生随机样本的建议分布q(x)（通常选用均匀分布或者正态分布）,让密度函数q(x)的k倍kq（x）,刚刚好覆盖密度函数p(x)（如下图所示）。

微信图片_20210611112640.png

当构建好q(x)和解出k后，我们可以用q(x)产生一些随机样本 $z_0,z_1....z_n$ ，但是这些样本不一定满足p(x)，所以此时，我们要设置一个条件，满足这个条件的样本刚好满足p(x),我们就接受这些样本，不满足这些条件的样本不满足p(x),我们就拒绝这些样本。这样我们就得到了满足p(x)的随机样本。

所以问题的关键是设置的这个条件是什么。可以证明，设置的这个条件规则如下：

针对每一个样本 $z_i$ ，在Unifrom[0,1]中生存随机数U，如果：
$U<= \frac{p(z_i)}{kq(z_i)}$
则样本 $z_i$ 落在红线下的白色区域，我们接受这个样本，样本是p(x)产生的随机数。
如果：
$U> \frac{p(z_i)}{kq(z_i)}$
则样本 $z_i$ 落在红线上方的灰色区域，我们拒接这个样本，样本不是p(x)产生的随机数。

这样我们就生成了满足p(x)的随机样本。

对于为什么满足上面条件规则就是p(x)的样本，不满足就不是，本文不打算给出证明。

马尔科夫链蒙特卡洛方法(MCMC)

上面介绍了蒙特卡洛方法，为了理解MCMC，需要先理解马尔科夫链及平稳分布。

马尔科夫链

对于一个随机过程 $X_t$ ,如果存在一下公式：
$P(X_{t+1}=x|X_t,X_{t-1}.... )=P(X_{t+1}=x|X_t)$

即状态转移的概率只依赖于前一个状态。那么这个随机过程就是一个马尔科夫链。

有了马氏链，那马氏链的平稳分布如何刻画呢。有以下定理可以给出任意非周期且两两状态联通的马氏链的平稳分布。

马尔科夫链的平稳分布

（马氏链收敛定理）如果一个非周期的马氏链具有状态转移矩阵P，且他的任何两个状态是联通的，那么对于P的第i行第j列元素有：
$\lim_{n->+\infty}P_{i,j}^n$

存在且与i无关。我们记：
$\pi(j)=\lim_{n->+\infty}P_{i,j}^n$

则得到：
1、
$\lim_{n->\infty}P^n= \left[ \begin{matrix} \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \vdots & \vdots & \ddots & \vdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \end{matrix} \right]$

2、
$\pi(j)= \sum_{i=1}^{+\infty} \pi(i)P_{i,j}$

3、 $\pi$ 是方程 $\pi P=\pi$ 的唯一解。

其中： $\pi=( \pi(1),\pi(2),\pi(3)... )$ ， $\sum_{i=1}^{+\infty} \pi(i)=1$

上面的 $\pi$ 称为马氏链的平稳分布。

上面就是马氏链的收敛定理，这个定理给出了，马氏链可以在有限步之后收敛到其平稳分布，之一性质导致了MCMC方法的诞生，是MCMC方法的基础。

原始MCMC

有了上面马氏链和马氏链收敛的定理，Metropolis在1953年在处理玻尔兹曼分布采样问题时，想到一个绝妙的想法，构造了一个超级牛逼的算法，就是MCMC方法。

对于高维或者没有分布显示表达式的分布如何采样，Metropolis提出，可以构造一个马氏链，让马氏链的平稳分布刚好是p(x)，那么此马氏链在第n步之后产生的状态就是满足分布p(x)的随机样本。但是这个过程随之产生另一个问题，就是怎么构造这个马氏链，让这个马氏链的平稳分布刚好是p(x)。接下来我们一步步找到这个马氏链。

（细致平稳条件定理）如果非周期的马氏链的状态转移矩阵P和分布 $\pi(x)$ 满足
$\pi(i)P_{i,j}=\pi(j)P_{j,i}$
对任意i,j成立，则 $\pi(x)$ 是马氏链的平稳分布，上式被称为细致平稳条件。

此定理的证明不是很难，在此不展开，直接使用。

我们一开始随便找一个马氏链,对应的状态转移矩阵记为Q，很显然我们随便找的这个Q不一定满足马氏链收敛到给定p(x)的细致平稳条件，即：
$p(i)Q_{i,j} \neq p(j)Q_{j,i}$

对于以上不等式，我们引入 $\alpha(i,j)$ ,使得：
$\alpha(i,j)= p(j)Q_{j,i}$

所以有：
$p(i)Q_{i,j} \alpha(i,j)= p(j)Q_{j,i} \alpha(j,i)$

令：
$Q_{i,j}^/=Q_{i,j} \alpha(i,j)$

所以有：
$p(i)Q_{i,j}^/= p(j)Q_{j,i}^/$

所以构造出了一个新的马尔科夫链，这个马氏链状态转移矩阵是 $Q^/$ ,且满足马氏链收敛到给定p(x)的细致平稳条件。

这样我们就可以使用MCMC方法生成满足p(x)的随机样本。

MCMC的算法如下：

1、输入平稳分布p(x),输入任意的一个状态转移矩阵Q, 输入状态转移次数阈值 $n_1$ ,输入需要的样本个数 $n_2$ .

2、随机生成一个初始状态 $x_0$ .

3、for t=0 to $n_1+n_2-1:$
　　应用状态转移矩阵Ｑ生成下一步的状态 $x_*$ 。
　　在均匀分布Uniform[0,1]中生成随机数U。
　　如果 $U<\alpha(*,t)=p(x_*)Q(x_*,X_t)$ ,则接受转移 $x_t->x_*$ ,即 $x_{t+1}=x_*$ ,否则不接受转移， $x_{t}=x_*$

4、样本 $(x_{n1},x_{n1+1},x_{n1+2},.... )$ 即为满足分布p(x)的随机样本。

在MCMC原始的算法中，会存在一个问题，那就是上面第三步中 $p(x_*)Q(x_*,X_t)$ 过小，导致随机产生的U比它大，从而经常拒绝转移，导致状态原地踏步。

为解决这一问题，Hastings改进了上面的算法，也就是Metropolis-Hastings采样。

Metropolis-Hastings采样

上面MCMC方法问题所在就是 $p(x_*)Q(x_*,X_t)$ 过小，而 $p(x_*)Q(x_*,X_t)$ 刚好就是 $\alpha(*,t)$ ,而由上面的 $alpha$ 的由来可以看出，只要 $\alpha$ 满足以下平稳细致条件即可。
$p(i)Q_{i,j} \alpha(i,j)= p(j)Q_{j,i} \alpha(j,i)$

我们同时在上面等式两边乘以一个数，使得等式两边的 $\alpha$ 都等倍数的放大，就可以解决以上问题，但是不能无限放大，因为
$\alpha(i,j)= p(j)Q_{j,i}$

所以 $\alpha$ 最大只能到1，所以，我们在与U比较使用接受拒绝规则时， $\alpha$ 取以下值：
$\alpha(*,t)=min\{ \frac{ p(t)Q_{t,*} }{p(*)Q_{*,t}} ,1\}$

所以我们只需要将MCMC算法中的 $\alpha(*,t)$ 修改为上式即可。

Gibbs 采样

Gibbs采样是MCMC在高维中的形式。虽然在高维里面，但是其形式要比Metropolis-Hastings和原始的MCMC简单，因为在高维里面，满足细致平稳条件的状态转移矩阵Q可以直接找到。

以二维为例，我们推导Gibbs采样算法。

设p(x,y)是我们待采样的分布，我们直接用p(x,y)构造状态转移矩阵Q.

我们构造以下马氏链： $\{(x_1,y_1),(x_1,y_2),(x_2,y_2)... \}$

如果按以上方式构造马氏链，那么设 $(x_1,y_1)$ 为点A ， $(x_1,y_2)$ 为点B

则：
$P(A)P(y_2|x_1)=P(x_1,y_1)P(y_2|x_1)=P(x_1)P(y_1|x_1) P(y_2|x_1)$

$P(B)P(y_1|x_1)=P(x_1,y_2)P(y_1|x_1)=P(x_1)P(y_2|x_1) P(y_1|x_1)$

所以：
$P(A)P(y_2|x_1)=P(B)P(y_1|x_1)$

所以，我们如下构造状态转移矩阵Q:
$Q( (x_1,y_1)->(x_1,y_2) )= P( y_2|x_1 )$
$Q( (x_1,y_1)->(x_2,y_1) )= P( x_2|y_1 )$
$Q( (x_1,y_1)->(x_i,y_j) )= 0 对任意i,j都不为1$

所以：
$P(A)Q( (x_1,y_1)->(x_1,y_2) )=P(B)Q( (x_1,y_1)->(x_2,y_1) )$

上式就是状态转移矩阵Q的细致平稳条件。所以用这个规则(下一个状态和上一个状态有一个坐标是一样的)生成的马氏链的平稳分布是p(x,y).

具体算法如下：

1、随机初始化 $X_0=x_0,Y_0=y_0$ .
2、对t=1,2,3.... 循环采样
　　 $y_{t+1}～ p(y|x_t)$
　　 $x_{t+1}～ p(x|y_{t+1})$

最后编辑于：2021.07.08 10:57:59

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 203,271评论 5赞 476
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,275评论 2赞 380
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,151评论 0赞 336
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,550评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,553评论 5赞 365
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,559评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,924评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,580评论 0赞 257
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,826评论 1赞 297
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,578评论 2赞 320
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,661评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,363评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,940评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,926评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,156评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,872评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,391评论 2赞 342

Machine Learning-贝叶斯规则分类器(上)

贝叶斯规则分类器

目录

一、简介

二、模型相关基础知识

2.1关联分析

2.2Beta分布和二项分布（共轭分布）

2.3Dirichlet分布和多项式分布（共轭分布）

2.4马尔科夫链蒙特卡洛方法(MCMC)

三、贝叶斯规则分类器

四、模拟实验

五、模型的评估与对比

六、结束语

一、简介

二、模型相关基础知识

2.1关联分析

事物：

事物库：

项：

项集：

项集的支持度计数：

项集的支持度：

频繁项集：

关联规则：

关联规则的支持度：

关联规则的置信度：

Apriori算法

1、使用项集的支持度寻找频繁项集

2、使用关联规则的置信度寻找关联规则

FP-growth算法

FP 树

条件模式基

条件树

寻找频繁项集

2.2Beta分布和多项式分布（共轭分布）

Beta分布

Beta分布的推导

Beta分布的均值和方差

Beta分布的性质

似然（likelihood）与概率（probability）

贝叶斯推断

共轭分布与共轭先验

Beta分布是二项分布的共轭先验

棒球击球率问题

2.3Dirichlet分布和多项式分布（共轭分布）

Dirichlet分布

Dirichlet分布的均值

Dirichlet分布是多项式分布的共轭先验(以二维为例)

2.4马尔科夫链蒙特卡洛方法(MCMC)

蒙特卡洛方法（随机模拟）

拒绝接受采样

马尔科夫链蒙特卡洛方法(MCMC)

马尔科夫链

马尔科夫链的平稳分布

原始MCMC

Metropolis-Hastings采样

Gibbs 采样

推荐阅读更多精彩内容