频率主义线性回归和贝叶斯线性回归

整体概述

频率主义（Frequentist）线性回归和贝叶斯（Bayesian）线性回归是统计学中用于数据分析和预测的两种主要方法，特别是在建模关于因变量和自变量之间线性关系的上下文中。尽管两种方法都用于线性回归分析，但它们在解释参数、估计方法和结果的不确定性方面有所不同。

频率主义线性回归

频率主义线性回归建立在经典统计学的原则上，尤其是参数估计的频率主义概念。在频率主义视角下，模型参数被认为是固定的未知常数，我们的目标是基于样本数据去估计这些常数。频率主义线性回归的目标是找到一个线性模型来描述因变量Y和一组自变量X之间的关系：
$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\dots +{\beta }_p{X}_p+ϵ$
其中，${\beta }_0,{\beta }_1,...,{\beta }_p$是模型参数，而 $ϵ$是误差项，通常假设为独立同分布的正态分布。
参数估计常用最小二乘法（OLS），目标是最小化误差项的平方和。频率主义的置信区间和假设检验（如t检验、F检验）用于评估模型参数的不确定性和显著性。

贝叶斯线性回归

贝叶斯线性回归遵循贝叶斯统计的原则，模型的参数被视为随机变量而非固定未知常数。这种方法将先验知识或者信念和样本数据结合起来，以估计参数。贝叶斯回归的基本形式与频率主义线性回归相同，但在估计参数时采用的是贝叶斯定理：
$P(\Theta |Data)=\frac{P(Data|\Theta )\cdot P(\Theta )}{P(Data)}$
这里的 $P(\Theta |Data)$ 是给定数据后参数$\Theta$ 的后验概率，$P(Data|\Theta )$是给定参数下数据的似然，$P(\Theta )$ 是参数的先验概率，而$P(Data)$ 是证据（或边缘概率），作为归一化常数。
在贝叶斯回归中，不是寻找一组固定的参数估计，而是利用计算或模拟方法（如马尔可夫链蒙特卡洛，MCMC）来得到参数的后验分布。这允许我们直接从后验分布中得到参数估计的不确定性（表现为分布的形状），也能通过分布的特征（如均值、中位数、百分位数区间）来总结参数的估计。

比较

• 参数解释：频率主义视角下，参数是固定的，而在贝叶斯视角下，参数是具有分布的随机变量。
• 先验知识：贝叶斯方法可以整合先验知识，频率主义方法则通常不考虑先验知识。
• 不确定性评估：贝叶斯方法通过后验分布直接提供参数的不确定性评估，而频率主义方法使用置信区间和假设检验来进行不确定性评估。
• 计算复杂性：贝叶斯方法通常计算上更为复杂，需要使用MCMC等模拟技术，而频率主义方法如OLS在计算上更简单直接。

每种方法都有其优势和局限性，应用哪种方法取决于具体的问题、数据和分析目标。

频率主义线性回归

频率主义线性回归是一种统计方法，用于估计和推断一个或多个自变量（解释变量）与因变量（响应变量）之间的线性关系。在频率主义框架中，参数被认为是固定的，且数据是随机的。这一观点反映了一个信念，即在重复实验的条件下，某事件的频率趋近于其真实概率。

线性回归模型

线性回归模型可以表示为：
$Y=X\beta +ϵ$
其中，

• $Y$ 是一个 $n\times 1$的向量，表示因变量的观测值；
• $X$ 是一个$n\times p$ 的矩阵，其中 $n$是样本大小，$p$是参数的数量（包括截距），表示自变量；
• $\beta$是一个 $p\times 1$ 的向量，表示未知的回归系数；
• $ϵ$ 是一个 $n\times 1$的向量，表示误差项，假设$ϵ\sim N(0,{\sigma }^{2}I)$，即误差项独立同分布，且呈正态分布。

最小二乘估计

