共计 5995 个字符，预计需要花费 15 分钟才能阅读完成。
原文链接：http://tecdat.cn/?p=27267
咱们将钻研两种对散布进行抽样的办法：回绝抽样和应用 Metropolis Hastings 算法的马尔可夫链蒙特卡洛办法 (MCMC)。像平常一样，我将提供直观的解释、实践和一些带有代码的示例。
背景
在咱们进入主题之前，让咱们将马尔可夫链蒙特卡罗（MCMC）这个术语合成为它的根本组成部分：蒙特卡罗办法和马尔可夫链。
马尔可夫链
Markov Chain 是“在状态空间上经验从一种状态到另一种状态的转换的随机过程”。
正如你所看到的，它看起来就像一个无限状态机，只是咱们用概率正文了状态转换。例如，咱们能够查看明天是否晴天，今天晴天的概率为 0.9，下雨的概率为 0.1。同样在雨天开始。应该分明的是，从给定的状态，所有传出的转换应该总计 1.0，因为它是一个适当的散布。
示意此信息的另一种办法是通过转移矩阵 P：
将其示意为矩阵的乏味之处在于，咱们能够通过矩阵乘法来模仿马尔可夫链。例如，假如咱们从阳光明媚的状态开始，咱们能够将其示意为行向量：x(0)=[10]。这隐含地示意咱们处于晴天状态的概率为 1，因而处于下雨状态的概率为 0。当初，如果咱们执行矩阵乘法，咱们能够在一步之后找出处于每个状态的概率：
咱们能够看到今天有 0.9 的机会晴天（依据咱们的简略模型），有 0.1 的机会下雨。实际上，咱们能够持续将转换矩阵相乘，以在 k 步之后找到太阳 / 雨的机会：
咱们能够很容易地计算 x(k) 的各种 k 值，应用 numpy：
import numpy as np
Pa = nps.ardraasy(\[\[0.9, 0.1\], \[0.5, 0.5\]\])
istsdatea = np.arasdray(\[1, 0\])
siasdmulatase_markasdov(istaasdte, P, 10)
咱们能够在这里看到一个乏味的景象，当咱们在状态机中采取更多步骤时，晴天或下雨的概率仿佛会收敛。您可能认为这与咱们所处的初始状态无关，但实际上并非如此。如果咱们将初始状态初始化为随机值，咱们将失去雷同的后果：``````
siasmdulasteds_marksov(nap.arsdray(\[r, 1 - r\]), P, 10)
x^(0) = \[0.3653 0.6347\]
x^(1) = \[0.6461 0.3539\]
x^(2) = \[0.7584 0.2416\]
x^(3) = \[0.8034 0.1966\]
x^(4) = \[0.8214 0.1786\]
x^(5) = \[0.8285 0.1715\]
x^(6) = \[0.8314 0.1686\]
x^(7) = \[0.8326 0.1674\]
x^(8) = \[0.8330 0.1670\]
x^(9) = \[0.8332 0.1668\]
这种稳态散布称为 stationary distribution 通常用 π 示意。能够通过多种形式找到该稳态向量 π。最间接的办法是在 nn 靠近无穷大时取极限。
下一种办法就是求解方程。因为依据定义 q 是稳态，因而乘以 P 应该返回雷同的值：
其中 I 是单位矩阵。如果您扩大咱们的向量 / 矩阵符号，您会发现这只是一个方程组以及 π1,π2,…,πn 总和为 1 (即 π 造成概率分布）。在咱们的例子中只有两个状态：π1+π2=1。
然而，并不是每个马尔可夫链都有一个安稳的散布，甚至是惟一的 然而，如果咱们向马尔可夫链增加两个额定的束缚，咱们能够保障这些属性：
-
_不可约_：咱们必须可能最终从任何其余状态达到任何一种状态（即冀望步数是无限的）。
-
_非周期性_：零碎永远不会返回到具备固定周期的雷同状态（例如，不会每 5 步确定性地返回开始“晴天”）。
一个重要的定理说，如果马尔可夫链是遍历的，那么它有一个惟一的稳态概率向量 ππ。在 MCMC 的上下文中，咱们能够从任何状态跳转到任何其余状态（以肯定的无限概率），轻松满足不可约性。
咱们将应用的另一个有用的定义是 detailed balance and reversible Markov Chains. 如果存在满足此条件的概率分布 π，则称马尔可夫链是可逆的（也称为具体平衡条件）：
换句话说，从久远来看，你从状态 i 转移到状态 j 的次数比例，与你从状态 j 转移到状态 i 的次数比例雷同。事实上，如果马尔可夫链是可逆的，那么咱们就晓得它具备安稳散布（这就是咱们应用雷同符号 π 的起因）。
马尔可夫链蒙特卡罗办法
马尔可夫链蒙特卡罗 (MCMC) 办法只是一类应用马尔可夫链从特定概率分布（蒙特卡罗局部）中采样的算法。他们通过创立一个马尔可夫链来工作，其中限度散布（或安稳散布）只是咱们想要采样的散布。
这是一张可能有助于形容该过程的图片. 设想一下，咱们正在尝试制作一个 MCMC 来尝试应用 PDF f(x)对任意一维散布进行采样。在这种状况下，咱们的状态将是沿 x 轴的点，而咱们的转换概率将是从一种状态到另一种状态的机会。这是状况的简化图：
该图显示了咱们试图用粗黑线迫近的密度函数，以及应用从橙色状态过渡的蓝线的马尔可夫链的一部分的可视化。特地是，对于 i=-3,-2,-1,1,2,3，只是从状态 X0 到 Xi 的转换。然而，x 轴线上的每个点实际上都是这个马尔可夫链中的一个势态。请留神，这意味着咱们有一个有限的状态空间，因而咱们不能再将转换很好地示意为矩阵。MCMC 办法的真正“窍门”是咱们想要设计状态（或 x 轴上的点）之间的转换概率，以便咱们将大部分工夫花在 f(x) 很大的区域中，并且在它较小的区域中的工夫绝对较少（即与咱们的密度函数成准确比例）。
就咱们的人物而言，咱们心愿将大部分工夫花在核心四周，而较少工夫花在里面。事实上，如果咱们模仿咱们的马尔可夫链足够长的工夫，状态的限度散布应该近似于咱们试图采样的 PDF。因而，应用 MCMC 办法进行采样的根本算法为：
-
从任意点 x 开始。
-
以肯定的转移概率跳转到点 x’（这可能意味着放弃雷同的状态）。
-
转到第 2 步，直到咱们转换了 T 工夫。
-
记录以后状态 x′，进行步骤 2。
当初，咱们在每个点 x 轴上破费的比例次数应该是咱们试图模仿的 PDF 的近似值，即如果咱们绘制 x 值的直方图，咱们应该失去雷同的形态。
回绝抽样
当初，在咱们进入 MCMC 办法的具体算法之前，我想介绍另一种对概率分布进行采样的办法，咱们稍后将应用它，称为 rejection sampling. 次要思维是，如果咱们试图从散布 f(x) 中采样，咱们将应用另一个工具散布 g(x) 来帮忙从 f(x) 中采样。惟一的限度是对于某些 M>1，f(x)<Mg(x)。它的主要用途是当 f(x) 的模式难以间接采样时（但依然能够在任何点 xx 对其进行评估）。
以下是算法的细分：
-
从 g(x) 中采样 x。
-
从 U(0,Mg(x)) 中采样 y（均匀分布）。
-
如果 y<f(x)，则承受 x 作为 f(x) 的样本，否则转到步骤 1。
另一种对待它的办法是咱们采样点 x0 的概率。这与从 g 中采样 x0 的概率乘以咱们承受的次数的比例成正比，它简略地由 f(x0) 和 Mg(x0) 之间的比率给出：
等式通知咱们对任意点进行采样的概率与 f(x0) 成正比。在对许多点进行采样并找到咱们看到 x0 的次数比例之后，常数 M 被归一化，咱们失去了 PDF f(x) 的正确后果。
让咱们通过一个例子更直观地看一下它。咱们要从中采样的指标散布 f(x) 是 double gamma 散布，基本上是一个双边伽马散布。咱们将应用正态分布 g(x) 作为咱们的包络散布。上面的代码向咱们展现了如何找到缩放常数 M，并为咱们描述了回绝抽样在概念上是如何工作的。
# 指标 = 双伽马散布
dsg = stats.dgamma(a=1)
# 生成 PDF 的样本
x = np.linspace
# 绘图
ax = df.plot(style=\['--', '-'\]
从图中，一旦咱们找到 g(x)的样本（在这种状况下 x=2），咱们从范畴等于 Mg(x) 高度的均匀分布中绘制. 如果它在指标 PDF 的高度内，咱们承受它（绿色），否则回绝（回绝）。
施行回绝抽样
上面的代码为咱们的指标双伽马散布实现回绝采样。它绘制标准化直方图并将其与咱们应该失去的实践 PDF 匹配。
# 从回绝采样算法生成样本
sdampales = \[rejeasdctioan_samplaing for x in range(10000)\]
# 绘制直方图与指标 PDF
df\['Target'\].plot
总的来说，咱们的回绝采样器非常适合。与实践散布相比，抽取更多样本会改善拟合。
回绝抽样的很大一部分是它很容易实现（在 Python 中只需几行代码），但有一个次要毛病：它很慢。
Metropolis-Hastings 算法
这 Metropolis-Hastings Algorithm (MH) 是一种 MCMC 技术，它从难以间接采样的概率分布中抽取样本。与回绝抽样相比，对 MH 的限度实际上更加宽松：对于给定的概率密度函数 p(x)，咱们只要求咱们有一个 与 p_成正比_的函数 f(x)f(x) (x)p(x)！这在对后验散布进行采样时十分有用。
Metropolis-Hastings 算法的推导
为了推导出 Metropolis-Hastings 算法，咱们首先从最终目标开始：创立一个马尔可夫链，其中稳态散布等于咱们的指标散布 p(x)。就马尔可夫链而言，咱们曾经晓得状态空间将是什么：概率分布的反对，即 x 值。因而（假如马尔可夫链的结构正确）咱们最终失去的稳态散布将只是 p(x)。剩下的是确定这些 x 值之间的转换概率，以便咱们能够实现这种稳态行为。
马尔可夫链的具体平衡条件，这里用另一种形式写成：
这里 p(x)是咱们的指标散布，P(x→x′) 是从点 x 到点 x′ 的转移概率。所以咱们的指标是确定 P(x→x′)的模式。既然咱们要构建马尔可夫链，让咱们从应用等式 5 作为该构建的根底开始。请记住，具体的平衡条件保障咱们的马尔可夫链将具备安稳散布（它存在）。此外，如果咱们也包含遍历性（不以固定距离反复状态，并且每个状态最终都可能达到任何其余状态），咱们将建设一个具备惟一安稳散布 p(x)的马尔可夫链.
咱们能够将方程重新排列为：
这里咱们应用 f(x)来示意一个 与 p(x)_成正比的函数。_这是为了强调咱们并不明确须要 p(x)，只是须要与它成比例的货色，这样比率能力达到雷同的成果。当初这里的“技巧”是咱们将把 P(x→x′)合成为两个独立的步骤：一个提议散布 g(x→x′) 和承受散布 A(x→x′)（相似于回绝抽样的工作原理）。因为它们是独立的，咱们的转移概率只是两者的乘积：
此时，咱们必须弄清楚 g(x)和 A(x) 的适合抉择是什么。因为 g(x) 是“倡议散布”，它决定了咱们可能采样的下一个点。因而，重要的是它具备与咱们的指标散布 p(x)（遍历性条件）雷同的反对。这里的一个典型抉择是以以后状态为核心的正态分布。当初给定一个固定的提议散布 g(x)，咱们心愿找到一个匹配的 A(x)。
尽管不显著，但满足公式 的 A(x) 的典型抉择是：
咱们能够通过思考 f(x′)g(x′→x)小于等于 1 和大于 1 的状况。当小于等于 1 时，它的倒数大于 1，因而 LHS 的分母，A(x′→ x), 等式 8 为 1，而分子等于 RHS。或者，当 f(x′)g(x′→x)是大于 1 LHS 的分子是 1，而分母正好是 RHS 的倒数，导致 LHS 等于 RHS。
这样，咱们曾经证实，咱们创立的马尔可夫链的稳固状态将等于咱们的指标散布 (p(x))，因为具体的平衡条件通过结构失去满足。所以整体算法将是（与下面的 MCMC 算法十分匹配）：
-
通过抉择一个随机 x 来初始化初始状态。
-
依据 g(x→x′)找到新的 x′。
-
依据 A(x→x′) 以平均概率承受 x′。如果承受到 x’ 的转换，否则放弃状态 x。
-
进行第 2 步，T 次。
-
将状态 x 保留为样本，转至步骤 2 对另一个点进行采样。
预烧和相干样本
在咱们持续实现之前，咱们须要探讨对于 MCMC 办法的两个十分重要的话题。第一个主题与咱们抉择的初始状态无关。因为咱们随机抉择 xx 的值，它很可能位于 p(x) 十分小的区域（想想咱们的双伽马散布的尾部）。如果从这里开始，它可能会破费不成比例的工夫来遍历低密度的 x 值，从而谬误地给咱们一种感觉，即这些 x 值应该比它们更频繁地呈现。所以解决这个问题的办法是“预烧”采样器通过生成一堆样本并将它们扔掉。样本的数量将取决于咱们试图模仿的散布的细节。
咱们下面提到的第二个问题是两个相邻样本之间的相关性。因为依据咱们的转换函数 P(x→x′)的定义，绘制 x′ 取决于以后状态 x。因而，咱们失去了样本的一项重要属性：独立性。为了纠正这一点，咱们抽取 Tth 个样本，并且只记录最初抽取的样本。假如 T 足够大，样本应该是绝对独立的。与预烧一样，T 的值取决于指标和提议散布。
实现 Metropolis-Hastings 算法
让咱们应用下面的双伽马散布示例。让咱们将咱们的提议散布定义为以 x 为核心、标准差为 2、N(x, 2) 的正态分布（记住 x 是以后状态）：
给定 f(x) 与咱们的根底散布 p(x) 成比例，咱们的承受散布如下所示：
因为正态分布是对称的，因而正态分布的 PDF 在其各自点进行评估时会互相对消。当初让咱们看一些代码：
# 模仿与双伽马散布成比例的 f(x)
f = ambd x: g.df(x* mat.i
采样器 = mhspler()
# 样本
sames = \[nex(saper) for x in range(10000)\]
# 绘制直方图与指标 PDF
df\['Target'\].plot
来自咱们的 MH 采样器的样本很好地近似于咱们的双伽马散布。此外，查看自相干图，咱们能够看到它在整个样本中十分小，表明它们是绝对独立的。如果咱们没有为 T 抉择一个好的值或没有预烧期，咱们可能会在图中看到较大的值。
论断
我心愿你喜爱这篇对于应用回绝抽样和应用 Metropolis-Hastings 算法进行 MCMC 抽样的简短文章。
最受欢迎的见解
1. 应用 R 语言进行 METROPLIS-IN-GIBBS 采样和 MCMC 运行
2.R 语言中的 Stan 概率编程 MCMC 采样的贝叶斯模型
3.R 语言实现 MCMC 中的 Metropolis–Hastings 算法与吉布斯采样
4.R 语言 BUGS JAGS 贝叶斯剖析 马尔科夫链蒙特卡洛办法（MCMC）采样
5.R 语言中的 block Gibbs 吉布斯采样贝叶斯多元线性回归
6.R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真剖析
7.R 语言用 Rcpp 减速 Metropolis-Hastings 抽样预计贝叶斯逻辑回归模型的参数
8.R 语言应用 Metropolis- Hasting 抽样算法进行逻辑回归
9. R 语言中基于混合数据抽样 (MIDAS) 回归的 HAR-RV 模型预测 GDP 增长回归的 HAR-RV 模型预测 GDP 增长 ”)