I have tried to replicate a number of examples using pymc3 and compared the results. Below is the example for estimating 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)
The summary of the example shows hpd intervals (posterior's HDI) - the image is a print screen from the example website:
az.summary(trace)
When I try the same command, however, I get HDI:
I am not sure the parameters represent the same thing; the version of pymc3 used in the example is 3.10.0 whereas the version I used to run the example is 3.11.5.
Does anyone know if the naming convention has been changed or is there something else to be adapted in order to get the actual HPD values in newer version?
They are the same. hpd
was renamed to hdi
to be clear in that it can be used to compute highest density intervals of any quantity, not only of the posterior. See this GitHub PR if interested in the renaming.