在频率主义线性回归中，通常使用最小二乘法（Ordinary Least Squares, OLS）来估计参数$\beta$。最小二乘估计旨在最小化误差项的平方和：
$S(\beta )=(Y-X\beta {)}^{\prime }(Y-X\beta )$
通过对 $S(\beta )$求导并设其为零，可以得到 $\beta$ 的OLS估计：
$\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }Y$
其中$\hat{\beta }$是参数的估计值。

置信区间

置信区间是用来估计模型参数不确定性的一种方法。对于每个回归系数 ${\beta }_i$，我们可以计算一个置信区间，该区间以一定的置信水平（如95%）覆盖真实参数值。置信区间基于回归系数的抽样分布，通常假设为正态分布，可以表示为：
$\widehat{{\beta }_i}\pm t_{\alpha /2,n-p}\cdot SE(\widehat{{\beta }_i})$
其中，

• $t_{\alpha /2,n-p}$ 是自由度为 $n-p$ 的t分布的临界值；
• $SE(\widehat{{\beta }_i})$ 是估计系数 $\widehat{{\beta }_i}$的标准误差。

假设检验

假设检验用于评估模型系数是否显著不同于零（或其他指定的值）。最常用的是t检验，用于单个系数，以及F检验，用于模型中的多个系数。t检验的零假设是系数等于零$(H_0:{\beta }_i=0)$，备择假设是系数不等于零$(H_1:{\beta }_i\ne 0)$。
t统计量计算为：
$t=\frac{\widehat{{\beta }_i}-0}{SE(\widehat{{\beta }_i})}$
如果 (t) 的绝对值大于t分布的临界值，我们拒绝零假设，认为这个系数显著不为零。
F检验通常用于同时检验多个系数是否显著。它比较模型中所有解释变量与只有截距项的简化模型的拟合度。F统计量是比例，反映了两个模型中未解释变化的比例减少：
$F=\frac{(SS{E}_0-SS{E}_1)/(p-1)}{SS{E}_1/(n-p)}$
其中，

• $SS{E}_0$ 是没有解释变量的模型的误差平方和；
• $SS{E}_1$ 是完整模型的误差平方和；
• $p$ 是完整模型中参数的数量；
• $n$ 是样本大小。

如果F统计量的值大于F分布的临界值，我们拒绝零假设，认为模型中至少有一个系数显著不为零。
通过这些方法，频率主义线性回归不仅为我们提供了对数据的最佳线性拟合，而且还允许我们进行推断，评估模型系数的显著性以及估计系数的可靠性。

贝叶斯线性回归

贝叶斯线性回归通过结合先验分布和数据的似然来更新关于回归参数的信念，得到参数的后验分布。由于后验分布通常难以解析计算，尤其是在参数空间较大或模型复杂时，马尔可夫链蒙特卡洛（MCMC）方法被广泛用于生成后验分布的近似样本。这些样本可以用于参数估计、预测以及计算可信区间。
以下是贝叶斯线性回归估计和预测的详细步骤，包括使用MCMC方法和计算可信区间的说明。

模型设定

首先，设定一个贝叶斯线性回归模型：
$Y=X\beta +ϵ,ϵ\sim N(0,{\sigma }^{2}I)$

• $Y$是观测的响应变量；
• $X$ 是设计矩阵，包含自变量；
• $\beta$ 是回归系数，需要估计；
• $ϵ$是误差项，假设服从以0为均值、${\sigma }^2$为方差的正态分布。

先验分布设定

接着，给定参数$\beta$ 和 ${\sigma }^2$ 的先验分布。这些先验可以是无信息先验，例如正态先验或者均匀先验，也可以是包含先前知识的信息先验。

后验分布

根据贝叶斯定理，后验分布与似然函数和先验分布成正比：
$p(\beta ,{\sigma }^2|Y,X)\propto p(Y|X,\beta ,{\sigma }^2)p(\beta )p({\sigma }^2)$
其中 $p(\beta ,{\sigma }^2|Y,X)$ 是后验分布，$p(Y|X,\beta ,{\sigma }^2)$ 是似然函数，$p(\beta )$ 和 $p({\sigma }^2)$是先验分布。

