PyMC3：概率统计建模的得力助手！
大家好！今天我要介绍一个特别强大的统计建模工具 - PyMC3。它是一个专门用于贝叶斯统计分析的Python库，能让我们轻松构建概率模型。如果你觉得统计很难，别担心！PyMC3就像是一位贴心的向导，用简单的代码就能帮我们完成复杂的统计分析。
初识PyMC3
首先，让我们安装并导入必要的包：
# 安装PyMC3
pip install pymc3
# 导入必要的包
import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
简单的概率模型
让我们从一个简单的例子开始 -估计硬币的正面概率：
# 生成示例数据：假设我们抛了100次硬币，其中60次正面
n_trials = 100 # 试验次数
n_heads = 60 # 正面次数
# 构建模型
with pm.Model() as coin_model:
# 设置先验分布（我们对硬币正面概率的初始猜测）
theta = pm.Beta('theta', alpha=1, beta=1)
# 设置似然函数（观察到的数据如何产生）
y = pm.Binomial('y', n=n_trials, p=theta, observed=n_heads)
# 进行推断
trace = pm.sample(2000)
# 查看结果
pm.plot_posterior(trace['theta'])
plt.show()
小贴士： Beta分布是描述概率的理想选择，因为它的取值范围正好在0到1之间，完美符合概率的定义！
线性回归模型
接下来，让我们用PyMC3构建一个贝叶斯线性回归模型：
# 生成一些示例数据
np.random.seed(42)
X = np.random.randn(100, 1)
true_slope = 2.0
true_intercept = 1.0
y = true_intercept + true_slope * X + np.random.randn(100, 1) * 0.5
# 构建贝叶斯线性回归模型
with pm.Model() as linear_model:
# 设置先验分布
intercept = pm.Normal('intercept', mu=0, sd=10)
slope = pm.Normal('slope', mu=0, sd=10)
sigma = pm.HalfNormal('sigma', sd=1)
# 设置似然函数
mu = intercept + slope * X
y_obs = pm.Normal('y_obs', mu=mu, sd=sigma, observed=y)
# 进行推断
trace = pm.sample(2000)
# 绘制结果
pm.plot_trace(trace)
plt.show()
实际应用：时间序列分析
PyMC3在时间序列分析中也非常有用。让我们来预测一些数据的趋势：
# 创建示例时间序列数据
days = np.arange(0, 100)
trend = 0.1 * days
seasonal = 2 * np.sin(2 * np.pi * days / 7)
y = trend + seasonal + np.random.normal(0, 0.2, 100)
# 构建时间序列模型
with pm.Model() as time_series_model:
# 趋势项
slope = pm.Normal('slope', mu=0, sd=1)
intercept = pm.Normal('intercept', mu=0, sd=10)
# 季节性项
amplitude = pm.HalfNormal('amplitude', sd=2)
phase = pm.Uniform('phase', lower=0, upper=2*np.pi)
# 组合模型
trend = intercept + slope * days
seasonal = amplitude * np.sin(2*np.pi * days/7 + phase)
sigma = pm.HalfNormal('sigma', sd=1)
# 观测值
y_obs = pm.Normal('y_obs', mu=trend + seasonal, sd=sigma, observed=y)
# 采样
trace = pm.sample(1000)
实践小练习
试试看能不能用PyMC3解决这个问题：假设你在测量一批产品的重量，得到了以下数据：
weights = np.array([10.2, 10.5, 9.8, 10.1, 10.3, 10.0, 10.4, 10.2, 9.9, 10.3])
# 问题：
# 1. 估计这批产品的真实平均重量
# 2. 估计测量的误差范围
# 3. 计算重量在10.0到10.5之间的概率
# 提示：
# 1. 使用Normal分布描述测量值
# 2. 使用HalfNormal分布描述标准差
# 3. 可以用trace来计算后验概率
实用小贴士：
在使用PyMC3时，先画出概率图模型可以帮助我们理清思路 多次运行采样来确保结果的稳定性 使用 pm.summary()可以快速查看posterior分布的统计特征
今天的Python学习之旅就到这里啦！记得动手实践，亲自体验贝叶斯分析的魅力。统计分析其实可以很有趣，尤其是有PyMC3这样强大的工具相助。祝大家学习愉快，Python学习节节高！
