# State space model inference by HMC

**URL:** <https://discourse.edwardlib.org/t/state-space-model-inference-by-hmc/764>\
**Category:** General\
**Created:** [April 27, 2018, 1:04am UTC](https://discourse.edwardlib.org/t/state-space-model-inference-by-hmc/764 "2018-04-27T01:04:21Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![xiangze](https://avatars.discourse-cdn.com/v4/letter/x/fbc32d/32.png) [@xiangze](https://discourse.edwardlib.org/u/xiangze)\
**Post date:** [April 27, 2018, 1:04am UTC](https://discourse.edwardlib.org/t/state-space-model-inference-by-hmc/764/1 "2018-04-27T01:04:21Z")

</div>

I tried state space model by using Other guy’s code  
([http://www.ie110704.net/2018/04/21/edward、pystan、pymc3で状態空間モデルを実装してみた/](http://www.ie110704.net/2018/04/21/edward%E3%80%81pystan%E3%80%81pymc3%E3%81%A7%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%82%92%E5%AE%9F%E8%A3%85%E3%81%97%E3%81%A6%E3%81%BF%E3%81%9F/) in Japanese).

```auto
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  
 ![10](https://cdck-file-uploads-canada1.s3.dualstack.ca-central-1.amazonaws.com/flex035/uploads/edwardlib/original/1X/2c7689be485e10dac170e03b18751ba43460e95b.jpg)  
[http://nbviewer.jupyter.org/gist/xiangze/b88f8bba67db598d45ab51bfdc6323bc#Edward](http://nbviewer.jupyter.org/gist/xiangze/b88f8bba67db598d45ab51bfdc6323bc#Edward)

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

[http://nbviewer.jupyter.org/gist/xiangze/b88f8bba67db598d45ab51bfdc6323bc#Stan](http://nbviewer.jupyter.org/gist/xiangze/b88f8bba67db598d45ab51bfdc6323bc#Stan)  
priors are different.

Is there any problems in this code or HMC?

Similar problem  
[https://discourse.edwardlib.org/t/hmc-and-klpq-fail-on-edward-while-pymc3-works-fine/740](https://discourse.edwardlib.org/t/hmc-and-klpq-fail-on-edward-while-pymc3-works-fine/740)