利用马尔可夫链蒙特卡洛（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 probability
    return 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下要取连加，这里不是log
        for 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*probability

        return likelihood_value
    
    return likelihood

# 有了似然和先验，可以定义未归一化的后验分布（密度函数）
# 采样只需要未归一化的后验，因为显然概率大的地方样本就会多，
# 只和概率相对大小有关，和绝对大小无关
def posterior_builder(likelihood, prior):
    # 后验分布~似然x先验,返回后验概率密度函数
    def posterior(theta):
        post_prob = likelihood(theta)*prior(theta)
        return post_prob
    return 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))
        # 计算接受率alpha
        alpha = np.min([prob_func(theta_new)/prob_func(theta_old),1])
        # 随机生成0-1随机数，小于接受率，则将新样本接收
        if np.random.rand() < alpha:
            theta_old = theta_new
            samples[:,i] = theta_new.reshape(2,)
        # 大于接受率，则接受旧样本（注意不是舍弃）
        else:
            theta_old = theta_old
            samples[:,i] = theta_old.reshape(2,)

    return samples

N_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$$
则采样结果为：