porting pymc2 bayesian a/b testing example to pymc3

You were indeed close. However, this line:

obs = pm.Bernoulli('obs', p_true, observed = occurrences)

is wrong as you are just setting a constant value for p (p_true == 0.05). Thus, your random variable p defined above to have a uniform prior is not constrained by the likelihood and your plot shows that you are just sampling from the prior. If you replace p_true with p in your code it should work. Here is the fixed version:

import pymc as pm

import random
import numpy as np
import matplotlib.pyplot as plt

with pm.Model() as model:
   # Prior is uniform: all cases are equally likely
p = pm.Uniform('p', lower = 0, upper = 1)

# set constants
p_true = 0.05 # remember, this is unknown.
N = 1500

# sample N Bernoulli random variables from Ber(0.05).
# each random variable has a 0.05 chance of being a 1.
# this is the data - generation step
occurrences = [] # pm.rbernoulli(p_true, N)
for i in xrange(N):
   occurrences.append((random.uniform(0.0, 1.0) <= p_true))
occurrences = np.array(occurrences)
obs = pm.Bernoulli('obs', p, observed = occurrences)

start = pm.find_MAP()
step = pm.Metropolis()
trace = pm.sample(18000, step, start)

pm.traceplot(trace);

This worked for me. I generated the observations before initiating the model.

true_p_A = 0.05
true_p_B = 0.04
N_A = 1500
N_B = 750

obs_A = np.random.binomial(1, true_p_A, size = N_A)
obs_B = np.random.binomial(1, true_p_B, size = N_B)

with pm.Model() as ab_model:
   p_A = pm.Uniform('p_A', 0, 1)
p_B = pm.Uniform('p_B', 0, 1)
delta = pm.Deterministic('delta', p_A - p_B)
obs_A = pm.Bernoulli('obs_A', p_A, observed = obs_A)
osb_B = pm.Bernoulli('obs_B', p_B, observed = obs_B)

with ab_model:
   trace = pm.sample(2000)

pm.traceplot(trace)

You were very close - you just need to unindent the last two lines, which produce the traceplot. You can think of plotting the traceplot as a diagnostic that should occur after you finish sampling. The following works for me:

import pymc as pm

import random
import numpy as np
import matplotlib.pyplot as plt

with pm.Model() as model:
   # Prior is uniform: all cases are equally likely
p = pm.Uniform('p', lower = 0, upper = 1)

# set constants
p_true = 0.05 # remember, this is unknown.
N = 1500

# sample N Bernoulli random variables from Ber(0.05).
# each random variable has a 0.05 chance of being a 1.
# this is the data - generation step
occurrences = [] # pm.rbernoulli(p_true, N)
for i in xrange(N):
   occurrences.append((random.uniform(0.0, 1.0) <= p_true))
occurrences = np.array(occurrences)
obs = pm.Bernoulli('obs', p_true, observed = occurrences)

start = pm.find_MAP()
step = pm.Metropolis()
trace = pm.sample(18000, step, start)

#Now plot
pm.traceplot(trace)
plt.show()

Suggestion : 2

A walkthrough of implementing a Conditional Autoregressive (CAR) model in PyMC3, with WinBugs/PyMC2 and STAN code as references.,So lets try to implement the CAR model using theano.scan. First we create a theano function with theano.scan and check if it really works by comparing its result to the for-loop.,Another thing to keep in mind is that sometimes your model is sensitive to prior choice: for example, you can have a bad experiences using Normals with a large sd as prior. Gelman often recommends Cauchy or StudentT, and more heuristic on prior could be found on the Stan wiki.,Notice that since the model parameterization is different than in the WinBugs model, the alpha doesn’t bear the same interpretation.

In[1]:

% pylab inline

import numpy as np
import scipy.stats as stats
import pymc3 as pm
from theano
import shared
import theano
import theano.tensor as tt
floatX = "float32"

   %
   config InlineBackend.figure_format = 'retina'
plt.style.use('ggplot')

Populating the interactive namespace from numpy and matplotlib

In[2]:

model {
   for (i in 1: regions) {
      O[i] ~dpois(mu[i])
      log(mu[i]) < -log(E[i]) + beta0 + beta1 * aff[i] / 10 + phi[i] + theta[i]
      theta[i] ~dnorm(0.0, tau.h)
   }
   phi[1: regions] ~car.normal(adj[], weights[], Wplus[], tau.c)

   beta0~dnorm(0.0, 1.0E-5) # vague prior on grand intercept
   beta1~dnorm(0.0, 1.0E-5) # vague prior on covariate effect

   tau.h~dgamma(3.2761, 1.81)
   tau.c~dgamma(1.0, 1.0)

   sd.h < -sd(theta[]) # marginal SD of heterogeneity effects
   sd.c < -sd(phi[]) # marginal SD of clustering(spatial) effects

   alpha < -sd.c / (sd.h + sd.c)
}

Suggestion : 3

2022-01-29

% matplotlib inline
   %
   config InlineBackend.figure_format = 'retina'

from warnings
import filterwarnings

from aesara
import pprint
from matplotlib
import pyplot as plt, ticker
import numpy as np
import pandas as pd
import scipy as sp
import seaborn as sns
from sklearn.preprocessing
import StandardScaler

filterwarnings('ignore', category = RuntimeWarning, message = "overflow encountered in exp")
filterwarnings('ignore', category = UserWarning, module = 'pymc',
   message = "Unable to validate shapes: Cannot sample from flat variable")

FIG_SIZE = np.array([8, 6])
plt.rc('figure', figsize = FIG_SIZE)

dollar_formatter = ticker.StrMethodFormatter("${x:,.2f}")
pct_formatter = ticker.StrMethodFormatter("{x:.1%}")

sns.set(color_codes = True)

SEED = 123456789 #
for reproducibility

rng = np.random.default_rng(SEED)