使用 pymc3 拟合 lomax 模型答案

【问题标题】：Using pymc3 to fit lomax model使用 pymc3 拟合 lomax 模型
【发布时间】：2018-05-22 22:39:30
【问题描述】：

我有一个非常简单的示例，但似乎不起作用。我的目标是建立一个 Lomax 模型，并且由于 PyMC3 没有 Lomax 分布，我使用了一个事实，即与 Gamma 混合的指数是 Lomax（请参阅here）：

import pymc3 as pm
from scipy.stats import lomax

# Generate artificial data with a shape and scale parameterization
data = lomax.rvs(c=2.5, scale=3, size=1000)

# if t ~ Exponential(lamda) and lamda ~ Gamma(shape, rate), then t ~ Lomax(shape, rate)
with pm.Model() as hierarchical:
    shape = pm.Uniform('shape', 0, 10)
    rate = pm.Uniform('rate', 0 , 10)
    lamda = pm.Gamma('lamda', alpha=shape, beta=rate)
    t = pm.Exponential('t', lam=lamda, observed=data)
    trace = pm.sample(1000, tune=1000)

总结是：

>>> pm.summary(trace)
           mean        sd  mc_error   hpd_2.5  hpd_97.5   n_eff      Rhat
shape  4.259874  2.069418  0.060947  0.560821  8.281654  1121.0  1.001785
rate   6.532874  2.399463  0.068837  2.126299  9.998271  1045.0  1.000764
lamda  0.513459  0.015924  0.000472  0.483754  0.545652  1096.0  0.999662

我预计形状和速率估计值分别接近 2.5 和 3。我尝试了各种形状和速率的非信息性先验，包括pm.HalfFlat() 和pm.Uniform(0, 100)，但两者都导致更差的拟合。有什么想法吗？

【问题讨论】：

标签： bayesian pymc3

【解决方案1】：

想通了：要从指数伽马混合物中导出 lomax，我需要为数据集中的每个示例指定一个 lamda (lamda = pm.Gamma('lamda', alpha=shape, beta=rate, shape=len(data))。这是因为模型假设数据中的每个主题都有自己的lamda_i，其中lamda_i ~ Gamma(shape, rate) 对应每个i。

【讨论】：