Visualizing highest posterior density for multiple conditions using arviz plot_hpd

Hi,

There are a couple of small modifications to your code that should fix the problem, however, I’d recommend taking advantage of ArviZ and xarray as it is shown in this notebook. If you have any trouble with the notebook please ask here below, I hope the model will be similar enough and the explanations clear enough.

for i in range(groups):
    
    alpha = trace_rsp_rt['α'][:, i]
    beta = trace_rsp_rt['β'][:, i]
    mu = alpha + beta * rt  
    # there may be broadcasting issues requiring to use rt[None, :]
    # xarray would handle broadcasting automatically ass seen in the notebook
    
    ax_rsp_rt[i].plot(rt, mu.mean(), c='k', label= f'rsp = {alpha:.2f} + {beta:.2f} * rt')
    az.plot_hpd(rt, mu, credible_interval=0.98, color='k', ax=ax_rsp_rt[i])
    ax_rsp_rt[i].legend()
    # combining pyplot and object based commands can yield unexpected results
    ax.set_xlim(1.2, 1.8)  
    ax.set_ylim(0.6, 1) 

On a somewhat tangential topic, note that μ in trace_rsp_rt is the distribution of the mean, not the posterior predictive distribution which is what should be used for predicting and for posterior predictive checks (comparing the predictions of the model with observed data).

EDIT: The notebook is not yet available in https://docs.pymc.io/ but will be updated soon.

3 Likes