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
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
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