# Prior predictive sampling with transformed RV

**URL:** <https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420>\
**Category:** Questions\
**Created:** [January 27, 2020, 7:15pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420 "2020-01-27T19:15:01Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [January 27, 2020, 7:15pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/1 "2020-01-27T19:15:01Z")

</div>

I was playing with a Dirichlet distribution with an uninformed prior, as follows:

```auto
with pm.Model() as histogram_model:
    weights = pm.HalfNormal('weights', shape=10, testval=np.array([0.1] * 10), transform=pm.transforms.sum_to_1)
    pm.Dirichlet('obs', a=weights, shape=10)

```

When I sample from the prior predictive, it seems like the transform is not applied:

```auto
with histogram_model:
    prior_trace = pm.sample_prior_predictive(var_names=['weights'])
prior_trace['weights'][0]

```

returns, in one example:

```auto
array([0.32092762, 0.80871371, 0.14739512, 0.4890916 , 0.31236516,
       0.88778043, 1.18220203, 2.7092601 , 1.45828804, 0.76331364])

```

Am I doing something wrong? Or is this to be expected – that the transform won’t be applied?

And if the transform is not going to be applied, is there some way to evaluate a model like this to debug it?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [February 2, 2020, 10:51am UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/2 "2020-02-02T10:51:30Z")

</div>

> [@rpgoldman](#):
>
> Am I doing something wrong? Or is this to be expected – that the transform won’t be applied?

You are correct, the transform only applied on the latent “free” variable, but not on the random method, as the random method is assume to produce random sample already in the right domain.

> [@rpgoldman](#):
>
> And if the transform is not going to be applied, is there some way to evaluate a model like this to debug it?

I think the only way is to do `pm.sample()` directly.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 2, 2020, 6:27pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/3 "2020-02-02T18:27:57Z")

</div>

Should we treat this as a bug in `sample_prior_predictive` and `sample_posterior_predictive`? It really seems wrong to me.

If you agree, I will make an issue for this.

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [February 3, 2020, 4:33pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/4 "2020-02-03T16:33:44Z")

</div>

No, I wouldn’t treat this as a bug from sample prior predictive’s perspective. I think this is a consequence of the sometimes contradictory behavior of `transform`. I think that we should document somewhere the fact that passing transforms like this will only change `sample` but won’t affect prior predictives.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 3, 2020, 4:46pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/5 "2020-02-03T16:46:25Z")

</div>

@lucianopaz But doesn’t this make the prior predictive wrong as a tool for debugging a model? It is going to give results that are outside the support of the distribution, correct?

Even worse, it seems to me, is what will happen if we use posterior predictive sampling for out-of-sample predictions (machine learning applications). In this case we are going to return predictions that are simply wrong, aren’t we?

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [February 3, 2020, 5:14pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/6 "2020-02-03T17:14:43Z")

</div>

It will return wrong samples in this case. But it will return the correct samples for many other cases, like the halfnormal, halfcauchy or even the uniform distribution. In my opinion, transform is at fault here.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 3, 2020, 6:15pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/7 "2020-02-03T18:15:13Z")

</div>

I see your point about the transforms being at fault, but I still think that not applying the transform works against the idea of Bayesian modeling and against the principle of understandable software.

1. Bayesian modeling is based on using generative models and inverting them to reason from evidence. Having generative models that don’t generate properly seems wrong to me.

2. Programmers have every reason to expect that a transformed random variable will yield transformed samples under all circumstances. Having them generate transformed samples under only _some_ circumstances is akin to having the multiplication operator sometimes doing addition instead.

Changing it might be difficult, but I still think that what PyMC3 does here is wrong, and should be fixed. I can’t think of a principled argument – as opposed to an expedient argument – for the current behavior.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [February 3, 2020, 7:13pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/8 "2020-02-03T19:13:51Z")

</div>

I dont think we should change this, FWIW what you are doing in

```python
weights = pm.HalfNormal('weights', shape=10, testval=np.array([0.1] * 10), transform=pm.transforms.sum_to_1)

```

is incorrect, as `HalfNormal` does not produce a vector that sum to 1. What transform does is transferring a (vector) of value in \text{R} into something that are in the domain that defined by the distribution.

In fact, the transformation in PyMC3 is defined in this direction only. It is different than package like TFP, where transformation is defined in both direction and you can apply both forward and inverse transformation.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 3, 2020, 7:38pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/9 "2020-02-03T19:38:11Z")

</div>

> What transform does is transferring a (vector) of value in R into something that are in the domain that defined by the distribution.

> In fact, the transformation in PyMC3 is defined in this direction only. It is different than package like TFP, where transformation is defined in both direction and you can apply both forward and inverse transformation.

I’m sorry: I don’t follow this argument. I’m asking for the transforms to be applied in `sample_posterior_predictive` and `sample_prior_predictive`. This is only applying the transform forward, as you say, so why would it be wrong to do this?

It seems like the argument _for_ doing this is clear: a programmer that builds a model with transformed variables, and then asks PyMC3 to generate samples from that model, will expect to see transformed variables (e.g., no truncated normal samples below the lower bound).

