[Bayesian spatiotemporal model]

Questions
1

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)
	}
	d[1:n] ~ car.normal(adj[],weights[],num[],taud)
	U[1:n] ~ car.normal(adj[],weights[],num[],omega)
	
	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)

adj = np.array([[11,8,7],
[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]])

for i in range(len(adj)):
for j in range(len(adj[i])):
adj[i][j] = adj[i][j]-1
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’)

for i, a in enumerate(adj):
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))

trace2 = pm.sample(1000, tune=500, cores=4,init='advi')
2

I’m not trying to be rude, but I think you should reformat your post so your code looks cleaner. You gave headings to some parts which is good, but there is a lot going on here, especially with some code formatted in the special way, and other parts not.

3

Sorry, I’ve re edited it~~