I wrote a blog post on cmdstanpy and arviz integration that uses xarray-einstats for posterior predictive sampling:
Once you have an InferenceData object, nothing depends on coming from pymc, stan or even turing in julia so a large part of the post is also useful as an example of using xarray-einstats on pymc results.