原文链接：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 npPa = 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)\]# 绘制直方图与指标 PDFdf\['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)\]# 绘制直方图与指标 PDFdf\['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增长")