Dimension Mismatch for az.plot_hdi

Questions version agnostic
1

Hi there, I am beginning level leaner of pymc. I tried to build the following model to learn Prior and Posterior Predictive Checks. I ran the code to plot the predicted relationship. However, I got the error (TypeError: Dimension mismatch for x: (500, 1) and hdi: (1, 500, 2). Check the dimensions of y andhdi_kwargs to make sure they are compatible)

I checked the shape of posterior_predictive[‘y’]. It shows (4, 1000, 1, 500). It looks different from the example in Posterior Predictive Checks. I can’t figure out what went wrong. Can someone help me?

X = (df1[["temp_actual"]]).values
obs = (df1[['rides']]).values.T

with pm.Model() as model_1:
    
    # Intercepts Priors
    intercept = pm.Normal("intercept", mu=3000, sigma=1500)
    
    # Slope Priors
    slope = pm.Normal("slope", mu=100, sigma=40)

    # Expected value
    mu = intercept + slope.dot(X.T)
    sigma = pm.Exponential("sigma", 0.0008)

    # Data likelihood
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=obs)

    idata_1 = pm.sample_prior_predictive(samples=100)
    idata_1.extend(pm.sample(1000, tune=2000))
    pm.sample_posterior_predictive(idata_1, extend_inferencedata=True)

# Plot the predicted relationship

post = idata_1.posterior
mu_pp = post["intercept"] + post["slope"] * xr.DataArray(df1['temp_actual'].values)

_, ax = plt.subplots()

ax.plot(
    X, mu_pp.mean(("chain", "draw")), label="Mean outcome", color="C1", alpha=0.6
)
ax.scatter(X, idata_1.observed_data["y"])
az.plot_hdi(X, idata_1.posterior_predictive["y"])

ax.set_xlabel("Temp Actual")
ax.set_ylabel("Rides");
2

You might check package versions. I am not able to replicate your issue. Granted I don’t have your data but try this and see if it works.

import pymc as pm
import pymc.sampling.jax as pmjax
import arviz as az
import xarray as xr
import matplotlib.pyplot as plt
import seaborn as sns

# Load data
df = sns.load_dataset("penguins").dropna()

X = df["flipper_length_mm"].values
obs = df["body_mass_g"].values

with pm.Model() as model_1:
    
    # Intercepts Priors
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    
    # Slope Priors
    slope = pm.Normal("slope", mu=0, sigma=10)

    # Expected value
    mu = intercept + slope*X #slope.dot(X.T)
    sigma = pm.Exponential("sigma", .5)

    # Data likelihood
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=obs)

    idata_1 = pm.sample_prior_predictive(samples=100)
    idata_1.extend(pmjax.sample_numpyro_nuts(1000,chains=4))
    pm.sample_posterior_predictive(idata_1, extend_inferencedata=True)

# Plot the predicted relationship
post = idata_1.posterior
mu_pp = post["intercept"] + post["slope"] * xr.DataArray(X)

_, ax = plt.subplots()

ax.plot(
    X, mu_pp.mean(("chain", "draw")), label="Mean outcome", color="C1", alpha=0.6
)
ax.scatter(X, idata_1.observed_data["y"])
az.plot_hdi(X, idata_1.posterior_predictive["y"])
1 Like
3

Hi Justin, thank you so much for your help. It works perfectly. I may have overly complicated the model. After simplifying it by your suggestions, it looks good.

1 Like