I deleted this since I though I understood my mistake, but now I’m confused again. It behaves as expected if you change the code to
import numpy as np
import pymc3 as pm
X = np.random.randint(0,2,10000)
Y = np.random.randint(0,2,10000)
# setup model
model = pm.Model()
with model:
w = pm.Normal('weight', mu = 0, sd = 1)
pm.Potential('coupling', w*(2*X-1)*(2*Y-1)/10000)
trace = pm.sample()
pm.traceplot(trace)
The only change made is that I divided the coupling by the number of datapoints. I am confused as to why this is correct. Can someone help?