火车等待时间的贝叶斯建模:模型定义

机器算法验证 贝叶斯 pymc
2022-03-04 10:05:32

这是我第一次尝试让来自常客阵营的人进行贝叶斯数据分析。我阅读了 A. Gelman 的 Bayesian Data Analysis 的许多教程和几章。

作为我选择的第一个或多或少独立的数据分析示例是火车等待时间。我问自己:等待时间的分布是怎样的?

该数据集在博客上提供,在 PyMC 之外进行了略微不同的分析。

我的目标是根据这 19 个数据条目估计预期的火车等待时间。

我建立的模型如下:

μN(μ^,σ^)

σ|N(0,σ^)|

λΓ(μ,σ)

ρPoisson(λ)

其中是数据均值,是数据标准差乘以 1000。μ^σ^

我使用泊松分布将预期等待时间建模为该分布的速率参数使用 Gamma 分布建模,因为它是泊松分布的共轭分布。超先验分别用正态分布和半正态分布建模。标准偏差尽可能宽泛,尽可能不明确。ρμσσ

我有一堆问题

  • 这个模型对于任务是否合理(几种可能的建模方式?)?
  • 我犯了任何初学者的错误吗?
  • 模型可以简化吗(我倾向于把简单的事情复杂化)?
  • 如何验证速率参数 ( ) 的后验是否实际拟合数据?ρ
  • 如何从拟合的泊松分布中抽取一些样本来查看样本?

5000 Metropolis 步后的后验是这样的: 跟踪图

我也可以发布源代码。在模型拟合阶段,我使用 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()
1个回答

我会先告诉你我会做什么,然后我会回答你的具体问题。

我会做什么(至少在最初)

这是我从您的帖子中收集到的信息,您有 19 次观察的训练等待时间,并且您有兴趣推断出预期的等待时间。

我将定义为的等待时间我认为这些等待时间没有理由应该是整数,所以我假设它们是正连续量,即我假设实际上观察到了所有的等待时间。Wii=1,,19iWiR+

有几种可能的模型假设可以使用,并且通过 19 次观察,可能很难确定哪个模型更合理。一些例子是对数正态、伽马、指数、威布尔。

作为第一个模型,我建议建模,然后假设 有了这个选择,您就可以获得大量现有的正规理论,例如共轭先验。共轭先验是正态逆伽马分布,即 其中是逆伽马分布。或者,您可以使用默认的先验在这种情况下,后验也是正态逆伽马分布。Yi=log(Wi)

YiindN(μ,σ2).
μ|σ2N(m,σ2C)σ2IG(a,b)
IGp(μ,σ2)1/σ2

由于的联合样本来回答有关预期等待时间的问题,这是一个正态逆分布-gamma 分布,然后计算每个样本的这从后部对预期的等待时间进行采样。E[Wi]=eμ+σ/2μσ2eμ+σ/2

回答您的问题

  • 这个模型对于任务是否合理(几种可能的建模方式?)?

泊松似乎不适用于可能是非整数值的数据。您只有一个,因此您无法学习分配给的伽马分布参数。另一种说法是,您已经建立了层次模型,但数据中没有层次结构。λλ

  • 我犯了任何初学者的错误吗?

见之前的评论。

此外,如果您的数学和您的代码一致,这将非常有帮助,例如,您的 MCMC 结果中的您的代码中的 sd 和 rate 是什么?λ

您的先验不应依赖于数据。

  • 模型可以简化吗(我倾向于把简单的事情复杂化)?

是的,它应该。请参阅我的建模方法。

  • 如何验证速率参数 ( ) 的后验是否实际拟合数据?ρ

应该是你的数据吗?你的意思是吗?要检查的一件事是确保样本平均等待时间相对于平均等待时间的后验分布有意义。除非你有一个奇怪的先验,否则样本平均值应该接近后验分布的峰值。ρλ

  • 如何从拟合的泊松分布中抽取一些样本来查看样本?

我相信你想要一个后验预测分布。对于 MCMC 中的每次迭代,您插入该迭代的参数值并进行采样。