Convert WINBUGS model to PyMC3

Question

I am currently taking a class on Bayesian statistics, we are allowed to use any package to computationally solve the models but all of the examples are provided in WINBUGS. I would prefer to use Python and PyMC3. I don't have much experience using PyMC3 and could use some help on how to convert this simple WINBUGS model into a PyMC3 model.

The example WINBUGS code is below. It is a simple Binomial comparing two options with a different number of observations per sample. The model also tests 5 different priors.

model{
for(i in 1:5){
   n1[i] <- Tot1   #100
   n2[i] <- Tot2   #3
   y1[i] <- Positives1
   y2[i] <- Positives2
   y1[i] ~ dbin(p1[i],n1[i])
   y2[i] ~ dbin(p2[i],n2[i])
  diffps[i] <- p1[i]-p2[i]   #100seller - 3seller
}

# Uniform priors
   p1[1] ~ dbeta(1, 1); p2[1] ~ dbeta(1, 1)

# Jeffreys'  priors
   p1[2] ~ dbeta(0.5,  0.5);  p2[2] ~ dbeta(0.5, 0.5)

# Informative priors centered at about 93% and 97%
   p1[3] ~ dbeta(30,2);  p2[3] ~ dbeta(2.9,0.1)

# Zellner priors prop to 1/(p * (1-p))
  logit(p1[4]) <- x[1]
    x[1] ~ dunif(-10000, 10000)  # as  dflat()
   logit(p2[4]) <- x[2]
   x[2] ~ dunif(-10000, 10000)   # as dflat()

#Logit centered at 3 gives mean probs close to 95%
  logit(p1[5]) <- x[3]
    x[3] ~  dnorm(3, 1)  
   logit(p2[5]) <- x[4]
    x[4] ~ dnorm(3, 1) 
}

DATA
list(Tot1=100,  Tot2=3, Positives1=95, Positives2=3)

INITS
list(p1=c(0.9, 0.9, 0.9, NA,NA), 
p2=c(0.9, 0.9, 0.9, NA,NA), x=c(0,0,0,0))

list(p1=c(0.5, 0.5, 0.5, NA,NA), 
p2=c(0.5, 0.5, 0.5, NA,NA), x=c(0,0,0,0))

list(p1=c(0.3, 0.3, 0.3, NA,NA), 
p2=c(0.3, 0.3, 0.3, NA,NA), x=c(0,0,0,0))

In PyMC3 I attempted to implement the first of the 5 priors on a single sample (I am not sure how to do both) with the following code:

import np as np
import pymc3 as pm

sample2 = np.ones(3)

with pm.Model() as ebay_example:
    prior = pm.Beta('theta', alpha = 1, beta = 1)
    likelihood = pm.Bernoulli('y', p = prior, observed = sample2)
    trace = pm.sample(1000, tune = 2000, target_accept = 0.95)

The above model ran, but the results don't not align with the BUGS results, I am not sure if it because I didn't do a burn in or some other larger issue. Any guidance would be great.

Answer 1

we are taking the same class right now. Following are my codes and the difference of means is close to BUGS results.

na = 100
nb = 3
pos_a = 95
pos_b = 3

with pm.Model() as model:
    # priors
    p0a = pm.Beta('p0a', 1, 1)

    # likelihood
    obs_a = pm.Binomial("obs_a", n=na, p=p0a, observed=pos_a)

    # sample
    trace1_a = pm.sample(1000)


with pm.Model() as model:
    # priors
    p0b = pm.Beta('p0b', 1, 1)

    # likelihood
    obs_b = pm.Binomial("obs_b", n=nb, p=p0b, observed=pos_b)

    # sample
    trace1_b = pm.sample(1000)

pm.summary(trace1_a)["mean"][0] - pm.summary(trace1_b)["mean"][0]
OUT: 0.1409999999999999