利用马尔可夫链蒙特卡洛（MCMC）进行贝叶斯线性回归和非线性回归的python代码（不调包）

1.利用MCMC进行线性回归

本文的特点是不利用任何市面上的贝叶斯推断的包，将全过程自己实现，利用的是M-H采样算法，从而让读者对整个过程有深刻理解。

本文呢不介绍任何数学原理。

关于线性回归数学原理的解释请看：

1. 一般的线性回归，最小二乘和最大似然估计、最大后验估计视角:
   https://www.bilibili.com/video/BV1hW41167iL?spm_id_from=333.999.0.0
2. 贝叶斯线性回归:
   https://www.bilibili.com/video/BV1St411m7XJ?spm_id_from=333.999.0.0

更加详细和全面的推导请看：《PRML》第三章。

关于MCMC的原理，请看我上一篇博文：
https://blog.csdn.net/RSstudent/article/details/122636064?spm=1001.2014.3001.5502
或查阅《PRML》等书籍的相关章节即可。

贝叶斯线性回归得到完整的后验分布，并可以给出后验分布的期望，从而避免了在最大后验估计情形下可能会出现的一些问题。

代码

import numpy as np
import scipy
import seaborn
import matplotlib.pyplot as plt"""
从概率视角来看线性回归，包括频率派视角和贝叶斯视角
本文是贝叶斯视角的，因此用的是MCMC采样，期望得到完整的后验分布而不是点估计
数学细节：bilibili上的机器学习白板推导系列
或者更详细的细节和数学推导参看：《PRML》
视频链接：
https://www.bilibili.com/video/BV1hW41167iL?spm_id_from=333.999.0.0
https://www.bilibili.com/video/BV1St411m7XJ?spm_id_from=333.999.0.0"""
# 高斯分布函数
# 默认参数设置为了二维标准高斯分布，高维也能用，只要维度对就行了
# 按道理应该检查一下用户输入的均值和方差矩阵是不是维度相符合
def guassian(x,mean=np.array([[0],[0]]),covariance = np.array([[1,0],[0,1]])):dimension = x.shape[0]# 高维if dimension > 1:guassian_kernel =\np.exp((-1/2)*np.dot(np.dot((x-mean).T,np.linalg.inv(covariance)),(x-mean)))probability = guassian_kernel/(np.power(2*np.pi,dimension/2)*np.sqrt(np.linalg.det(covariance)))# 一维，这里的协方差其实退化为方差else:guassian_kernel =\np.exp((-1/2)*(1/covariance)*(x-mean)**2)probability = guassian_kernel/(np.sqrt(2*np.pi*covariance))return probability#先验概率密度是高斯
def prior_builder(mean = np.array([[0],[0]]), covariance = np.array([[1,0], [0,1]])):def prior(parameters):if parameters.shape != (2,1):raise Exception("Wrong dimension for parameters.")probability = guassian(parameters, mean, covariance)return probabilityreturn prior# 定义似然函数
# 二维线性回归
def likelihood_builder(x,y):"""返回线性回归模型的似然函数"""def likelihood(theta):def model(x,theta):"""线性回归模型的似然函数y = theta.T*x theta是二维的，也就是有两个参数x也是二维的"""return np.dot(theta.T, x)likelihood_value = 1 # 初始化似然函数的值n_samples = x.shape[0] # 获得样本的个数，也就是数据矩阵的行数# 连乘，log下要取连加，这里不是logfor i in range(n_samples):current_x = x[i,:].reshape(2,1)current_y = y[i]guassian_mean = model(cuurent_x ,theta)probability = guassian(current_y, guassian_mean, 1)likelihood_value = likelihood_value*probabilityreturn likelihood_valuereturn likelihood# 有了似然和先验，可以定义未归一化的后验分布（密度函数）
# 采样只需要未归一化的后验，因为显然概率大的地方样本就会多，
# 只和概率相对大小有关，和绝对大小无关
def posterior_builder(likelihood, prior):# 后验分布~似然x先验,返回后验概率密度函数def posterior(theta):post_prob = likelihood(theta)*prior(theta)return post_probreturn posterior# M-H采样
def metropolis_hastings(prob_func, n_burn_in, n_samplings):"""n_burn_in: 预烧期采样次数n_smaplings: 总采样次数prob_func: 需采样的概率密度"""# 因为是二维的，所以就直接初始化一个二维向量# 这样写不能通用了，但是比较简单theta_old = np.array([[0],[0]]) # 任意初始化一个样本点# 生成一个数组，存放样本点，一共采样10000次，有10000个样本点# 但是最后这10000个样本的前1000个不需要，是预烧期的samples = np.zeros((2,n_samplings), np.float32)for i in range(n_samplings):# 进行一次随机游走# 这里原本是以旧样本为均值，方差为1的高斯分布来采样，模拟随机游走# 等价于在原来的样本上加上一个0均值，方差为1的高斯分布采样theta_new = np.random.normal(loc=theta_old, scale = 1, size = (2,1))# 计算接受率alphaalpha = np.min([prob_func(theta_new)/prob_func(theta_old),1])# 随机生成0-1随机数，小于接受率，则将新样本接收if np.random.rand() < alpha:theta_old = theta_newsamples[:,i] = theta_new.reshape(2,)# 大于接受率，则接受旧样本（注意不是舍弃）else:theta_old = theta_oldsamples[:,i] = theta_old.reshape(2,)return samplesN_BURN_IN = 3000
N_SAMPLINGS = 20000"""
模拟一批真实数据！
假设有100个样本
"""N = 100
# 真实参数
theta_real = np.array([[2.5],[6.5]],np.float32)#模拟数据
simu_x = np.random.rand(N,2) # 模拟数据矩阵x，用均匀分布生成
simu_y = np.dot(simu_x, theta_real) + np.random.randn(N,1)*0.2 #产生模拟的y# 构建在模拟数据上的似然和先验分布函数
likelihood = likelihood_builder(simu_x, simu_y)
# 给了一个方差为100的高斯，表示较弱的先验信息
para_prior = prior_builder(np.array([[3],[3]]), np.array([[100, 0], [0, 100]]))
# 得到后验分布函数
para_posterior = posterior_builder(likelihood, para_prior)# 采样
samples = metropolis_hastings(para_posterior, N_BURN_IN, N_SAMPLINGS)
seaborn.jointplot(samples[0,N_BURN_IN:], samples[1,N_BURN_IN:])# 采样结果均值，根据大数定律是后验分布的均值
print(np.mean(samples[0,N_BURN_IN:]), np.mean(samples[1,N_BURN_IN:]))
plt.show()

采样结果图：

输出的均值：

2.503646 6.5010796

2. 非线性回归

非线性回归就是非线性模型的回归，因此只需要将likelihood_builder函数中的likelihood函数中的model函数中修改为你所需要进行回归的非线性函数即可。另外后面模拟的数据也需要模拟非线性模型的数据。我们如果使用如下的非线性模型：
$$y={\theta }_1^2{x}_1+{\theta }_2{x}_2$$
则采样结果为：