这是我第一次尝试让来自常客阵营的人进行贝叶斯数据分析。我阅读了 A. Gelman 的 Bayesian Data Analysis 的许多教程和几章。
作为我选择的第一个或多或少独立的数据分析示例是火车等待时间。我问自己:等待时间的分布是怎样的?
该数据集在博客上提供,在 PyMC 之外进行了略微不同的分析。
我的目标是根据这 19 个数据条目估计预期的火车等待时间。
我建立的模型如下:
其中是数据均值,是数据标准差乘以 1000。
我使用泊松分布将预期等待时间建模为该分布的速率参数使用 Gamma 分布建模,因为它是泊松分布的共轭分布。超先验和分别用正态分布和半正态分布建模。标准偏差尽可能宽泛,尽可能不明确。
我有一堆问题
- 这个模型对于任务是否合理(几种可能的建模方式?)?
- 我犯了任何初学者的错误吗?
- 模型可以简化吗(我倾向于把简单的事情复杂化)?
- 如何验证速率参数 ( ) 的后验是否实际拟合数据?
- 如何从拟合的泊松分布中抽取一些样本来查看样本?
我也可以发布源代码。在模型拟合阶段,我使用 NUTS和然后在第二步中,我为速率参数做 Metropolis 。最后,我使用内置工具绘制了轨迹。
我将非常感谢任何能让我掌握更多概率编程的评论和评论。可能还有更多值得尝试的经典例子吗?
这是我使用 PyMC3 在 Python 中编写的代码。数据文件可以在这里找到。
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import pymc3
from scipy import optimize
from pylab import figure, axes, title, show
from pymc3.distributions import Normal, HalfNormal, Poisson, Gamma, Exponential
from pymc3 import find_MAP
from pymc3 import Metropolis, NUTS, sample
from pymc3 import summary, traceplot
df = pd.read_csv( 'train_wait.csv' )
diff_mean = np.mean( df["diff"] )
diff_std = 1000*np.std( df["diff"] )
model = pymc3.Model()
with model:
# unknown model parameters
mu = Normal('mu',mu=diff_mean,sd=diff_std)
sd = HalfNormal('sd',sd=diff_std)
# unknown model parameter of interest
rate = Gamma( 'rate', mu=mu, sd=sd )
# observed
diff = Poisson( 'diff', rate, observed=df["diff"] )
with model:
step1 = NUTS([mu,sd])
step2 = Metropolis([rate])
trace = sample( 5000, step=[step1,step2] )
plt.figure()
traceplot(trace)
plt.savefig("rate.pdf")
plt.show()
plt.close()