MCMC采样

MCMC方法允许我们从复杂的后验分布中抽取样本。最常用的MCMC算法是吉布斯采样（Gibbs sampling）和梅特罗波利斯-哈斯廷斯算法（Metropolis-Hastings algorithm）。这些算法通过构建一个随机游走的马尔可夫链，最终收敛到后验分布。

• 吉布斯采样在每个步骤中条件于其他所有参数对每个参数进行采样。
• 梅特罗波利斯-哈斯廷斯算法可以采样整个参数向量，并利用接受-拒绝规则来确保样本来自正确的分布。

后验分析

当MCMC收敛后，我们可以通过分析生成的样本来估计后验分布。这包括：

• 参数估计：通过取MCMC样本的平均值或中位数作为参数的估计值。
• 预测：使用样本生成预测分布，通过在每个采样步骤上将新数据点的预测值 (X_{new}\beta) 累积起来并应用后验分布。
• 可信区间计算：计算参数的可信区间，这可以通过对MCMC样本进行排序并取出相应的百分位数来实现。例如，95%的可信区间由2.5%和97.5%的分位数确定。

可信区间（Credible Intervals）

与频率主义的置信区间不同，贝叶斯的可信区间提供了一个区间，使得有高概率（如95%）参数的真实值会落在这个区间内。计算步骤如下：

• 对每个参数的MCMC样本进行排序。
• 对于95%的可信区间，选取2.5%的分位数作为下限，97.5%的分位数作为上限。

诊断和验证

最后，对MCMC算法进行诊断，确保它已经收敛到目标后验分布。可以使用迹图、自相关图以及收敛诊断统计量，如Gelman-Rubin统计量进行检查。
综上所述，贝叶斯线性回归利用MCMC方法生成后验分布的样本，允许我们进行参数估计、预测以及可信区间的计算。这为包含不确定性和先验信念的统计建模提供了一种强大的工具。

示例：

例题

根据$x_0,x_1$值预测$y$值。我们划分为训练集、验证集和测试集。数据保存在train_set.csv和test_set.csv。使用以下评价指标：R²、MSE、RMSE、MAE。
R平方 (R²) 或决定系数：衡量模型解释变量能力的比例，值在0到1之间，1表示模型能够完美预测因变量。
均方误差 (MSE)：衡量模型预测误差的平均值的大小，计算公式是预测值与实际值差的平方的平均。
均方根误差 (RMSE)：MSE的平方根，它的量纲与因变量一致，因此更易于解释。
平均绝对误差 (MAE)：衡量预测值与实际值之间差的绝对值的平均大小，与MSE类似，但是对异常值的敏感度较低。

频率主义线性回归

首先我们使用scikit-learn的线性回归模型进行学习和推断。

