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

整体概述

频率主义（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​+β1​X1​+β2​X2​+…+βp​Xp​+ϵ
其中， β 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 Y=Xβ+ϵ
其中，

  • 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} = (X'X)^{-1}X'Y β^​=(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​:βi​=0)，备择假设是系数不等于零 ( H 1 : β i ≠ 0 ) (H_1: \beta_i \neq 0) (H1​:βi​=0)。
t统计量计算为：
t = β i ^ − 0 S E ( β i ^ ) t = \frac{\hat{\beta_i} - 0}{SE(\hat{\beta_i})} t=SE(β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)} F=SSE1​/(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) Y=Xβ+ϵ,ϵ∼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(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中的基础功能来计算线性回归模型的参数。这可以通过计算最小二乘估计来实现，即通过解析解找到系数β的估计值，使得残差平方和最小。
线性回归模型的解析解可通过以下公式计算得到：
β = ( 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(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 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_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 提供了用于评估预测误差的函数。
  1. 定义对数似然函数 log_likelihood，它接收参数 w（模型的权重）和数据 data，并根据正态分布返回给定数据下模型参数 w 的对数似然值。
  2. 定义对数先验分布 log_prior，它基于模型参数 w 的正态分布先验，返回参数的对数先验概率。
  3. 定义随机提案函数 random_proposal，它用于在 Metropolis-Hastings 算法中生成新的参数样本，通过添加正态分布噪声实现。
  4. 实现 Metropolis-Hastings 算法的 metropolis_hastings 函数。这个函数执行以下步骤：
  • 初始化参数值；
  • 在每次迭代中生成新的参数样本；
  • 根据新旧样本的对数后验概率差计算接受率；
  • 用一个随机数与接受率比较来确定是否接受新样本；
  • 重复上述过程直到完成指定次数的迭代；
  • 返回所有的样本。
  1. 定义 predict 函数，它使用后验样本来预测新数据点的响应值的均值。
  2. 加载训练和测试数据集。
  3. 设置随机种子以确保结果的可重复性。
  4. 从训练数据中提取输入特征 x0_data, x1_data 和目标变量 y_data。
  5. 运行 Metropolis-Hastings 算法以获得模型参数的样本。
  6. 定义并移除后验样本中的“烧入”期，这是模拟的开始部分，可能还没有收敛。
  7. 计算去除“烧入”样本后的后验样本的均值。
  8. 使用训练好的模型和后验样本对测试数据集进行预测。
  9. 计算并打印出模型预测的统计指标，包括均方误差（MSE）、均方根误差（RMSE）、平均绝对误差（MAE）和 R-squared 值。

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

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mzph.cn/news/694478.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

【LeetCode】746. 使用最小花费爬楼梯（简单）——代码随想录算法训练营Day38

【LeetCode】746. 使用最小花费爬楼梯（简单）——代码随想录算法训练营Day38

题目链接&#xff1a;746. 使用最小花费爬楼梯 题目描述 给你一个整数数组 cost &#xff0c;其中 cost[i] 是从楼梯第 i 个台阶向上爬需要支付的费用。一旦你支付此费用&#xff0c;即可选择向上爬一个或者两个台阶。 你可以选择从下标为 0 或下标为 1 的台阶开始爬楼梯。 …
阅读更多...
Git 客户端可视化工具tortoisegit

Git 客户端可视化工具tortoisegit

Git 使用教程 git一点通 (kdocs.cn) 二、Git 客户端可视化工具-推荐 1.常用工具 tortoisegit 官网 https://tortoisegit.org/ 推荐 sourcetree 官网 https://www.sourcetreeapp.com/ 2.tortoisegit安装 2.1 下载安装包 2.2 下载语言包 2.3 安装 2.4 安装语言包 5.使用 5.1 新建…
阅读更多...
C++ CRTP设计范式

C++ CRTP设计范式

CRTP&#xff08;Curiously Recurring Template Pattern&#xff09;奇异递归模板范式是一个相对少有人知的C设计范式&#xff0c;它可以实现基类指针调用派生类的函数&#xff0c;从而实现另类多态。 如&#xff1a; #include <iostream> // 基类模板&#xff0c;接受…
阅读更多...
Sora 一款文本转视频模型

Sora 一款文本转视频模型

**Sora** 是一个由美国人工智能研究机构 **OpenAI** 开发的 AI 视频模型。让我们一起了解一下这个令人兴奋的项目&#xff1a; 1. **名称和含义**&#xff1a; - 在日语中&#xff0c;**Sora** 是“天空”的意思&#xff0c;引申含义还有“自由”。 - **Sora** 的官方介绍页上展…
阅读更多...
js之事件循环

js之事件循环

JavaScript的事件循环是它的并发模型的核心部分&#xff0c;使得JavaScript能够在单线程中处理异步操作。事件循环允许JavaScript在执行代码时&#xff0c;同时进行非阻塞的I/O操作&#xff08;如网络请求、文件操作等&#xff09;。这个概念对于理解如何高效地构建交互式Web应…
阅读更多...
本地模拟发送、接收RabbitMQ数据

本地模拟发送、接收RabbitMQ数据

文章目录 前言一、相关文章二、相关代码1.模拟的 Channel 类2.接收消息3.模拟推送MQ数据前言 日常开发中，当线上RabbitMQ坏境还没准备好时，可在本地模拟发送、接收消息 一、相关文章 Docker安装RabbitMQ 【SpringCloud】整合RabbitMQ六大模式应用（入门到精通） Spring R…
阅读更多...
Spring学习笔记（三）--Spring中的Bean的管理

Spring学习笔记（三）--Spring中的Bean的管理

一、什么是Bean Bean是注册到Spring容器中的Java类&#xff0c;控制反转和依赖注入都是通过Bean实现的&#xff0c;任何一个Java类都可以是一个Bean。Bean由Spring进行管理&#xff0c;可以通过xml文件对bean进行配置和管理。 二、BeanFactory接口和ApplicationContext接口&a…
阅读更多...
利用Python实现科学式占卜

利用Python实现科学式占卜

一直以来&#xff0c;中式占卜都是基于算命先生手工实现&#xff0c;程序繁琐&#xff08;往往需要沐浴、计算天时、静心等等流程&#xff09;。准备工作复杂&#xff08;通常需要铜钱等道具&#xff09;&#xff0c;计算方法复杂&#xff0c;需要纯手工计算二进制并转换为最终…
阅读更多...
2023年12月 Python（六级）真题解析#中国电子学会#全国青少年软件编程等级考试

2023年12月 Python（六级）真题解析#中国电子学会#全国青少年软件编程等级考试

Python等级考试（1~6级）全部真题・点这里 一、单选题（共25题，共50分） 第1题 运行以下程序，输出的结果是？（ ） class A():def __init__(self,x):self.x=x
阅读更多...
手机连接电脑后资源管理器无法识别(识别设备但无法访问文件)

手机连接电脑后资源管理器无法识别(识别设备但无法访问文件)

问题描述 小米8刷了pixel experience系统,今天用电脑连接后无法访问手机文件,但是手机选择了usb传输模式为文件传输 解决办法 在设备和打印机页面中右键选择属性 点击改变设置 卸载驱动,注意勾选删除设备的驱动程序软件 卸载后重新连接手机,电脑弹出希望对设备进行什么操作时…
阅读更多...
Code-Audit(代码审计)习题记录4-5

Code-Audit(代码审计)习题记录4-5

4、习题4 题目内容如下&#xff1a; <?php error_reporting(0); show_source(__FILE__); $a $_REQUEST[hello]; eval("var_dump($a);"); 函数解释 $REQUEST — HTTP Request 变量&#xff0c;默认情况下包含了 [$GET]&#xff0c;[$POST] 和 [$COOKIE]的数…
阅读更多...
HTML知识点

HTML知识点

HTML 【一】HTML简介 【1】什么是HTML HTML是一种用于创建网页结构和内容的超文本标记语言&#xff0c;它是构建网页的基础。为了让浏览器正确渲染页面&#xff0c;我们必须遵循HTML的语法规则。浏览器在解析网页时会将HTML代码转换为可视化的页面&#xff0c;所以我们在浏览…
阅读更多...
python接口测试实践：参数化测试、数据驱动测试和断言的使用

python接口测试实践：参数化测试、数据驱动测试和断言的使用

在Python接口测试实践中&#xff0c;参数化测试、数据驱动测试和断言是常用的技术手段。 参数化测试 参数化测试是指将测试用例中的某些部分&#xff08;如输入数据或配置&#xff09;作为参数传递给测试函数&#xff0c;以便于复用和减少代码重复。 例如&#xff0c;使用uni…
阅读更多...
2024年最新腾讯云轻量4核8G12M服务器CPU内存性价比如何？

2024年最新腾讯云轻量4核8G12M服务器CPU内存性价比如何？

4核8G服务器支持多少人同时在线访问&#xff1f;阿腾云的4核8G服务器可以支持20个访客同时访问&#xff0c;关于4核8G服务器承载量并发数qps计算测评&#xff0c;云服务器上运行程序效率不同支持人数在线人数不同&#xff0c;公网带宽也是影响4核8G服务器并发数的一大因素&…
阅读更多...
设计模式--组合模式（Composite Pattern）

设计模式--组合模式（Composite Pattern）

组合模式（Composite Pattern）是一种结构型设计模式，它允许你将对象组合成树形结构，并且能像使用独立对象一样使用它们。 组合模式主要包含以下几个角色： Component：这是组合中对象声明接口，在适当的情况下，实现所有类共有接口的默认行为。声明一个接口用于访问和管理C…
阅读更多...
oracle数据库事务的四大特性与隔离级别与游标

oracle数据库事务的四大特性与隔离级别与游标

数据库事务的四大特性: 这里提到了 ACID 四个特性&#xff0c;分别是&#xff1a; A&#xff08;Atomicity&#xff09;&#xff1a; 原子性&#xff0c;确保事务中的所有操作要么全部执行成功&#xff0c;要么全部不执行&#xff0c;不存在部分执行的情况。 C&#xff08;…
阅读更多...
时序分析深入必学的时序模型详细讲解

时序分析深入必学的时序模型详细讲解

目录 一、前言 二、时序建模 2.1 反相器 2.2 线性时序模型 2.3 非线性时序模型 三、时序模块 3.1 组合单元 3.2 时序单元 3.3 同步检查&#xff1a;setup和hold 3.4 异步检查&#xff1a;recovery和removal 3.5 脉冲宽度检查Pulse width check 3.6 传输时延 3.7 …
阅读更多...
微信小程序配置文件

微信小程序配置文件

目录 小程序配置文件 1. 配置文件介绍 2. 全局配置 2.1 pages 2.2 window 2.3 tabBar 3. 页面配置 4. 项目配置文件 5. 支持使用 sass/less 6. sitemap.json 小程序配置文件 1. 配置文件介绍 JSON是一种轻量级的数据格式&#xff0c;常用于前后端数据的交互&#xf…
阅读更多...
LDRA Testbed软件静态分析_软件质量度量

LDRA Testbed软件静态分析_软件质量度量

系列文章目录 LDRA Testbed软件静态分析_操作指南 LDRA Testbed软件静态分析_自动提取静态分析数据生成文档 LDRA Testbed软件静态分析_Jenkins持续集成_(1)自动进行静态分析的环境搭建 LDRA Testbed软件静态分析_Jenkins持续集成_(2)配置邮件自动发送静态分析结果 LDRA Testb…
阅读更多...
python数据类型-集合set

python数据类型-集合set

1 集合&#xff08;set&#xff09;的定义 1.1 集合是一个无序且不重复元素的序列&#xff1a; 1&#xff09;无序&#xff1a;存储顺序和添加的顺序不一定相同&#xff0c;不支持索引、切片 2&#xff09;元素不重复&#xff1a;当添加重复元素时&#xff0c;集合会自动去重…
阅读更多...
最新文章

Copyright @ 2022~2023