With all due respect, I can’t see how any argument could be made _for_ the current behavior.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [February 3, 2020, 8:19pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/10 "2020-02-03T20:19:20Z")

</div>

What I am trying to say is that, the current behavior never intended to modify the `random` method of a distribution, as the random method is always assuming to produce random sample that satisfy the constrain. Thus user should not expect transform to be apply on `sample_posterior_predictive` and `sample_prior_predictive`.  
It is however a valid feature request to have transform apply on a distribution so that it does the correct forward and inverse on both samples and log\_prob computation, for that we could introduce a new transform distribution class like the one in TFP.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 3, 2020, 9:29pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/11 "2020-02-03T21:29:01Z")

</div>

I see. It seems like the prior and posterior predictive sampling are not “first class citizens” (a term that programming languages people use) in the PyMC3 language.

One thing remains unclear to me: if we have “internal” transformed variables (e.g., hyper-parameters in a hierarchical model), will predictive samples of the observed variables be only untransformed, or will they also be wrong?

E.g., if one has a truncated normal RV as a Free RV, and it is used to give the mean of a child variable, will that child variable’s value be wrong, also? Or does the theano network ensure that the value that goes to the child variable will be transformed?

Like this:  
I = [1, 2, 3]  
μ ~ TruncatedNormal( μ=I, σ=1)  
o ~ TruncatedNormal(μ=μ, σ=1)

[sorry for such a stupid example]  
Now even if I apply the transform to “o” myself, the answer would still be wrong if μ is not transformed.

I don’t know if the value that Theano transmits from μ to o is going to be the transformed value, or the raw value (in this case presumably a sample from an _un-truncated_ Normal). Which is it?

If theano transforms variables like μ, then this can be worked around relatively easily. But if it does not, then the results of sampling could be not just a little wrong, but arbitrarily wrong.

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [February 3, 2020, 9:49pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/12 "2020-02-03T21:49:12Z")

</div>

It’s kind of hard to answer just with writing, I think that a schematic will help explain what’s going on. In your example with the `HalfNormal`, the transform will get in the way and make a separate `FreeRV` called `__log_{your_halfnormal_name}`, and the `HalfNormal` itself will become a `TransformedRV`. The dependency graph from theano will look something like this:

 ![path843](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/b/bec550ce3220d39f13d85990c2b2436f5cf147ea.png)

and in general it would look like this:

 ![path1169](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/6/6c8f56a18e616f4eeac9cc30dca9de1b86f48a9b.png)

In your example with the `sum_to_1` transform you can ask theano to print the `weights` tensor computational graph like this:

```python
>>> theano.printing.debugprint(weights)
ViewOp [id A] 'weights'   
 |Join [id B] ''   
   |TensorConstant{0} [id C]
   |Subtensor{::} [id D] ''   
   | |weights_sumto1__ [id E]
   |Elemwise{sub,no_inplace} [id F] ''   
     |InplaceDimShuffle{x} [id G] ''   
     | |TensorConstant{1} [id H]
     |InplaceDimShuffle{x} [id I] ''   
       |Sum{axis=[0], acc_dtype=float64} [id J] ''   
         |Subtensor{::} [id K] ''   
           |weights_sumto1__ [id E]

```

So, as you can see, basically, pymc3 always puts an rv higher up the hierarchy. When you call `sample`, the transformed (`weights_sumto1__`) will be used as the `FreeRV`, and the `weights` value will be computed as a `Deterministic`. So in the end, the sampled trace will have the samples of both `weights_sumto1__` and `weights`. In sample posterior predictive, nothing very special happens, the values for each of the nodes are grabbed straight out of the trace.

The complicated part happens in sample prior predictive. What goes on is that we first try to draw a value from `weights_sumto1__`, that will raise a `MissingInputError` because we didn’t provide the computational graph’s input value. After that, we simply call the `weights` random method and return its value as it came out, because we assume that the samples from `weights` already lie in the proper sample space.

To get what you want (weights that are constrained to sum to one), you have to be a bit more explicit:

```python
import pymc3 as pm
from theano import tensor as tt

with pm.Model() as histogram_model:
    weights_raw = pm.HalfNormal('weights', shape=10)
    weights = weights_raw / tt.sum(weights_raw, axis=-1, keepdims=True)
    pm.Dirichlet('obs', a=weights, shape=10)

```

If you’re interested in storing the values of `weights` when you `sample` from your model (or draw prior predictive samples), you just have to wrap it with a `Deterministic`.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [February 3, 2020, 9:55pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/13 "2020-02-03T21:55:55Z")

</div>

