Sorry if this is a very easy solution - I searched for a few hours and wasn’t able to find a solution that I understood.
I was wondering if there was a way to introduce uncertainty in a parameter that is being used in the model, without adding it to the chain. For example, suppose I have a linear relationship between two parameters x and y, such that
y = Ax + B, where A = -0.006 +/- 0.001 and B = 1 +/- 0.04. I am sampling x and want to change A and B each chain sample so I can incorporate the uncertainty in my input parameters into the model uncertainty. Is there a way to do this? Thanks.
EDIT: Sorry all, it looks like my phrasing above was confusing. In pseudo-code, here’s what I want:
import numpy as np
import pymc3 as pm
# Dummy data
Y = np.arange(0, 50) + np.random.randn(50)
for sample in int(nburnin + nsamples):
with pm.Model() as model:
A = np.random.normal(-0.006, 0.001)
B = np.random.normal(1.0, 0.04)
y_true = pm.Data('y', Y)
x = pm.Flat('x')
y_syn = A*x + B
#... augment likelihood to look at difference between y_true and y_syn ...
trace = pm.sample(...)
Of course the above code wouldn’t work the way it is intended with PyMC3, but this demonstrates my point visually, I think.
I don’t want to use pm.Normal('A', mu=-0.006, sigma=0.001 and pm.Normal('B', mu=1, sigma=0.04) because (my understanding) is that this would sample those parameters to solve them in the MCMC. I just want the values and their uncertainties added into my model solving for x.
it depends on your statistical model. for example, if A is normally distributed, uncertainty is the ‘sigma’ parameter. (Note: I’m not sure if I understand your question properly!)
Not that I know of. But note that 1 chain is just 1 sampling of the model, so you could just perturb A and B between calls to pm.sample something like
for j in range(nchains):
A += np.random.(...) # I *think* += operates on the underlying memory rather than creating a new object
B += np.random.(...)
x_chains.append(pm.sample(..., chains=1))
I think this is the closest thing to what I need, thanks. I also had the idea of perturbing the A & B each chain, but I thought I’d ask if there was a way to do it each sample.
Sorry for the late reply. Wouldn’t using A_noisy = pm.Normal("A_noisy", A, 1) allow the sampling to change the distribution of A over the chain? I would want my posterior distribution for A to be the same as its prior.
I am not sure what your goal is. You want a prior that’s fixed but random? The way I would think about it would be to have fixed parameters (if the approximation is not too bad). So you don’t have a fixed A with a special pdf but a distribution with fixed parameters. Say A_mean = 10 and A_std = 1. Then A_noisy = pm.Normal("A_noisy", A_mean, A_std). The latter would change obviously during inference but A itself would be fixed as those are just numbers. Then the prior for A_noisy2 = pm.Normal("A_noisy", A_mean, A_std) in another context is exactly the same as A_noisy1 was. It’s a different way of thinking about it but I feel it achieves the same result, no? If you don’t like to force a Normal prior you can use something else.
Basically, I have two parameters A and B that have known value and uncertainty. I don’t need to sample in the MCMC, since their value is already determined. I would, however, like to incorporate their uncertainties into my model. Ideally, I would accomplish this by changing the A and B parameters each sample, but not updating these parameters each step.
I think I understand the idea (although not the rationale). I really don’t know how wise it is to change things that pymc expects to be fixed during sampling.
I still feel that what you want is equivalent in practice to just resetting your prior of A every time you sample, but I might be missing something. How would the generative model for your data look like?
Actually I think you’re correct - resetting the prior each sample would achieve the same thing. I don’t see how to do that using pymc3.sample(), however.