State space model inference by HMC


I tried state space model by using Other guy’s code
(、pystan、pymc3で状態空間モデルを実装してみた/ in Japanese).


T = 1000
N = len(y)

muZero = Normal(loc=0.0, scale=1.0)
sigmaW = InverseGamma(concentration=1.0, rate=1.0)
mu = [0]*N
mu[0] = Normal(loc=muZero, scale=sigmaW)
for n in range(1, N):
    mu[n] = Normal(loc=mu[n-1], scale=sigmaW)
sigmaV = InverseGamma(concentration=1.0, rate=1.0)
y_pre   = Normal(loc=tf.stack(mu), scale=sigmaV)
qmuZero = Empirical(tf.Variable(tf.zeros(T)))
qsigmaW = Empirical(tf.Variable(tf.zeros(T)))
qmu = [Empirical(tf.Variable(tf.zeros(T))) for n in range(N)]
qsigmaV = Empirical(tf.Variable(tf.zeros(T)))
latent_vars = {m: qm for m, qm in zip(mu, qmu)}
latent_vars[muZero] = qmuZero
latent_vars[sigmaW] = qsigmaW
latent_vars[sigmaV] = qsigmaV

inference = ed.HMC(latent_vars, data={y_pre: y.values})

def plotseries(y,mu_samples,title,ax=0):
    #axis=0 stan, 1 pymc
    mu_lower, mu_upper = np.percentile(mu_samples, q=[2.5, 97.5], axis=ax)
    pred = mu_samples.mean(axis=ax)
    plt.plot(list(range(len(y)+1)), list(y)+[None], color='black')
    plt.plot(list(range(1, len(y)+1)), pred)
    plt.fill_between(list(range(1, len(y)+1)), mu_lower, mu_upper, alpha=0.3)

qmu_sample=np.array([ _qmu.sample(1000).eval()  for _qmu in qmu])

In edward 1.3.4 the result is

Same code in edward 1.3.5 outputs all 0.
As a reference stan’s result is the following one.

priors are different.

Is there any problems in this code or HMC?

Similar problem