使用pymc3拟合lomax模型
Using pymc3 to fit lomax model
我有一个非常简单的示例,但似乎不起作用。我的目标是构建一个 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)
,但两者都导致更差的拟合。有什么想法吗?
想出来了:要从指数-伽马混合中导出 lomax,我需要为数据集中的每个示例指定一个 lamda
(lamda = pm.Gamma('lamda', alpha=shape, beta=rate, shape=len(data)
)。这是因为该模型假设数据中的每个主题都有自己的 lamda_i
,其中 lamda_i ~ Gamma(shape, rate)
对应每个 i
。
我有一个非常简单的示例,但似乎不起作用。我的目标是构建一个 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)
,但两者都导致更差的拟合。有什么想法吗?
想出来了:要从指数-伽马混合中导出 lomax,我需要为数据集中的每个示例指定一个 lamda
(lamda = pm.Gamma('lamda', alpha=shape, beta=rate, shape=len(data)
)。这是因为该模型假设数据中的每个主题都有自己的 lamda_i
,其中 lamda_i ~ Gamma(shape, rate)
对应每个 i
。