[Bayesian spatiotemporal model]

I am currently making a Bayesian space-time model. I have already written a winbugs model. Now I have converted it to pymc3 code. Currently pymc3 can run, but I do n’t know if this conversion is correct. Please help me take a look；

Winbugs

``````model {
for (i in 1:n) {
for (j in 1:58) {
Y[i,j] ~ dnorm(mu[i,j],sigma2)
mu[i,j] <- alpha0 + alpha1*tf[j] + U[i] + V[i] + g[j] + locatime[i,j] + d[i]*tf[j] + beta1*prs[i,j] + beta2*tem[i,j]
}
locatime[i,1] ~ dnorm(0,taupsi)
for (j in 2:58){
locatime[i,j] ~ dnorm(locatime[i,j-1],taupsi) }
V[i] ~ dnorm(0,tauv)
}

for (k in 1:sumNumNeigh){
weights[k] <- 1
}
g[1] ~ dnorm(0,0.0001)
tf[1] <- time[1] - 1959
for (j in 2:58){
g[j] ~ dnorm(g[j-1], tau.g)
tf[j] <- time[j] - 1959}

alpha0 ~ dflat()
alpha1 ~ dnorm(0,0.0001)
omega ~ dgamma(0.5,0.0005)
tauv <- 1/pow(sdv,2)
taupsi <- 1/pow(sdpsi,2)
taud <- 1/pow(sdd,2)
tau.g <- 1/pow(sdg,2)
sdg ~ dunif(0,10)
sdpsi ~ dunif(0,10)
sdd ~ dunif(0,10)
sdv ~ dunif(0,10)
sigma ~ dunif(0,10)
sigma2 <- sigma*sigma
beta1 ~ dnorm(0.0, tau.beta1)
tau.beta1 <- 1/pow(sdbeta1,2)
sdbeta1 ~ dunif(0,10)
beta2 ~ dnorm(0.0, tau.beta2)
tau.beta2 <- 1/pow(sdbeta2,2)
sdbeta2 ~ dunif(0,10)
}
``````

pymc3

def MaxMinNormalization(x):
“”"[0,1] normaliaztion"""
x = (x - np.min(x)) / (np.max(x) - np.min(x))
return x

rx1 = MaxMinNormalization(rx1)
prs = MaxMinNormalization(prs)
tem = MaxMinNormalization(tem)

N = len(rx1)

[11],
[7,4],
[5,6,3],
[6,4],
[4,5],
[3,9,10,1,8],
[1,10,7],
[7,10],
[7,8,9],
[1,2]])

weights = np.array([[1,1,1],
[1],
[1,1],
[1,1,1],
[1,1],
[1,1],
[1,1,1,1,1],
[1,1,1],
[1,1],
[1,1,1],
[1,1]])

Wplus = np.asarray([sum(w) for w in weights])
value = np.asarray(np.random.randn(N,), dtype=theano.config.floatX)

maxwz = max([sum(w) for w in weights])
N = len(weights)
wmat = np.zeros((N, N))
amat = np.zeros((N, N), dtype=‘int32’)

amat[i, a] = 1
wmat[i, a] = weights[i]

value = np.asarray(np.random.randn(N,), dtype=theano.config.floatX)

print(np.sum(value*amat, axis=1)/np.sum(wmat, axis=1))

def mu_phi(value):
N = len(weights)
# Calculate mu based on average of neighbours
mu = np.array([np.sum(weights[i]*value[adj[i]])/Wplus[i] for i in range(N)])
return mu
sdpsi = np.random.uniform(0,10)
taupsi = 1/pow(sdpsi,2)
locatime = np.zeros((58,))
locatime[0] = np.random.normal(0,taupsi,1)
for j in range(1,58):
locatime[j] = np.random.normal(locatime[j-1],taupsi,1)
sdg = np.random.uniform(0,10)
taug = 1/pow(sdg,2)
g = np.zeros((58,))
g[0] = np.random.normal(0,0.0001,1)
for j in range(1,58):
g[j] = np.random.normal(g[j-1],taug,1)
tf = list(range(1,59))

with pm.Model() as model:

``````sdbeta1 = pm.Uniform('sdbeta1',0,10)
sdbeta2 = pm.Uniform('sdbeta2',0,10)
sigma1 = pm.Uniform('sigma1',0,10)
sigma2 = sigma1**2
alpha0 = pm.Normal('alpha0', mu=0.2, tau=0.07)
alpha1 = pm.Normal('alpha1', mu=0.4, tau=0.048)
taubeta1 = 1/pow(sdbeta1,2)
taubeta2 = 1/pow(sdbeta2,2)
tau_u = pm.Gamma('tau_u', alpha=2.5, beta=1.4)
tau_d = pm.Gamma('tau_d', alpha=2.5, beta=1.4)
tau_v = pm.Gamma('tau_v', alpha=3.2761, beta=1.81)
mu_U = CAR.CAR2('mu_U', w=wmat, a=amat, tau=tau_u, shape=N)
mu_d = CAR.CAR2('mu_d', w=wmat, a=amat, tau=tau_d, shape=N)

V = pm.Normal('V', mu=0.0, tau=tau_v, shape=N)
U = pm.Deterministic('U', mu_U-tt.mean(mu_U))

d= pm.Deterministic('d', mu_d-tt.mean(mu_d))
beta0 = pm.Normal('beta0', mu=0.0, tau=taubeta1)
beta1 = pm.Normal('beta1', mu=0.0, tau=taubeta2)
mu = pm.Deterministic('mu', tt.exp(alpha0 + alpha1*tf[j] + U[i] + V[i] + g[j]  + d[i]*tf[j] + beta0*prs.iloc[i,j] + beta1*tem.iloc[i,j] ))   ###  I don't know if it's correct when I using g[j] indicate timeseries variables;;

Yi = pm.Normal('Yi', mu=mu,sigma=sigma2,observed=rx1,shape=(11,58))