建设网站外国人可搜到,扬州抖音推广公司,无锡网站排名团队,做网站实训报告总结整体概述
频率主义#xff08;Frequentist#xff09;线性回归和贝叶斯#xff08;Bayesian#xff09;线性回归是统计学中用于数据分析和预测的两种主要方法#xff0c;特别是在建模关于因变量和自变量之间线性关系的上下文中。尽管两种方法都用于线性回归分析#xff…整体概述
频率主义Frequentist线性回归和贝叶斯Bayesian线性回归是统计学中用于数据分析和预测的两种主要方法特别是在建模关于因变量和自变量之间线性关系的上下文中。尽管两种方法都用于线性回归分析但它们在解释参数、估计方法和结果的不确定性方面有所不同。
频率主义线性回归
频率主义线性回归建立在经典统计学的原则上尤其是参数估计的频率主义概念。在频率主义视角下模型参数被认为是固定的未知常数我们的目标是基于样本数据去估计这些常数。频率主义线性回归的目标是找到一个线性模型来描述因变量Y和一组自变量X之间的关系 Y β 0 β 1 X 1 β 2 X 2 … β p X p ϵ Y \beta_0 \beta_1 X_1 \beta_2 X_2 \ldots \beta_p X_p \epsilon Yβ0β1X1β2X2…βpXpϵ 其中 β 0 , β 1 , . . . , β p \beta_0, \beta_1, ..., \beta_p β0,β1,...,βp是模型参数而 ϵ \epsilon ϵ是误差项通常假设为独立同分布的正态分布。 参数估计常用最小二乘法OLS目标是最小化误差项的平方和。频率主义的置信区间和假设检验如t检验、F检验用于评估模型参数的不确定性和显著性。
贝叶斯线性回归
贝叶斯线性回归遵循贝叶斯统计的原则模型的参数被视为随机变量而非固定未知常数。这种方法将先验知识或者信念和样本数据结合起来以估计参数。贝叶斯回归的基本形式与频率主义线性回归相同但在估计参数时采用的是贝叶斯定理 P ( Θ ∣ D a t a ) P ( D a t a ∣ Θ ) ⋅ P ( Θ ) P ( D a t a ) P(\Theta | Data) \frac{P(Data | \Theta) \cdot P(\Theta)}{P(Data)} P(Θ∣Data)P(Data)P(Data∣Θ)⋅P(Θ) 这里的 P ( Θ ∣ D a t a ) P(\Theta | Data) P(Θ∣Data) 是给定数据后参数 Θ \Theta Θ 的后验概率 P ( D a t a ∣ Θ ) P(Data | \Theta) P(Data∣Θ)是给定参数下数据的似然 P ( Θ ) P(\Theta) P(Θ) 是参数的先验概率而 P ( D a t a ) P(Data) P(Data) 是证据或边缘概率作为归一化常数。 在贝叶斯回归中不是寻找一组固定的参数估计而是利用计算或模拟方法如马尔可夫链蒙特卡洛MCMC来得到参数的后验分布。这允许我们直接从后验分布中得到参数估计的不确定性表现为分布的形状也能通过分布的特征如均值、中位数、百分位数区间来总结参数的估计。
比较
参数解释频率主义视角下参数是固定的而在贝叶斯视角下参数是具有分布的随机变量。先验知识贝叶斯方法可以整合先验知识频率主义方法则通常不考虑先验知识。不确定性评估贝叶斯方法通过后验分布直接提供参数的不确定性评估而频率主义方法使用置信区间和假设检验来进行不确定性评估。计算复杂性贝叶斯方法通常计算上更为复杂需要使用MCMC等模拟技术而频率主义方法如OLS在计算上更简单直接。
每种方法都有其优势和局限性应用哪种方法取决于具体的问题、数据和分析目标。
频率主义线性回归
频率主义线性回归是一种统计方法用于估计和推断一个或多个自变量解释变量与因变量响应变量之间的线性关系。在频率主义框架中参数被认为是固定的且数据是随机的。这一观点反映了一个信念即在重复实验的条件下某事件的频率趋近于其真实概率。
线性回归模型
线性回归模型可以表示为 Y X β ϵ Y X\beta \epsilon YXβϵ 其中 Y Y Y 是一个 n × 1 n \times 1 n×1的向量表示因变量的观测值 X X X 是一个 n × p n \times p n×p 的矩阵其中 n n n是样本大小 p p p是参数的数量包括截距表示自变量 β \beta β是一个 p × 1 p \times 1 p×1 的向量表示未知的回归系数 ϵ \epsilon ϵ 是一个 n × 1 n \times 1 n×1的向量表示误差项假设 ϵ ∼ N ( 0 , σ 2 I ) \epsilon \sim N(0, \sigma^2I) ϵ∼N(0,σ2I)即误差项独立同分布且呈正态分布。
最小二乘估计
在频率主义线性回归中通常使用最小二乘法Ordinary Least Squares, OLS来估计参数 β \beta β。最小二乘估计旨在最小化误差项的平方和 S ( β ) ( Y − X β ) ′ ( Y − X β ) S(\beta) (Y - X\beta)(Y - X\beta) S(β)(Y−Xβ)′(Y−Xβ) 通过对 S ( β ) S(\beta) S(β)求导并设其为零可以得到 β \beta β 的OLS估计 β ^ ( X ′ X ) − 1 X ′ Y \hat{\beta} (XX)^{-1}XY β^(X′X)−1X′Y 其中 β ^ \hat{\beta} β^是参数的估计值。
置信区间
置信区间是用来估计模型参数不确定性的一种方法。对于每个回归系数 β i \beta_i βi我们可以计算一个置信区间该区间以一定的置信水平如95%覆盖真实参数值。置信区间基于回归系数的抽样分布通常假设为正态分布可以表示为 β i ^ ± t α / 2 , n − p ⋅ S E ( β i ^ ) \hat{\beta_i} \pm t_{\alpha/2, n-p} \cdot SE(\hat{\beta_i}) βi^±tα/2,n−p⋅SE(βi^) 其中 t α / 2 , n − p t_{\alpha/2, n-p} tα/2,n−p 是自由度为 n − p n-p n−p 的t分布的临界值 S E ( β i ^ ) SE(\hat{\beta_i}) SE(βi^) 是估计系数 β i ^ \hat{\beta_i} βi^的标准误差。
假设检验
假设检验用于评估模型系数是否显著不同于零或其他指定的值。最常用的是t检验用于单个系数以及F检验用于模型中的多个系数。t检验的零假设是系数等于零 ( H 0 : β i 0 ) (H_0: \beta_i 0) (H0:βi0)备择假设是系数不等于零 ( H 1 : β i ≠ 0 ) (H_1: \beta_i \neq 0) (H1:βi0)。 t统计量计算为 t β i ^ − 0 S E ( β i ^ ) t \frac{\hat{\beta_i} - 0}{SE(\hat{\beta_i})} tSE(βi^)βi^−0 如果 (t) 的绝对值大于t分布的临界值我们拒绝零假设认为这个系数显著不为零。 F检验通常用于同时检验多个系数是否显著。它比较模型中所有解释变量与只有截距项的简化模型的拟合度。F统计量是比例反映了两个模型中未解释变化的比例减少 F ( S S E 0 − S S E 1 ) / ( p − 1 ) S S E 1 / ( n − p ) F \frac{(SSE_0 - SSE_1)/(p-1)}{SSE_1/(n-p)} FSSE1/(n−p)(SSE0−SSE1)/(p−1) 其中 S S E 0 SSE_0 SSE0 是没有解释变量的模型的误差平方和 S S E 1 SSE_1 SSE1 是完整模型的误差平方和 p p p 是完整模型中参数的数量 n n n 是样本大小。
如果F统计量的值大于F分布的临界值我们拒绝零假设认为模型中至少有一个系数显著不为零。 通过这些方法频率主义线性回归不仅为我们提供了对数据的最佳线性拟合而且还允许我们进行推断评估模型系数的显著性以及估计系数的可靠性。
贝叶斯线性回归
贝叶斯线性回归通过结合先验分布和数据的似然来更新关于回归参数的信念得到参数的后验分布。由于后验分布通常难以解析计算尤其是在参数空间较大或模型复杂时马尔可夫链蒙特卡洛MCMC方法被广泛用于生成后验分布的近似样本。这些样本可以用于参数估计、预测以及计算可信区间。 以下是贝叶斯线性回归估计和预测的详细步骤包括使用MCMC方法和计算可信区间的说明。
模型设定
首先设定一个贝叶斯线性回归模型 Y X β ϵ , ϵ ∼ N ( 0 , σ 2 I ) Y X\beta \epsilon, \quad \epsilon \sim N(0, \sigma^2I) YXβϵ,ϵ∼N(0,σ2I) Y Y Y是观测的响应变量 X X X 是设计矩阵包含自变量 β \beta β 是回归系数需要估计 ϵ \epsilon ϵ是误差项假设服从以0为均值、 σ 2 \sigma^2 σ2为方差的正态分布。
先验分布设定
接着给定参数 β \beta β 和 σ 2 \sigma^2 σ2 的先验分布。这些先验可以是无信息先验例如正态先验或者均匀先验也可以是包含先前知识的信息先验。
后验分布
根据贝叶斯定理后验分布与似然函数和先验分布成正比 p ( β , σ 2 ∣ Y , X ) ∝ p ( Y ∣ X , β , σ 2 ) p ( β ) p ( σ 2 ) p(\beta, \sigma^2 | Y, X) \propto p(Y | X, \beta, \sigma^2) p(\beta) p(\sigma^2) p(β,σ2∣Y,X)∝p(Y∣X,β,σ2)p(β)p(σ2) 其中 p ( β , σ 2 ∣ Y , X ) p(\beta, \sigma^2 | Y, X) p(β,σ2∣Y,X) 是后验分布 p ( Y ∣ X , β , σ 2 ) p(Y | X, \beta, \sigma^2) p(Y∣X,β,σ2) 是似然函数 p ( β ) p(\beta) p(β) 和 p ( σ 2 ) p(\sigma^2) p(σ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 x_0,x_1 x0,x1值预测 y y 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(fIntercept: {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(fMean Squared Error (MSE): {mse})
print(fRoot Mean Squared Error (RMSE): {rmse})
print(fMean Absolute Error (MAE): {mae})
print(fR-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中的基础功能来计算线性回归模型的参数。这可以通过计算最小二乘估计来实现即通过解析解找到系数β的估计值使得残差平方和最小。 线性回归模型的解析解可通过以下公式计算得到 β ( X T X ) − 1 X T y \beta (X^TX)^{-1}X^Ty β(XTX)−1XTy 这里的 β \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(fIntercept (beta_0): {beta_hat[0]})
print(fCoefficient of x0 (beta_1): {beta_hat[1]})
print(fCoefficient 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(fMean Squared Error (MSE): {mse})
print(fRoot Mean Squared Error (RMSE): {rmse})
print(fMean Absolute Error (MAE): {mae})
print(fR-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, yNone):# 定义先验分布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), obsy)return mean# 运行Hamiltonian Monte Carlo (HMC)算法NUTS是HMC的一个优化版本
nuts_kernel NUTS(model)
mcmc MCMC(nuts_kernel, num_samples500, warmup_steps300)
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(fMean Squared Error (MSE): {mse})
print(fRoot Mean Squared Error (RMSE): {rmse})
print(fMean Absolute Error (MAE): {mae})
print(fR-squared (R2): {r2})
torch.Size([100]) torch.Size([100])
Sample: 100%|██████████████████████████████████████████| 800/800 [04:39, 2.86it/s, step size1.76e-02, acc. prob0.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上。代码的目的是通过MCMCMarkov Chain Monte Carlo方法中的NUTSNo-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次并预热warmup300步预热步骤帮助算法找到合适的参数开始采样这有助于改进采样质量。mcmc.run运行MCMC算法来采样模型的后验分布。mcmc.get_samples()获取后验分布的样本。
后验分析与预测
利用后验分布的样本和pyro.infer.Predictive来预测测试集的y值。计算了均方误差MSE、均方根误差RMSE、平均绝对误差MAE和决定系数R^2来评估模型的预测性能。
输出结果
输出的是后验分布的平均值这些平均值表示模型参数的估计值。输出各种统计指标来评价模型的预测效果包括MSE、RMSE、MAE和R^2。
最后注释中的输出结果显示了训练过程中得到的参数估计和评估指标的特定值。这些值表明模型拟合得相当好因为 R 2 R^2 R2值接近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_in1:]
posterior_samples_mean np.mean(posterior_samples,axis0)
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(fMean Squared Error (MSE): {mse})
print(fRoot Mean Squared Error (RMSE): {rmse})
print(fMean Absolute Error (MAE): {mae})
print(fR-squared (R2): {r2})这段代码是一个 Python 脚本它实现了贝叶斯线性回归模型并通过 Metropolis-Hastings 算法一种马尔可夫链蒙特卡洛方法来进行参数估计。脚本中用到了 NumPy、SciPy、Pandas、和 scikit-learn 这些库。 以下是代码的具体步骤和组成部分的解释
导入必要的库
numpy 用于数组和矩阵运算matplotlib.pyplot 通常用于绘制图形但在这段代码中并未使用scipy.stats 用于统计操作这里主要用于正态分布的概率密度函数pandas 用于数据操作和分析特别是用于读取 CSV 文件和数据处理sklearn.metrics 提供了用于评估预测误差的函数。
定义对数似然函数 log_likelihood它接收参数 w模型的权重和数据 data并根据正态分布返回给定数据下模型参数 w 的对数似然值。定义对数先验分布 log_prior它基于模型参数 w 的正态分布先验返回参数的对数先验概率。定义随机提案函数 random_proposal它用于在 Metropolis-Hastings 算法中生成新的参数样本通过添加正态分布噪声实现。实现 Metropolis-Hastings 算法的 metropolis_hastings 函数。这个函数执行以下步骤
初始化参数值在每次迭代中生成新的参数样本根据新旧样本的对数后验概率差计算接受率用一个随机数与接受率比较来确定是否接受新样本重复上述过程直到完成指定次数的迭代返回所有的样本。
定义 predict 函数它使用后验样本来预测新数据点的响应值的均值。加载训练和测试数据集。设置随机种子以确保结果的可重复性。从训练数据中提取输入特征 x0_data, x1_data 和目标变量 y_data。运行 Metropolis-Hastings 算法以获得模型参数的样本。定义并移除后验样本中的“烧入”期这是模拟的开始部分可能还没有收敛。计算去除“烧入”样本后的后验样本的均值。使用训练好的模型和后验样本对测试数据集进行预测。计算并打印出模型预测的统计指标包括均方误差MSE、均方根误差RMSE、平均绝对误差MAE和 R-squared 值。
代码的主要目的是展示如何通过贝叶斯方法和 MCMC 算法来进行线性回归模型的参数估计并对新数据进行预测最后评估模型的预测性能。这是机器学习和统计模型中的一种常见方法特别是在涉及到不确定性估计和概率编程时。
