第3章 Python概率编程入门

第3章：Python概率编程入门：PyMC库安装、基本模型构建、MCMC采样初体验

好，咱们开始动手了。

前面两章聊了不少贝叶斯的理论，什么先验、后验、似然函数，听着挺唬人。但说实话，我当年学的时候，光看公式差点劝退。直到我真正打开Python，装好PyMC，跑通第一个MCMC采样——那种「原来如此」的感觉，才让我彻底爱上这套工具。

这一章，咱们就干三件事：装好PyMC，搭一个最简单的贝叶斯模型，然后让MCMC帮我们采样。你跟着我走一遍，会发现这东西比想象中简单得多。

3.1 安装PyMC：别踩我踩过的坑

PyMC是目前Python生态里最主流的概率编程库之一。我个人的习惯是，所有新项目都先建一个虚拟环境，免得依赖打架。

⚠️ 我曾经踩过的坑： 直接在base环境里装PyMC，结果跟TensorFlow的旧版本冲突，折腾了一下午。所以，请务必用虚拟环境。

安装命令很简单：

# 推荐用conda创建环境
conda create -n bayesian python=3.10
conda activate bayesian

# 安装PyMC
pip install pymc

# 验证安装
python -c "import pymc as pm; print(pm.__version__)"

如果你用的是M1/M2芯片的Mac，可能会遇到一些底层编译问题。我建议直接用conda安装，它会自动帮你处理底层依赖：

conda install -c conda-forge pymc

装完之后，咱们先画个知识地图，看看这一章要学什么：

3.2 第一个贝叶斯模型：抛硬币问题

咱们从最经典的例子开始——抛硬币。你想知道这枚硬币是不是均匀的。

假设我们抛了10次，得到7次正面、3次反面。用频率学派的方法，点估计是0.7。但贝叶斯学派会说：「别急，我们先给θ一个先验，再看看数据怎么说。」

核心思路： 先验(θ) × 似然(数据|θ) → 后验(θ|数据)

在PyMC里，模型构建分三步走：

定义先验：用 pm.Beta 分布表示我们对θ的初始信念
定义似然：用 pm.Binomial 分布描述数据生成过程
运行采样：用MCMC从后验中抽取样本

看代码：

import pymc as pm
import arviz as az
import numpy as np

# 模拟数据：10次抛掷，7次正面
n_trials = 10
n_heads = 7

# 构建模型
with pm.Model() as coin_model:
    # 第一步：定义先验 - Beta(2,2) 表示我们认为硬币大概率是均匀的
    theta = pm.Beta('theta', alpha=2, beta=2)
    
    # 第二步：定义似然 - 二项分布
    y = pm.Binomial('y', n=n_trials, p=theta, observed=n_heads)
    
    # 第三步：MCMC采样
    trace = pm.sample(draws=2000, tune=1000, random_seed=42)

# 查看结果
az.summary(trace, var_names=['theta'])

跑完之后，你会看到类似这样的输出：

参数	均值	标准差	HDI 3%	HDI 97%	R-hat
theta	0.64	0.12	0.41	0.86	1.001

看到了吗？后验均值是0.64，不是0.7。为什么？因为先验Beta(2,2)把结果往0.5拉了一点。这就是贝叶斯推断的「收缩效应」——数据量少的时候，先验会起到平滑作用。

💡 我的经验： 刚开始用PyMC时，我总担心采样不收敛。后来养成一个习惯：每次跑完都看一眼R-hat值。如果R-hat接近1.0（一般<1.01），说明采样没问题。如果大于1.05，嗯...你得增加采样步数或者检查模型了。

3.3 MCMC采样到底在干嘛？

你可能会问：「MCMC采样到底在干嘛？为什么能帮我们得到后验分布？」

说白了，MCMC就像一个聪明的探索者。它在我们关心的参数空间里随机游走，但又不是完全瞎走——它倾向于停留在后验概率高的区域。走的时间越长，留下的「脚印」就越能反映后验分布的形状。

PyMC默认用的是NUTS采样器（No-U-Turn Sampler），这是HMC（Hamiltonian Monte Carlo）的改进版。我个人觉得，NUTS最大的好处就是不用手动调步长，省心不少。

咱们可以画一下采样轨迹，看看采样器是怎么工作的：

import matplotlib.pyplot as plt

# 绘制采样轨迹
az.plot_trace(trace, var_names=['theta'])
plt.show()

你会看到两张图：左边是后验分布的密度图，右边是采样轨迹。轨迹图里那条「毛毛虫」一样的线，就是MCMC在参数空间里游走的路径。如果它看起来像白噪声，没有明显的趋势或漂移，那就说明采样收敛了。

⚠️ 注意： 如果轨迹图有明显的「爬坡」或「下坡」趋势，说明还没收敛。我遇到过这种情况，当时把采样步数从1000加到5000才稳定下来。别慌，这是新手常遇到的问题。

3.4 换个先验试试看

贝叶斯分析里，先验的选择很关键。咱们试试不同的先验，看看后验会怎么变：

先验分布	含义	后验均值
Beta(1,1)	无信息先验，认为所有θ等可能	0.67
Beta(2,2)	弱信息先验，偏向均匀	0.64
Beta(10,10)	强信息先验，强烈认为硬币均匀	0.59

你看，先验越强，后验越被拉向0.5。这就是贝叶斯统计的「先验主导」现象。在量化投资里，如果你对某个参数有很强的先验信念（比如夏普比率不会超过2），那就用强先验；如果没什么把握，就用弱先验，让数据说话。

  实战建议： 我一般会跑多个先验做敏感性分析。如果后验对先验不敏感，说明数据信息足够强；如果很敏感，那就得小心了——你的结论可能更多来自假设，而不是数据。

3.5 保存和加载采样结果

MCMC采样通常比较耗时，尤其是复杂模型。我习惯把采样结果存下来，方便后续分析：

# 保存
az.to_netcdf(trace, 'coin_trace.nc')

# 加载
trace_loaded = az.from_netcdf('coin_trace.nc')

这样下次打开Jupyter Notebook就不用重新采样了。嗯，这个小习惯帮我省了不少时间。

3.6 本章小结

咱们这一章做了三件事：

装好了PyMC，避开了环境冲突的坑
用抛硬币的例子，搭了第一个贝叶斯模型
跑了MCMC采样，看到了后验分布长什么样

说白了，概率编程就是把你的统计假设写成代码，然后让计算机帮你算后验。你不需要手动推导复杂的积分公式——PyMC和MCMC帮你搞定了。

下一章，咱们会把这个框架用到更实际的量化场景里，比如估计股票的收益率分布。到时候你会发现，贝叶斯方法在处理小样本和不确定性方面，确实有独到之处。

公众号:蓝海资料掘金营，微信deep3321

📚 课程目录