Prevent prior from updating?

Hi all,

I am currently setting up a model in which I intend to use a normally distributed parameter which interacts with other parameters to come to certain observed results.

However, for this known normally distributed parameter, I would like to be able to set a prior that does not update during sampling, but just samples from the initially set distribution for this prior. I don’t have observed samples from this prior but could easily generate some, but then I run into shape mismatch errors. Is there no simple option where you could just set observed=True or sth similar which shows pymc3 that this should be an observed variable rather than a FreeRV?

So you want a stochastic node in your model? Check out http://deeplearning.net/software/theano/library/tensor/shared_randomstreams.html

Hi @junpenglao ,

Thank you for your response.
I would like to set a parameter used in my model to a certain normal distribution with known mu and sigma, which should not change due to the observed data outcome. Now this tends to happen if I naively set this prior within the pymc3 model I set up.

Are there any examples where someone uses theano shared random streams to do this within a pymc3 model?

Thanks!

A simple example goes like this:

import numpy as np
import pymc3 as pm
import theano.tensor as tt

true_mu, true_sigma = 5., 2.
y_obs = np.random.randn(50) * true_sigma + true_mu

with pm.Model() as m:
    mu = pm.Normal('mu', 0., 100.)
    sigma = pm.HalfCauchy('sigma', 5.)
    y = pm.Normal('y', mu, sigma, observed=y_obs)
    trace = pm.sample()

pm.summary(trace)

srng = tt.shared_randomstreams.RandomStreams(seed=234)

with pm.Model() as m:
    mu = pm.Deterministic('mu', srng.normal(avg=4.9, std=.1))
    sigma = pm.HalfCauchy('sigma', 5.)
    y = pm.Normal('y', mu, sigma, observed=y_obs)
    trace = pm.sample()

pm.summary(trace)

In the second model, there is only 1 free parameter sigma, and mu always follows the same distribution.

3 Likes

Building off of this question - how would I modify if, instead of using a known normal distribution as a stochastic node (as in the case of mu above), I wanted to use some posterior traces from some other model? e.g. if I had a trace containing mu, but it is not necessarily normally distributed for conversion to a normal dist.

I’ve read this tutorial: https://docs.pymc.io/notebooks/updating_priors.html for creating priors from pm.Interpolated() of a kde of the posterior, which addresses creation of an RV from arbitrarily distributed traces, but I’m not sure how to make it do so in a way that it isn’t a free parameter in the model.

Thanks in advance.

@Boilermaker14 Have you had any chance with that?

I used that approach successfully for a categorical distribution. However, predictive prior sampling does not work as expected. Can someone explain why?

Below an example of what I am doing and why. I have 1000 samples from a random variable from a previous sampling (representing global-mean sea level), which needs to be combined with other parameters (representing local effects at the coast). These local parameters need to undergo a Bayesian update, but the original ensemble should stay put, i.e. conserving its original statistical properties.

Here is what I started from:

# generate integer to index a pre-existing 1000-member ensemble
ensemble_i = pm.Categorical("global", p=np.ones(1000)/1000) 

which is now replaced by the following, following the advice in this post:

srng = tt.random.utils.RandomStream(seed=234)
ensemble_i = pm.Deterministic("global", srng.categorical(np.ones(1000)/1000)) 

The sampling then works beautifully (takes 16 seconds instead of 4 seconds, but that’s OK). However, I was surprised to find that prior predictive sampling does not work as expected any more:

with model:
    prior = pm.sample_prior_predictive(samples=1000)

Where I used to have an array for my “global” variable, I now have a scalar !
However, I could verify in the posterior that “global” is sampled from 0 to 999, as expected. The easy workaround is to have two versions of the model, one with “fixed” R.V. to generate the posterior, and one with “free” R.V. to generate the prior, but I’d love to understand this weird behavior.

Does it make sense to anyone?