# Sample posterior predictive with a vengeance
**URL:** https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926
**Category:** Questions
**Created:** [October 7, 2020, 2:24pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926 "2020-10-07T14:24:41Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@bridgeland](https://discourse.pymc.io/u/bridgeland)
#### Post date: [October 7, 2020, 2:24pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/1 "2020-10-07T14:24:41Z")
My model gives some nice results. Now I would like to investigate what would happen if the value of some RVs were to be changed, e.g. by an outside intervention. In particular I would like a function that is much like `sample_posterior_predictive()`, but has an additional parameter, something like this:
```auto
sample_posterior_intervention(
trace, var_names=['results'],
intervention={'mu': lambda x: x+1, 'sigma': lambda x: x * 0.5}
)
```
meaning run the posterior intervention, but instead of using the value of `mu` in the trace, use `mu + 1`, and instead of using the value of `sigma`, use `sigma * 0.5`.
I appreciate that `pymc3.sampling` has nothing like my imagined `sample_posterior_intervention()` already. How could I implement this behavior using the existing `pymc3.sampling` functions?
Or is this a nonsensical thing to do?
---
### Author:  [@RavinKumar](https://discourse.pymc.io/u/RavinKumar)
#### Post date: [October 7, 2020, 11:32pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/2 "2020-10-07T23:32:58Z")
What you can do is use shared variables in Theano to run inference with certain values, and then adjust them after inference, and then run posterior predictive simulations.
There is an example in the docs linked below. Let us know if you have questions!
> **[Advanced usage of Theano in PyMC3 — PyMC3 3.9.3 documentation](https://docs.pymc.io/Advanced_usage_of_Theano_in_PyMC3.html)**
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [October 8, 2020, 1:31am UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/3 "2020-10-08T01:31:30Z")
I understand that you want to modify posterior samples for some of the variables in the posterior. In that case you can edit the `posterior` xarray Dataset of your result (given it’s stored as inferencedata) and then pass it to `sample_posterior_predictive`
You can also choose to also combine this with the shared objects and `pm.Data` to modify predictors instead of posterior samples, but doing only one of the 2 is also possible.
---
### Author:  [@simeoncarstens](https://discourse.pymc.io/u/simeoncarstens)
#### Post date: [October 8, 2020, 2:39pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/4 "2020-10-08T14:39:29Z")
> Or is this a nonsensical thing to do?
You most certainly are, but be aware that the resulting “posterior predictive” samples don’t have any clear statistical meaning - after all, they are _not_ drawn from the posterior predictive distribution.
---
### Author:  [@bridgeland](https://discourse.pymc.io/u/bridgeland)
#### Post date: [October 8, 2020, 2:54pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/5 "2020-10-08T14:54:01Z")
Can a variable be defined as both a shared variable and an RV that is fit to the data in the original sample? I don’t see any examples of this in the documentation.
Maybe a shared variable can be defined for the intervention, with the original RV as a term in a newly created deterministic. So instead of this
```auto
with pm.Model() as model:
mu = pm.Normal('mu', 0, 10)
```
an intervention-ready `mu` is defined like this:
```auto
with pm.Model() as model:
mu = pm.Deterministic('mu', mu_original + mu_intervention)
mu_original = pm.Normal('mu_original', 0, 10)
mu_intervention = theano.shared(0.0)
```
Then after the original sampling, the intervention is run and the posterior is re-sampled, like this:
```auto
with model:
mu_intervention.set_value(1.0)
intervention_trace = sample_posterior_predictive(
trace, var_names=['results'])
```
---
### Author:  [@bridgeland](https://discourse.pymc.io/u/bridgeland)
#### Post date: [October 8, 2020, 3:29pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/6 "2020-10-08T15:29:41Z")
Your understanding is correct @OriolAbril.
I think you are suggesting I change the MultiTrace object that is provided by `sample()`. (If not, please correct my misunderstanding.)
How can the MultiTrace object be changed? `trace['mu']` provides an numpy array of the values of `mu`, one element in the array for each draw. But changing those values seems to have no effect on the original `trace`. It appears that `trace['mu']` creates a new copy of the array.
---
### Author:  [@OriolAbril](https://discourse.pymc.io/u/OriolAbril)
#### Post date: [October 8, 2020, 3:52pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/7 "2020-10-08T15:52:45Z")
You have to use `arviz.InferenceData` instead of `pm.MultiTrace`. I’d also recommend trying to switch from MultiTrace to InferenceData as it should be more robust and flexible. PyMC3 is progressively working towards integration with InferenceData, you can now use `return_inferencedata=True` in `pm.sample` (which will become the default in the long run) and [`pm.sample_posterior_predictive`](https://docs.pymc.io/api/inference.html#pymc3.sampling.sample_posterior_predictive) already accepts both xarray dataset and InferenceData in addition to multitrace.
Here are some references that could be useful:
- Introduction to InferenceData in ArviZ docs: [https://arviz-devs.github.io/arviz/notebooks/XarrayforArviZ.html](https://arviz-devs.github.io/arviz/notebooks/XarrayforArviZ.html)
- PyMC3 example notebook showcasing some of the features of using InferenceData and xarray: [https://docs.pymc.io/notebooks/multilevel\_modeling.html](https://docs.pymc.io/notebooks/multilevel_modeling.html)
- A blog post (written by me) about PyMC3\<\>ArviZ integration and the use of labeled coords and dims with your models:
> **[PyMC3 with labeled coords and dims](https://oriolabril.github.io/oriol_unraveled/python/arviz/pymc3/xarray/2020/09/22/pymc3-arviz.html)**
>
> Go crazy with your virtual label-maker!
---
### Author:  [@bridgeland](https://discourse.pymc.io/u/bridgeland)
#### Post date: [October 9, 2020, 8:16pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/8 "2020-10-09T20:16:51Z")
For someone reading this in the future, here’s the resulting code. (I have not actually run this. My actual model is more complex than this simple example. So there may be a bug in my simplified example here. Caveat emptor.)
```auto
import xarray as xr
def change_mu_and_sigma_in_posterior(posterior_ds):
"""Create new posterior dataset with mu and sigma altered"""
new_post_ds = posterior_ds.copy(deep=False)
new_post_ds['mu'] = xr.apply_ufunc(
lambda mu_da: mu_da + 1,
new_post_ds['mu'])
new_post_ds['sigma'] = xr.apply_ufunc(
lambda sigma_da: sigma_da * 0.5,
new_post_ds['sigma'])
return new_post_ds
with model:
trace_with_intervention = trace.map(
change_mu_and_sigma_in_posterior,
groups='posterior')
posterior_post_intervention_results = pm.sample_posterior_predictive(
trace_with_intervention)
```
---
### Author:  [@bridgeland](https://discourse.pymc.io/u/bridgeland)
#### Post date: [October 13, 2020, 1:24pm UTC](https://discourse.pymc.io/t/sample-posterior-predictive-with-a-vengeance/5926/9 "2020-10-13T13:24:58Z")
There is a possible third approach, in addition to the two suggested by @RavinKumar and @OriolAbril. In a [post about applying an existing model to new data](https://discourse.pymc.io/t/how-do-we-predict-on-new-unseen-groups-in-a-hierarchical-model-in-pymc3/2571/2), @lucianopaz suggested creating a new model with a model factory. I think his model factory approach could also work for solving this problem: create a new model using the existing data and fit to newly altered posterior samples of the old model.
Note: I haven’t tried this approach.