Convert WinBUGS model to PyMC3

Questions
1

I have WinBUGS model:

model{
	for (i in 1:n)
	{
		for(j in 1:J)
		{
			delta[i,j] <- pow(theta[i,1], Q[j,1]) * pow(theta[i,2], Q[j,2])
			delta_plus_1[i,j] <- delta[i,j] + 1
			x[i,j] ~ dbern(pi_plus_1[delta_plus_1[i,j],j])
		}
	}

	for(j in 1:J)
	{  
		pi_plus_1[1,j] ~ dbeta(3.5,23.5) 
		pi_plus_1[2,j] ~ dbeta(23.5,3.5) 

		pi_0[j] <- pi_plus_1[1,j]
		pi_1[j] <- pi_plus_1[2,j]
	}
} 

DATA

list(
n=1000, 
J=11, 
Q=structure(.Data=c(
1,0,1,0,1,1,1,0,1,1,1,0,1,0,1,0,1,0,1,0,1,0),.Dim=c(11, 2)), 
x=structure(.Data= c(
1,1,1,0,1,1,1,0,0,0,0,
0,1,1,1,1,0,0,1,0,0,0,
1,1,0,0,0,0,0,0,0,0,1,
0,1,0,1,0,1,0,0,1,0,0,
1,1,1,1,1,1,1,1,0,1,0,
1,1,1,1,1,1,0,0,0,0,0,
1,1,1,1,1,1,1,0,0,1,0,
1,1,1,0,1,0,0,0,0,1,0,
0,1,1,1,0,1,0,0,0,0,0,
1,1,1,1,0,0,1,0,0,0,0),.Dim=c(1000, 11))

I’m trying to convert this script to PyMC3:

import numpy as np
import pymc3 as pm
from theano import theano, tensor as tt
import matplotlib.pyplot as plt

seed = 314159
np.random.seed(seed)

# Data
data = [
	[1,1,1,0,1,1,1,0,0,0,0],
	[0,1,1,1,1,0,0,1,0,0,0],
	[1,1,0,0,0,0,0,0,0,0,1],
	[0,1,0,1,0,1,0,0,1,0,0],
	[1,1,1,1,1,1,1,1,0,1,0],
	[1,1,1,1,1,1,0,0,0,0,0],
	[1,1,1,1,1,1,1,0,0,1,0],
	[1,1,1,0,1,0,0,0,0,1,0],
	[0,1,1,1,0,1,0,0,0,0,0],
	[1,1,1,1,0,0,1,0,0,0,0]
]

# Constants
n = 10
J = 11
Q = [
	[1,0],
	[1,0],
	[1,1],
	[1,0],
	[1,1],
	[1,0],
	[1,0],
	[1,0],
	[1,0],
	[1,0],
	[1,0]
]

with pm.Model() as model:
	gamma_1 = pm.Beta('gamma_1', 6, 3)
	gamma_2 = pm.Beta('gamma_2', 6, 3)

	theta = pm.Bernoulli('theta', p=tt.stack([gamma_1, gamma_2]), shape=(n,2))

	pi_0 = pm.Beta('pi_0', 3.5, 23.5, shape=J)
	pi_1 = pm.Beta('pi_1', 23.5, 3.5, shape=J)

	# pData = ???

	x = pm.Bernoulli('x', p=pData, shape=(n,J), observed=data)

	# Inference
	step = pm.Metropolis([gamma_1, gamma_2, pi_0, pi_1, theta])
	trace = pm.sample(1000, step=step, cores=1, chains=1)

	print(pm.summary(trace, varnames=['gamma_1', 'gamma_2', 'pi_0', 'pi_1', 'theta']))

	pm.traceplot(trace, ['gamma_1', 'gamma_2'])
	plt.show()

I understand that pData is this piece of code in WinBUGS:

for (i in 1:n)
{
	for(j in 1:J)
	{
		delta[i,j] <- pow(theta[i,1], Q[j,1]) * pow(theta[i,2], Q[j,2])
		delta_plus_1[i,j] <- delta[i,j] + 1
		x[i,j] ~ dbern(pi_plus_1[delta_plus_1[i,j],j])
	}
}

but I don’t know how to implement it…
Please, help me.