使用python手写Metropolis-Hastings算法的贝叶斯线性回归

在学习贝叶斯计算的解马尔可夫链蒙特卡洛(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中完成:

1. import numpy as np
2. import matplotlib.pyplot as plt
3. import scipy.stats as sc
5. # sample x
6. np.random.seed(2022)
7. x = np.random.rand(100)*30
9. # set parameters
10. a = 3
11. b = 20
12. sigma = 5
14. # obtain response and add noise
15. y = a*x+b
16. noise = np.random.randn(100)*sigma
18. # create a matrix containing the predictor in the first column
19. # and the response in the second
20. data = np.vstack((x,y)).T + noise.reshape(-1,1)
22. # plot data
23. plt.scatter(data[:,0], data[:,1])
24. plt.xlabel("x")
25. plt.ylabel("y")

合成数据如下图所示：

数据的准备已经完成了，下一节将涉及定义 Metropolis Hastings 算法的函数和一组迭代次数的循环。

完整文章：

https://avoid.overfit.cn/post/fd35419dff344202af60885ce263e3cc