点击上方“Deephub Imba”,关注公众号,好文章不错过 !
在学习贝叶斯计算的解马尔可夫链蒙特卡洛(MCMC)模拟时，最简单的方法是使用PyMC3，构建模型，调用Metropolis优化器。但是使用别人的包我们并不真正理解发生了什么，所以本文通过手写Metropolis-Hastings来深入的理解MCMC的过程，再次强调我们自己实现该方法并不是并不是为了造轮子，而是为了更好的通过代码理解该概念。
贝叶斯线性回归包含了几十个概念和定义，这使得我们的整个研究成为一种折磨，并且真正发生的事情。在本文中，我将通过常见Metropolis-Hastings 算法构建一个马尔可夫链，并提供一个实际的使用案例。我们将着重于推断简单线性回归模型的参数（但是这里说“简单”并不能代表它背后的原理简单）。我们还可以用任何其他条件分布(泊松/伽马/负二项或其他)替换正态分布，这样可以通过MCMC实现几乎相同的GLM（只更改4或5行代码）。由于推断的样本来自参数的后验分布，我们可以很自然地使用这些来构建参数的置信区间(在这种情况下，“可信区间”更正确)。
下面我们将简要描述为什么使用MCMC方法，提供一个线性回归模型的MH算法的实现，并将以一个可视化的方式显示当算法寻找生成数据的参数集时，真正发生了什么。
数据准备
设Y和X分别为模型的响应和输入。线性模型表明，给定输入的响应的条件分布是正态的。也就是:
对于合适的参数a(斜率)、b(偏差)和σ(噪声强度)。
我们的任务是推断a, b和σ。
所以我们首先要知道一些模型需要遵循的“基本规则”。在这个例子中，a, b几乎可以是任何数值，正的或负的，但σ必须是严格正的(因为从来没有听说过负标准差的正态分布，对吧?)除此之外，没有其他任何规则。也许你会说：“我们需要先了解这些是如何分布的”，但是后验分布的渐近正态性保证告诉我们，只要有足够的后验样本，这些样本无论如何都是正态分布(如果马尔可夫链达到其平稳分布)，所以分布不是我们考虑的必要因素。
现在让我们为回归生成合成数据，这里使用参数a=3, b=20和σ=5。可以通过以下代码在python中完成:
代码语言：javascript复制import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sc
# sample x
np.random.seed(2022)
x = np.random.rand(100)*30
# set parameters
a = 3
b = 20
sigma = 5
# obtain response and add noise
y = a*x b
noise = np.random.randn(100)*sigma
# create a matrix containing the predictor in the first column
# and the response in the second
data = np.vstack((x,y)).T noise.reshape(-1,1)
# plot data
plt.scatter(data[:,0], data[:,1])
plt.xlabel("x")
plt.ylabel("y")合成数据如下图所示：
数据的准备已经完成了，下一节将涉及定义 Metropolis Hastings 算法的函数和一组迭代次数的循环。
算法介绍
假设θ=[a,b，σ]是算法上面的参数向量，θ '是一组新参数的建议，MH比较参数(θ '和θ)的两个竞争假设之间的贝叶斯因子(似然和先验的乘积)，并通过条件建议分布的倒数缩放该因子。然后将该因子与均匀分布的随机变量的值进行比较。这给模型增加了随机性，使不可能的参数向量有可能被探索，也可能被丢弃(很少)。
这听起来有点复杂，让我们从头一步一步对它进行代码的实现。
Proposal Distribution
首先，我们定义了一个建议的分布（Proposal Distribution）g(θ′|θ):这是前一个时间步固定参数的分布。在我们的例子中，a和b可以是正的也可以是负的，因此一个自然的选择就是从一个以前一个迭代步骤为中心的多元正态分布中采样它们。我们可以从伽马分布中取样σ，这些分布的定义我们可以根据实际情况进行选择，但是一个更好的方法(这里我们将不涉及)是从逆伽马分布抽样σ。
因此，我们可以按照以下方式定义进行Proposal Distribution:
分布抽样σ为
参数 k 用于控制分布在其均值周围的“扩展”。
