How to (re-)sample from an observed random (mixture) variable

junpenglao · April 9, 2018, 10:32am

I tend to trust sampling a bit more as well. There really not enough theory around automatic VI yet unfortunately.

I think including generative RV in the model is currently unsafe in PyMC3. In theory it should not change the inference as the generative node is conditional on the other parent parameters and has no effect on the overall model inference. But currently it does (see a discussion here):

github.com/pymc-devs/pymc

Weird turbulence cause by generative discrete RV

opened 11:34AM - 04 Apr 17 UTC

closed 12:18PM - 22 Dec 18 UTC

junpenglao

In `Stan`, there is an option to write a `generated quantities` block for sample… generation. Doing the similar in pymc3, however, seems to introduce weird turbulence to the sampler, especially if the generated RV is discrete. Consider the following simple sample: ```python import numpy as np import pymc3 as pm import theano.tensor as tt # Data x = np.array([1.1, 1.9, 2.3, 1.8]) n = len(x) sigmoid = lambda x: 1/(1 + tt.exp(-x)) np.random.seed(42) with pm.Model() as model1: # prior mu = pm.Normal('mu', mu=0, tau=.001) sigma = pm.Uniform('sigma', lower=0, upper=10) # observed xi = pm.Normal('xi', mu=mu, tau=1/(sigma**2), observed=x) # inference trace = pm.sample(1000, njobs=5, tune=500, init=None) pm.traceplot(trace, varnames=['mu', 'sigma_interval_']) ``` The result is fairly normal: ``` Assigned NUTS to mu Assigned NUTS to sigma_interval_ 100%|██████████| 1000/1000 [00:01<00:00, 855.81it/s] ``` ![figure_1](https://cloud.githubusercontent.com/assets/12952641/24654452/f6234bdc-1939-11e7-8dd7-e6a671254d7c.png) If now I add an addition RV node to the graph: ```python np.random.seed(42) with pm.Model() as model1: # prior mu = pm.Normal('mu', mu=0, tau=.001) sigma = pm.Uniform('sigma', lower=0, upper=10) # observed xi = pm.Normal('xi', mu=mu, tau=1/(sigma**2), observed=x) # generation p = pm.Deterministic('p', sigmoid(mu)) count = pm.Binomial('count', n=10, p=p, shape=10) # inference trace = pm.sample(1000, njobs=5, tune=500, init=None) ``` The output trace is quite unstable and converge much slower: ``` Assigned NUTS to mu Assigned NUTS to sigma_interval_ Assigned Metropolis to count 100%|██████████| 1000/1000 [00:06<00:00, 153.00it/s] ``` ![figure_2](https://cloud.githubusercontent.com/assets/12952641/24654564/76a708f2-193a-11e7-832a-d997dc0fb1a0.png) if the added RV is continuous, the effect seems to be minimal: ```python np.random.seed(42) with pm.Model() as model1: # prior mu = pm.Normal('mu', mu=0, tau=.001) sigma = pm.Uniform('sigma', lower=0, upper=10) # observed xi = pm.Normal('xi', mu=mu, tau=1/(sigma**2), observed=x) # generation p = pm.Deterministic('p', sigmoid(mu)) xi_ = pm.Normal('xi_', mu=p, shape=10) # inference trace = pm.sample(1000, njobs=5, tune=500, init=None) ``` ``` Assigned NUTS to mu Assigned NUTS to sigma_interval_ Assigned NUTS to xi_ 100%|██████████| 1000/1000 [00:07<00:00, 128.04it/s] ``` ![figure_3](https://cloud.githubusercontent.com/assets/12952641/24654624/c6c91730-193a-11e7-9c25-16e3cd4e814b.png)

For sampling, the problem is that the auto compound step does not fully follow the order of the graph, which means you might sometimes end up sampling the generative RV before sampling its parent, which generated a bias because now the parent parameters are depending on the sample of the generative RV. But even when the generative RV is sampled after its parent, it still has an effect on the sampling in the next sample, as the model.logp includes the generative RV and thus the MH update is dependent on the value of the generative RV.

It might work in some model, overall you should always added the generative node after inference, using sample_ppc, eg

Topic		Replies	Views
4 lines basic mixture model outputs wrong results Questions	0	383	March 12, 2020
Manually sampling from prior predictive Questions	10	1001	July 4, 2022
Sample from posterior predictive after fixing latent varaibles Questions	1	574	October 23, 2018
Sampling from a model where the observed data are numerically transformed from a RV Questions	1	412	November 24, 2022
Mixture with multiple observations Questions	6	2294	December 14, 2017

How to (re-)sample from an observed random (mixture) variable

Related topics