> [@rpgoldman](#):
>
> I see. It seems like the prior and posterior predictive sampling are not “first class citizens” (a term that programming languages people use) in the PyMC3 language.

That is very correct, also main reason why we have so much headache of shape etc issue right now.

> [@rpgoldman](#):
>
> Or does the theano network ensure that the value that goes to the child variable will be transformed?

It will be correct, consider a subgraph in a Bayesian Graph (i.e., Markov blanket):  
… → A ~ distribution\_1 → B ~ distribution\_3(A, C) → …  
C ~ distribution\_2 -------------/

for log\_prob computation, you replace A in the graph so that:  
… + distribution\_1.log\_prob(input\_A) + distribution\_3(input\_A, input\_C).log\_prob(…) + …

As you can see, A is replace by user/programmatic input in both distribution and value, thus as long as it is transform to the proper constrained space that defined by distribution\_1 and distribution\_2, it is valid. Note how it is different than in the forward computation that distribution\_1 does not depending on the value of A, whereas in the inverse (log\_prob) computation it does depending on it.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 3, 2020, 9:58pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/14 "2020-02-03T21:58:54Z")

</div>

Thank you, @junpenglao! That is a relief to know!

So the work-around here is to take the predictive trace and for any transformed RV, apply the transform. Not ideal, but definitely feasible.

It would be nice to see that predictive sampling becomes a first class citizen in PyMC4.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 3, 2020, 10:08pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/15 "2020-02-03T22:08:07Z")

</div>

Thanks for the explanation, @lucianopaz.

I’m not sure how this relates to @junpenglao’s answer, though.

As I understood his answer, the only point at which we have to worry about getting un-transformed values is at the bottom of the theano graph.

But if I understand you correctly, that is not true – when sampling from the prior predictive, we could end up with random samples of untransformed RVs where we expect transformed ones, and those untransformed values could be used to compute values of children. If that is the case, then any use of a transformed RV, anywhere in a PyMC3 model, can cause prior predictive sampling to fail.

For _posterior_ predictive, things will be correct if the transformed values are captured in the trace generated by `pm.sample()` **But** sometimes this cannot be done for out of sample predictions, if variables’ shapes are functions of data.

If I am summarizing correctly, then any use of predictive sampling in PyMC3 is potentially risky.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 3, 2020, 10:34pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/16 "2020-02-03T22:34:05Z")

</div>

Oh, and now I think I understand @junpenglao better: The transform, as he says, acts only “backwards” from the Transformed RV up to the underlying RV in the theano graph.

TBH, I still don’t understand why it’s the Transformed RV that gets the random method instead of the Free RV – I guess that’s just the most convenient way to do things for sampling that is driven by observations.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 4, 2020, 4:08pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/17 "2020-02-04T16:08:30Z")

</div>

@lucianopaz @junpenglao – I’m left with two questions:

1. If a transformed RV does not support predictive sampling, why don’t we signal an exception when the programmer tries to sample from such a distribution? We could easily make it a special exception so that the user who understands the limitation could set a handler that accepts the un-transformed samples. I’m kind of disturbed that I have been reporting results that are inaccurate for months without realizing it.

2. Could this be fixed in `draw_values`? Could we check to see if a variable is transformed, and if so, collect the un-transformed value from the `random` method and transform it?

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [February 4, 2020, 4:27pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/18 "2020-02-04T16:27:36Z")

</div>

1. The `TransformedRV` has its random method, and it is used in predictive sampling. What’s wrong is the notion that passing a `transform` parameter to a distribution will change the values returned from said RV’s `random`.

2. It shouldn’t be `draw_values` responsibility to deal with this. We should fix the the way in which TransformedRVs are built and clearly distinguish between the two cases at hand: automatically unconstraining the values of an RV in `pm.sample` and transforming the values drawn from a given distribution. The first case is what we currently call `transform` and I think that’s wrong sand confusing. The second case is what you want to get, it would apply the transform during inference and also in predictive sampling.

---

<div class="post-metadata">

**Author:** ![rpgoldman](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/rpgoldman/32/1800_2.png) [@rpgoldman](https://discourse.pymc.io/u/rpgoldman)\
**Post date:** [February 4, 2020, 4:35pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/19 "2020-02-04T16:35:52Z")

</div>

So am I right in thinking that this is something that should be fixed by fixing the `random()` method of the transformed RV?

I assume right now we just pass the call to `random` through to a `TransformedRV`'s `distribution`, which is the un-transformed distribution. Is that correct?

Actually, I just looked and it looks like the `TruncatedNormal`, at least, has a special `random()` method that looks like it will do the right thing (at least it consults the `upper` and `lower` properties). So maybe at least _some_ of the transformed RVs work?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [February 4, 2020, 5:23pm UTC](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420/20 "2020-02-04T17:23:44Z")

</div>

I think there is still some concepts are getting confused here:

- The default transform does the right thing, because it is mapping values in the unconstrained space to the right domain
- The default random also does the right thing, as it is producing random sample regardless of what transformation in the inference time it is doing, as long as it is volume preserved.

Problem you are describing are: user specified transformation that change the domain of the distribution, for example, a sum to 1 transformation on Uniform distribution. In this case, the random method would be wrong as the transformation is modifying the domain instead of just unbounding the free parameter.

In a way, we are speaking about 2 operation here, and the current transform kwarg only deal with [1] (unbounding a constrained parameter), but not intended to solve [2] (transforming a distribution into another) although we could use it to partially solve [2].

[Next page](https://discourse.pymc.io/t/prior-predictive-sampling-with-transformed-rv/4420.md?page=2)