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

原创

已于 2022-02-02 11:34:29 修改 · 7.5k 阅读

65 ·

CC 4.0 BY-SA版权

文章标签：

#python #线性回归 #概率论

于 2022-01-26 17:35:57 首次发布

本文通过MCMC方法实现线性回归，不依赖现有贝叶斯包，自行完成采样过程，深入理解贝叶斯视角下的线性回归。采用Metropolis-Hastings算法获取后验分布，提供完整代码及非线性回归扩展方法。

1.利用MCMC进行线性回归

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

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

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

一般的线性回归，最小二乘和最大似然估计、最大后验估计视角:
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

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

关于MCMC的原理，请看我上一篇博文：
https://blog.youkuaiyun.com/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(

最低0.47元/天解锁文章