State space model inference by HMC


#1

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

data=pd.read_csv("https://raw.githubusercontent.com/statsmodels/statsmodels/master/statsmodels/datasets/nile/nile.csv",index_col='year',)

y=data["volume"]
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})
inference.run(n_iter=T)

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)
    plt.title(title)
    plt.show()

qmu_sample=np.array([ _qmu.sample(1000).eval()  for _qmu in qmu])
plotseries(y,qmu_sample,"Edward",1)

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