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()
```

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)
}
```

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)