基于贝叶斯估计的横波速度预测（回归模型）

一、 查看数据

import pandas as pd
import pymc as pm
import matplotlib.pyplot as plt
import numpy as np
# 读取数据
plt.rcParams['font.sans-serif'] = ['SimHei']  # 设置中文字体为黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示为方块的
file_path = 'data.csv'
data = pd.read_csv(file_path)
# 转换数据类型
data['VP'] = pd.to_numeric(data['VP'], errors='coerce')
data['VS'] = pd.to_numeric(data['VS'], errors='coerce')

import matplotlib.pyplot as plt

# 绘制直方图
data.hist(bins=50, figsize=(10, 2))
plt.show()
max_values = data.max()
min_values = data.min()
median_values = data.median()
mean_values = data.mean()

# Print the results
print("最大值:")
print(max_values)
print("\n最小值:")
print(min_values)
print("\n中位数:")
print(median_values)
print("\n均值:")
print(mean_values)

二、定义贝叶斯回归模型

# 定义贝叶斯线性回归模型
with pm.Model() as linear_regression:
    # 系数alpha的先验分布
    alpha = pm.Normal('alpha', mu=0, sigma=10000)
    # 系数beta的先验分布
    beta = pm.Normal('beta', mu=0, sigma=100)
    # 模型误差的先验分布
    sigma = pm.Normal('sigma',mu=0, sigma=1000)  
    # 预测值的计算
    mu = alpha + beta * data['VP']
    # 似然函数
    VS = pm.Normal('VS', mu=mu, sigma=sigma, observed=data['VS'])
    # 进行MCMC采样
    trace = pm.sample(30000, tune=2000, chains=2)
# 查看采样结果的摘要
pm.summary(trace)

三、绘制后验分布和回归结果

# 绘制后验分布的直方图
import seaborn as sns
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
plt.figure(figsize=(10, 6))
sns.distplot(trace.posterior['alpha'].values, kde=True, bins=50)
plt.title('后验分布的直方图')
plt.xlabel('alpha')
plt.ylabel('频率')
plt.show()
plt.figure(figsize=(10, 6))
sns.distplot(trace.posterior['beta'].values, kde=True, bins=50)
plt.title('后验分布的直方图')
plt.xlabel('beta')
plt.ylabel('频率')
plt.show()
plt.figure(figsize=(10, 6))
sns.distplot(trace.posterior['sigma'].values, kde=True, bins=50)
plt.title('后验分布的直方图')
plt.xlabel('sigma')
plt.ylabel('频率')
plt.show()

# 绘制回归结果
plt.scatter(data['VP'], data['VS'], label='实际数据')
alpha_samples = trace.posterior['alpha'].values
beta_samples = trace.posterior['beta'].values
x = np.linspace(data['VP'].min(), data['VP'].max(), 100)
y_mean = alpha_samples.mean() + beta_samples.mean() * x
y_upper = alpha_samples.mean() + beta_samples.mean() * x + 1.96 * trace.posterior['sigma'].values.mean()
y_lower = alpha_samples.mean() + beta_samples.mean() * x - 1.96 * trace.posterior['sigma'].values.mean()
plt.fill_between(x, y_lower, y_upper, color='gray', alpha=0.3, label='95%置信区间')
plt.plot(x, y_mean, color='red', label='回归线')

# 添加回归方程
equation = f"VS = {alpha_samples.mean():.2f} + {beta_samples.mean():.2f} * VP"
plt.text(0.5, 1.05, equation, transform=plt.gca().transAxes, fontsize=12, verticalalignment='center', horizontalalignment='center')
plt.legend()
plt.xlabel('VP(m/s)')
plt.ylabel('VS(m/s)')
plt.show()