我尝试使用pymc3复制一些示例,并对结果进行比较。下面是估算HPD的示例:
import pymc3
import arviz as az
import numpy as np
import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
import graphviz
import pymc3 as pm
from pymc3 import Model, Normal, HalfNormal
from pymc3 import find_MAP
basic_model = Model()
with basic_model:
# Priors for unknown model parameters
alpha = Normal('alpha', mu=0, sd=5)
beta = Normal('beta', mu=0, sd=5, shape=2)
sigma = HalfNormal('sigma', sd=4)
# Expected value of outcome
mu = alpha + beta[0]*X1 + beta[1]*X2
# Deterministic variable, to have PyMC3 store mu as a value in the trace use
# mu = pm.Deterministic('mu', alpha + beta[0]*X1 + beta[1]*X2)
# Likelihood (sampling distribution) of observations
Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
pm.model_to_graphviz(basic_model)
with basic_model:
# obtain starting values via MAP
start = find_MAP(fmin=optimize.fmin_powell)
# instantiate sampler - not really a good practice
step = NUTS(scaling=start)
# draw 2000 posterior samples
trace = sample(2000, step, start=start)该示例的摘要显示了hpd间隔(后验HDI) -图像是来自示例网站的打印屏幕:
az.summary(trace)

然而,当我尝试使用相同的命令时,我得到了HDI:

我不确定参数是否代表相同的内容;示例中使用的pymc3版本为3.10.0,而用于运行示例的版本为3.11.5。
是否有人知道命名约定是否已经更改,或者是否有其他需要修改的东西,以便在更新版本中获得实际的HPD值?
发布于 2022-07-20 16:38:59
他们是一样的。为了清楚地将hpd重命名为hdi,它可以用来计算任意数量的最高密度间隔,而不仅仅是后验。如果对重命名感兴趣,请参见这个GitHub PR。
https://stackoverflow.com/questions/72914530
复制相似问题