import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score# 加载数据
data = pd.read_csv('train_set.csv')  
test_data = pd.read_csv('test_set.csv')  # 将x0和x1两列转换为一个Numpy数组
train_X = data[['x0', 'x1']].values
test_X = test_data[['x0', 'x1']].values
# 将y列转换为一个Numpy数组
train_y = data['y'].values
test_y = test_data[['y']].values
# 现在X是一个包含x0和x1的Numpy数组，y是一个单独的Numpy数组
print(train_X.shape,train_y.shape)
print(test_X.shape,test_X.shape)
# 创建模型
model = LinearRegression()
# 拟合模型
model.fit(train_X, train_y)
# 输出回归系数
print(f"Intercept: {model.intercept_}, Coefficients: {model.coef_}")
# 使用模型进行预测
y_pred = model.predict(test_X)# 计算统计指标
mse = mean_squared_error(test_y, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(test_y, y_pred)
r2 = r2_score(test_y, y_pred)# 打印统计指标
print(f'Mean Squared Error (MSE): {mse}')
print(f'Root Mean Squared Error (RMSE): {rmse}')
print(f'Mean Absolute Error (MAE): {mae}')
print(f'R-squared (R2): {r2}')'''
(400, 2) (400,)
(100, 2) (100, 2)
Intercept: 0.06634575026290213, Coefficients: [1.15689086 0.88110244]
Mean Squared Error (MSE): 0.0002869146993285752
Root Mean Squared Error (RMSE): 0.016938556589289867
Mean Absolute Error (MAE): 0.013865979260745431
R-squared (R2): 0.8927043803619303
'''

不使用任何专门的库，我们可以直接用Python中的基础功能来计算线性回归模型的参数。这可以通过计算最小二乘估计来实现，即通过解析解找到系数β的估计值，使得残差平方和最小。
线性回归模型的解析解可通过以下公式计算得到：
$\beta =(X^{T}X{)}^{-1}{X}^{T}y$
这里的$\beta$ 是一个向量，包含了截距项和斜率项。
以下是一个用Python实现的简单线性回归模型的例子：

import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score# 加载数据
data = pd.read_csv('train_set.csv')  
test_data = pd.read_csv('test_set.csv')  # 将x0和x1两列转换为一个Numpy数组
train_X = data[['x0', 'x1']].values
test_X = test_data[['x0', 'x1']].values
# 将y列转换为一个Numpy数组
train_y = data['y'].values
test_y = test_data[['y']].values
# 现在X是一个包含x0和x1的Numpy数组，y是一个单独的Numpy数组
print(train_X.shape,train_y.shape)
print(test_X.shape,test_X.shape)# 添加截距项
train_X_with_intercept = np.column_stack((np.ones((train_X.shape[0], 1)), train_X))# 计算β的估计值
beta_hat = np.linalg.inv(train_X_with_intercept.T @ train_X_with_intercept) @ train_X_with_intercept.T @ train_yprint(f"Intercept (beta_0): {beta_hat[0]}")
print(f"Coefficient of x0 (beta_1): {beta_hat[1]}")
print(f"Coefficient of x1 (beta_2): {beta_hat[2]}")# 使用模型进行预测
def predict(X, beta):X_with_intercept = np.hstack((np.ones((X.shape[0], 1)), X))y_pred = X_with_intercept @ betareturn y_pred
y_pred = predict(test_X,beta_hat)# 计算统计指标
mse = mean_squared_error(test_y, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(test_y, y_pred)
r2 = r2_score(test_y, y_pred)# 打印统计指标
print(f'Mean Squared Error (MSE): {mse}')
print(f'Root Mean Squared Error (RMSE): {rmse}')
print(f'Mean Absolute Error (MAE): {mae}')
print(f'R-squared (R2): {r2}')'''
(400, 2) (400,)
(100, 2) (100, 2)
Intercept (beta_0): 0.06634575026264328
Coefficient of x0 (beta_1): 1.156890863218591
Coefficient of x1 (beta_2): 0.8811024418527595
Mean Squared Error (MSE): 0.00028691469932902837
Root Mean Squared Error (RMSE): 0.016938556589303246
Mean Absolute Error (MAE): 0.013865979260770434
R-squared (R2): 0.8927043803617608
'''

在上述代码中，np.column_stack用来合并截距项和自变量。这种方法没有考虑矩阵(X^TX)可能不可逆的情况。在实际应用中，通常会加入一些正则化或者使用伪逆（np.linalg.pinv）来避免这种问题。

贝叶斯线性回归

我们使用Pyro库实现贝叶斯线性回归

import torch
import pyro
import pyro.distributions as dist
from pyro.infer import MCMC, NUTS
import pandas as pd
import numpy as np 
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score# 加载数据
data = pd.read_csv('train_set.csv')  
test_data = pd.read_csv('test_set.csv')  # 将x0和x1两列转换为一个Numpy数组
train_X0 = torch.from_numpy(data['x0'].values) 
train_X1 = torch.from_numpy(data['x1'].values)
test_X0 = torch.from_numpy(test_data['x0'].values)
test_X1 = torch.from_numpy(test_data['x1'].values)
# 将y列转换为一个Numpy数组
train_y = torch.from_numpy(data['y'].values)
test_y = torch.from_numpy(test_data['y'].values)
print(test_X0.shape,test_X1.shape)def model(x0, x1, y=None):# 定义先验分布beta0 = pyro.sample("beta0", dist.Normal(0., 1.))beta1 = pyro.sample("beta1", dist.Normal(0., 1.))beta2 = pyro.sample("beta2", dist.Normal(0., 1.))sigma = pyro.sample("sigma", dist.Uniform(0., 10.))# 定义线性模型mean = beta0 + beta1 * x0 + beta2 * x1# 用于观测数据的似然函数with pyro.plate("data", len(x0)):obs = pyro.sample("obs", dist.Normal(mean, sigma), obs=y)return mean# 运行Hamiltonian Monte Carlo (HMC)算法，NUTS是HMC的一个优化版本
nuts_kernel = NUTS(model)
mcmc = MCMC(nuts_kernel, num_samples=500, warmup_steps=300)
mcmc.run(train_X0, train_X1, train_y)# 获取后验分布的样本
posterior_samples = mcmc.get_samples()# 打印后验分布的平均值
for name, value in posterior_samples.items():print(f"{name} : {value.mean(0)}")#########预测和评价############################
# 使用后验样本来计算预测的分布
predictive = pyro.infer.Predictive(model, posterior_samples)
y_pred = predictive(test_X0, test_X1)["obs"]mse = mean_squared_error(test_y, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(test_y, y_pred)
r2 = r2_score(test_y, y_pred)# 打印统计指标
print(f'Mean Squared Error (MSE): {mse}')
print(f'Root Mean Squared Error (RMSE): {rmse}')
print(f'Mean Absolute Error (MAE): {mae}')
print(f'R-squared (R2): {r2}')'''
torch.Size([100]) torch.Size([100])
Sample: 100%|██████████████████████████████████████████| 800/800 [04:39,  2.86it/s, step size=1.76e-02, acc. prob=0.919]beta0 : 0.06777436401333638
beta1 : 1.1556829075576296
beta2 : 0.8807723691323056
sigma : 0.01814367967435664
Mean Squared Error (MSE): 0.0002873891791732731
Root Mean Squared Error (RMSE): 0.016952556714940465
Mean Absolute Error (MAE): 0.013864298414787972
R-squared (R2): 0.8925269422276634
'''

这段代码是一个使用Pyro库进行贝叶斯线性回归的统计分析和预测的例子。Pyro是一个Python库，用于概率编程，它构建在PyTorch上。代码的目的是通过MCMC（Markov Chain Monte Carlo）方法中的NUTS（No-U-Turn Sampler）算法来估计模型参数的后验分布，并基于这些参数的后验分布进行预测。
以下是代码的详细解释：
导入必要的库：

• torch 是一个开源机器学习库，提供了多种高级功能来处理多维数组（张量）和执行自动微分。
• pyro 是一个基于PyTorch的概率编程库，它允许用户定义概率模型并进行推断。
• pyro.distributions 是Pyro中包含概率分布的模块，用于定义模型中的随机变量。
• pyro.infer 包含推断算法的模块，这里特别用到了MCMC和NUTS算法。
• pandas 是一个数据操作和分析的库，通常用于处理结构化数据。
• numpy是一个包含多种数学计算工具的库。
• sklearn.metrics 提供了一系列评估预测误差的函数。

数据加载和预处理：

• 使用pandas库从’train_set.csv’和’test_set.csv’文件中加载训练集和测试集数据。
• 将训练集和测试集中的x0和x1列转换为PyTorch张量，用于模型输入。
• 将训练集和测试集中的y列（目标变量）转换为PyTorch张量。

定义模型：

• model函数定义了一个贝叶斯线性回归模型，其中x0和x1是输入特征，y是可选的观测目标变量。
• beta0, beta1, beta2是线性模型的系数，它们被赋予了正态分布作为它们的先验。
• sigma是模型误差项的标准差，赋予了一个0到10的均匀分布作为先验。
• mean是线性回归模型的均值部分，根据系数和输入特征计算得出。
• 使用pyro.plate来指示y观测值的索引，它能够提高推断的效率。
• obs是观测数据的似然，它是以mean和sigma为参数的正态分布。

MCMC推断：

• 使用NUTS算法作为MCMC推断的内核，来估计模型参数的后验分布。
• 设置采样500次并预热（warmup）300步，预热步骤帮助算法找到合适的参数开始采样，这有助于改进采样质量。
• mcmc.run运行MCMC算法来采样模型的后验分布。
• mcmc.get_samples()获取后验分布的样本。

后验分析与预测：

• 利用后验分布的样本和pyro.infer.Predictive来预测测试集的y值。
• 计算了均方误差（MSE）、均方根误差（RMSE）、平均绝对误差（MAE）和决定系数（R^2）来评估模型的预测性能。

输出结果：

• 输出的是后验分布的平均值，这些平均值表示模型参数的估计值。
• 输出各种统计指标来评价模型的预测效果，包括MSE、RMSE、MAE和R^2。

最后，注释中的输出结果显示了训练过程中得到的参数估计和评估指标的特定值。这些值表明模型拟合得相当好，因为$R^2$值接近1，而MSE、RMSE、MAE都相对较低，表明预测误差小。

接下来，我们不使用现成的python库，直接实现贝叶斯线性回归。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import pandas as pd
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score# 使用对数似然函数
def log_likelihood(w, data):x0, x1, y = datay_pred = w[0] + w[1] * x0 + w[2] * x1# 注意我们这里是返回对数似然return np.sum(stats.norm.logpdf(y, y_pred, 1)) # 假设误差项标准差为 1# 使用对数先验分布
def log_prior(w):# 返回参数的对数先验，这里同样假设参数的先验分布为正态分布return stats.norm.logpdf(w[0], 0, 1) + stats.norm.logpdf(w[1], 0, 1) + stats.norm.logpdf(w[2], 0, 1)# 随机提案函数（用于生成新的参数采样点）
def random_proposal(w, step_size=[0.5, 0.5, 0.5]):return np.random.normal(w, step_size)# 函数: Metropolis-Hastings 算法实现
def metropolis_hastings(log_likelihood_computer, log_prior_computer, random_proposal, iterations, data, initial_values=[0,0,0], step_size=[0.5, 0.5, 0.5]):# 初始参数值w = initial_values# 保存所有采样值samples = [w]# 计算初始参数的对数先验和对数似然log_prior_current = log_prior_computer(w)log_likelihood_current = log_likelihood_computer(w, data)log_posterior_current = log_likelihood_current + log_prior_currentfor i in range(iterations):# 生成新的参数采样点w_new = random_proposal(w, step_size)# 计算新采样点的对数先验和对数似然log_prior_new = log_prior_computer(w_new)log_likelihood_new = log_likelihood_computer(w_new, data)log_posterior_new = log_likelihood_new + log_prior_new# 对数接受概率log_acceptance_probability = min(0, log_posterior_new - log_posterior_current)# 接受或拒绝新采样点accept = np.log(np.random.rand()) < log_acceptance_probabilityif accept:# 更新当前参数值w = w_new.copy()log_prior_current = log_prior_newlog_likelihood_current = log_likelihood_newlog_posterior_current = log_posterior_newsamples.append(w.copy())return np.array(samples)# 预测新数据点
def predict(posterior_samples, x0_new, x1_new):prediction_mean_list = []for item in range(len(x0_new)):# 计算每个样本的预测值predictions = posterior_samples[:, 0] + posterior_samples[:, 1] * x0_new[item] + posterior_samples[:, 2] * x1_new[item]prediction_mean = np.mean(predictions)prediction_mean_list.append(prediction_mean)return np.array(prediction_mean_list) # 加载数据
data = pd.read_csv('train_set.csv')  
test_data = pd.read_csv('test_set.csv')  
# 模拟数据
np.random.seed(42) # 设置随机种子
true_w0 = 1
true_w1 = 3
true_w2 = -2
x0_data = data['x0'].values
x1_data = data['x0'].values
y_data = data['y'].values
test_X0 = test_data['x0'].values
test_X1 = test_data['x1'].values
test_y = test_data['y'].values
print(x0_data.shape,x1_data.shape,y_data.shape)# 运行 Metropolis-Hastings 算法
iterations = 10000
samples = metropolis_hastings(log_likelihood, log_prior, random_proposal, iterations, data=(x0_data, x1_data, y_data))# 查看参数的后验分布
burn_in = 8000 # 预热期间的样本数
posterior_samples = samples[burn_in+1:]
posterior_samples_mean = np.mean(posterior_samples,axis=0)
print(posterior_samples_mean)
y_pred = predict(posterior_samples,test_X0,test_X1)mse = mean_squared_error(test_y, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(test_y, y_pred)
r2 = r2_score(test_y, y_pred)# 打印统计指标
print(f'Mean Squared Error (MSE): {mse}')
print(f'Root Mean Squared Error (RMSE): {rmse}')
print(f'Mean Absolute Error (MAE): {mae}')
print(f'R-squared (R2): {r2}')

这段代码是一个 Python 脚本，它实现了贝叶斯线性回归模型，并通过 Metropolis-Hastings 算法（一种马尔可夫链蒙特卡洛方法）来进行参数估计。脚本中用到了 NumPy、SciPy、Pandas、和 scikit-learn 这些库。
以下是代码的具体步骤和组成部分的解释：

1. 导入必要的库：

• numpy 用于数组和矩阵运算；
• matplotlib.pyplot 通常用于绘制图形，但在这段代码中并未使用；
• scipy.stats 用于统计操作，这里主要用于正态分布的概率密度函数；
• pandas 用于数据操作和分析，特别是用于读取 CSV 文件和数据处理；
• sklearn.metrics 提供了用于评估预测误差的函数。

2. 定义对数似然函数 log_likelihood，它接收参数 w（模型的权重）和数据 data，并根据正态分布返回给定数据下模型参数 w 的对数似然值。
3. 定义对数先验分布 log_prior，它基于模型参数 w 的正态分布先验，返回参数的对数先验概率。
4. 定义随机提案函数 random_proposal，它用于在 Metropolis-Hastings 算法中生成新的参数样本，通过添加正态分布噪声实现。
5. 实现 Metropolis-Hastings 算法的 metropolis_hastings 函数。这个函数执行以下步骤：

• 初始化参数值；
• 在每次迭代中生成新的参数样本；
• 根据新旧样本的对数后验概率差计算接受率；
• 用一个随机数与接受率比较来确定是否接受新样本；
• 重复上述过程直到完成指定次数的迭代；
• 返回所有的样本。

6. 定义 predict 函数，它使用后验样本来预测新数据点的响应值的均值。
7. 加载训练和测试数据集。
8. 设置随机种子以确保结果的可重复性。
9. 从训练数据中提取输入特征 x0_data, x1_data 和目标变量 y_data。
10. 运行 Metropolis-Hastings 算法以获得模型参数的样本。
11. 定义并移除后验样本中的“烧入”期，这是模拟的开始部分，可能还没有收敛。
12. 计算去除“烧入”样本后的后验样本的均值。
13. 使用训练好的模型和后验样本对测试数据集进行预测。
14. 计算并打印出模型预测的统计指标，包括均方误差（MSE）、均方根误差（RMSE）、平均绝对误差（MAE）和 R-squared 值。

代码的主要目的是展示如何通过贝叶斯方法和 MCMC 算法来进行线性回归模型的参数估计，并对新数据进行预测，最后评估模型的预测性能。这是机器学习和统计模型中的一种常见方法，特别是在涉及到不确定性估计和概率编程时。