Your model is unidentifiable as it is. Think of it this way: the variance in each observation in y
comes from 2 places: from the GaussianRandomWalk controlled by σ
(also the covariance with the previous time point) and the Normal distribution controlled by sd
. With uniform prior, y
could be generated completely from a GaussianRandomWalk with σ
=var(diff(y)), or a Normal distribution with sd
=var(y).
You can put more informative prior on the σ for the random walk (currently you use Exponential which put too much weight on the small value, thus pushing the variance of y
all loaded onto the Normal), but I like to instead model the total variance and distribute them onto σ
and sd
:
with pm.Model() as model:
smth_parm = pm.Uniform("alpha", lower=0, upper=1)
mu2 = pm.Normal("mu", sd=100)
tau2 = pm.Exponential("tau", 1.0/100)
z2 = pm.GaussianRandomWalk("f",
mu=mu2,
tau=tau2 / (1.0 - smth_parm),
shape=y.shape)
obs = pm.Normal("obs",
mu=z2,
tau=tau2 / smth_parm,
observed=y)
trace = pm.sample(1000, njobs=4, tune=1000)
Which variance contributed more is controlled by alpha
, you can interpreted it as a smoothing parameter.
You can find more explanations here: