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

已于 2024-02-20 22:50:18 修改
阅读量1.2k 收藏 23
点赞数 17
分类专栏： 深度学习自学之路 文章标签： 线性回归 算法 回归
于 2024-02-20 15:02:23 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_39699362/article/details/136190610
版权
博客主要介绍频率主义和贝叶斯线性回归两种方法。频率主义线性回归基于经典统计学，用最小二乘法估计参数，通过置信区间和假设检验评估；贝叶斯线性回归遵循贝叶斯统计，结合先验知识，用MCMC采样得到后验分布。还给出示例对比两者应用及性能。

摘要生成于 C知道 ，由 DeepSeek-R1 满血版支持， 前往体验 >

整体概述

频率主义（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_y

print(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 @ beta
    return 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 = data
    y_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_current

    for 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_probability
        if accept:
            # 更新当前参数值
            w = w_new.copy()
            log_prior_current = log_prior_new
            log_likelihood_current = log_likelihood_new
            log_posterior_current = log_posterior_new

        samples.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 算法来进行线性回归模型的参数估计，并对新数据进行预测，最后评估模型的预测性能。这是机器学习和统计模型中的一种常见方法，特别是在涉及到不确定性估计和概率编程时。

贝叶斯线性回归 - Bayes
11-24
介绍了极大似然、最大后验、贝叶斯估计以及贝叶斯线性回归方法。 介绍了极大似然、最大后验、贝叶斯估计以及贝叶斯线性回归方法。 介绍了极大似然、最大后验、贝叶斯估计以及贝叶斯线性回归方法。
贝叶斯线性回归(Bayesian Linear Regression)
热门推荐
DarkScope从这里开始
01-31 4万+
mlapp看到了第七章，跳了第六章没看，第七章主要是讲线性回归的，前面首先是最朴素的线性回归，接着是ridge线性回归(其实就是带惩罚参数的回归)，然后就是本文要总结的贝叶斯线性回归。把prml和mlapp综合起来看了，效果还不错，有些东西互有详略，可以互做补充。 1.引入      贝叶斯线性回归的引入主要是在最大似然估计中很难决定模型的复杂程度，ridge回归加入的惩罚参数其实也是解决这个
【回归算法解析系列11】贝叶斯回归（Bayesian Regression）
最新发布
Is code
03-19 1007
贝叶斯回归通过概率建模的方式量化不确定性，在小样本和动态场景中表现出色。MCMC和变分推断作为近似后验的两大主流方法，各有其优缺点和适用场景。预测区间的提供为风险评估提供了直观的支持，使我们能够更全面地理解模型的预测结果。
【bayes】贝叶斯学派和频率学派
weixin_34205826的博客
04-08 677
一、 “探测仪，如果我问一个贝叶斯学派的统计学家如果……”“[掷]我是一个中微子探测仪，不是迷宫守卫。老实说，你是不是脑子坏掉了。”“[掷]...yes” 迷宫守卫的梗：说迷宫里有2条路，分别通向目的地和陷阱，路口各有一个守卫，一个只说真话一个只说假话，他们都知道路后面是什么以及彼此说话的真假，这时你只能选择其中一个人，问一个问题，如何确保走对路。 http://www.xkcd.com...
频率派机器学习——线性回归
weixin_47206104的博客
08-28 231
频率派机器学习——线性回归 线性回归的基本模型f(w,b)=wTx+bf(w,b)=w^Tx+bf(w,b)=wTx+b，线性回归有三要数： 线性 属性线性 打破→\underrightarrow{打破}打破​ 即将属性改为非线性，即存在xix^ixi (i>1)(i>1)(i>1)。则为特征转换（多项式回归） 全局线性 打破→\underrightarrow{打破}打破​ 即将全局线性改为非线性。例如神经网络中的激活函数，使其输出为非线性，则为线性分类 系数线性 打
从线性回归到贝叶斯线性回归(Bayesian Linear Regression)
David's Tweet
06-30 1950
参考： 《回字的四种写法——从线性回归到贝叶斯线性回归》 《贝叶斯线性回归(Bayesian Linear Regression)》 《透彻理解协方差矩阵》 本文顺序 一、回忆线性回归 线性回归用最小二乘法，转换为极大似然估计求解参数W，但这很容易导致过拟合，由此引入了带正则化的最小二乘法（可证明等价于最大后验概率） 二、什么是贝叶斯回归？ ...
贝叶斯线性回归模型的推导与实现
### 1.1 贝叶斯线性回归模型概述 贝叶斯线性回归模型是一种基于贝叶斯统计方法的回归分析模型。与传统的线性回归模型相比，贝叶斯线性回归模型不仅考虑了样本数据的观测误差，还引入了参数的不确定性。通过使用...
最小二乘法与逻辑斯谛回归：线性模型基石
首先，最小二乘法（Least Squares Method）是本章的核心内容，它是通过最小化样本数据点到预测直线的平均欧几里得距离来建立线性回归模型。这种方法在寻找最优拟合直线时，确保了模型的简单性和有效性。 接着，极大...
第三章 线性模型1
08-08
在频率主义观点下，参数被视为固定的未知量，通过最大化似然函数来确定；而在贝叶斯学派中，参数被视为随机变量，有自身的先验分布，后验分布可通过贝叶斯公式计算。 3. **Sigmoid函数**：Sigmoid函数，又称对数...
Python贝叶斯框架统计推断实战：线性/逻辑回归与混合模型
1. **线性回归模型**：利用贝叶斯方法进行线性回归，可以引入正则化项（如Laplace先验或Gaussian先验），通过MCMC采样估计回归系数和噪声水平的后验分布。 2. **逻辑回归模型**：用于分类问题，特别是二分类问题。...
R语言贝叶斯MCMC：GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例...
拓端研究室TRL
11-20 192
原文链接：http://tecdat.cn/?p=23236在频率学派中，观察样本是随机的，而参数是固定的、未知的数量（点击文末“阅读原文”获取完整代码数据）。相关视频什么是频率学派？概率被解释为一个随机过程的许多观测的预期频率。有一种想法是 "真实的"，例如，在预测鱼的生活环境时，盐度和温度之间的相互作用有一个回归系数？什么是贝叶斯学派？在贝叶斯方法中，概率被解释为对信念的主观衡量。所有的变量-...
Bayesian Linear Regression：这是一组做贝叶斯线性回归的MATLAB函数-matlab开发
05-30
这是一组用于进行贝叶斯线性回归的 MATLAB 函数。 派生也包括在内。
贝叶斯线性回归
Suchunsv的博客
01-07 235
贝叶斯线性回归 输入 Data : {(xi,yi)}i=1N,xi∈Rp,yi∈R\{(x_i, y_i)\}_{i=1}^N,x_i\in R^p,y_i\in R{(xi​,yi​)}i=1N​,xi​∈Rp,yi​∈R Inference 首先要清楚贝叶斯线性回归和频率派的区别在于，对y=wTxy = w^Txy=wTx，贝叶斯线性回归是假设w是一个分布而不是具体的常数的值。 其次...
贝叶斯线性回归（贝叶斯学派）
qq_36013249的博客
05-06 2321
前言 极大似然线性回归中，我们使用线性高斯判别模型，似然函数为： p(Y∣X,θ,β)=∏i=1nN(yi∣θTΦ(xi),β−1)p(Y|X,\theta,\beta)=\prod_{i=1}^n\mathcal{N}(y_i|\theta^T\Phi (x_i),\beta^{-1})p(Y∣X,θ,β)=i=1∏n​N(yi​∣θTΦ(xi​),β−1) 其中Φ(xi)\Phi (x_i)...
线性回归和贝叶斯的线性回归
Anarion的博客
06-06 1633
线性回归和贝叶斯的线性回归 问题是什么 这个h2标题令人疑惑，不就是线性回归么，高中就会啦。 y=θTx+b y=\theta^Tx+b y=θTx+b 或者也写成这样： y=θTx,x0=1 y=\theta^Tx,x_0=1 y=θTx,x0​=1 用一个xxx的分量代表偏置项bbb，为了方便表述，以下都使用这个记号。 线性回归的目标是，用线性的函数描述任何数据，使得误差最小。误差就是这样： L(i)=∣∣y(i)−(θTx(i))∣∣22 L^{(i)}=||y^{(i)}-(\theta^Tx^{(
贝叶斯线性回归（Bayesian Linear Regression）
u013250861的博客
08-30 910
基于上面的讨论，这里就可以引出本文的核心内容：贝叶斯线性回归。贝叶斯线性回归不仅可以解决极大似然估计中存在的过拟合的问题。它对数据样本的利用率是100%，仅仅使用训练样本就可以有效而准确的确定模型的复杂度。在极大似然估计线性回归中我们把参数看成是一个未知的固定值，而贝叶斯学派则把看成是一个随机变量。inference（求后验）和prediction下面结合图来说明贝叶斯线性回归的过程.
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值