This rookie is investigating an apparent bug in the binomial distribution, adding print statements to **Binomial.random()** to understand how the binomial variables are sampled in this simple model:

```
import pymc3 as pm
import numpy as np
with pm.Model() as model:
scalar_binomial = pm.Binomial('scalar_binomial', n=1, p=0.25)
arrayed_binomial = pm.Binomial('arrayed_binomial', n=1, p=0.25, shape=10)
arrayed_binomial_alt = pm.Binomial(
'arrayed_binomial_alt', n=np.ones(10),
p=np.full(10, 0.25), shape=10)
trace = pm.sample(draws=600, chains=2, tune=500)
```

Much to my surprise, this simple model never makes a call to **Binomial.random()**; nothing is ever printed. Adding print statements to **Binomail.logp()** reveals that calls are in fact made to that function (as expected), but only during model definition, not during the call to **pm.sample()**:

How are random samples for the binomials generated, if not via **Binomial.random()** or **Binomial.logp()**?

I doubt if this question depends on the version of PyMC3, as it appears to be something basic I am not understanding. But just in case, I am running commit bf734fc2, the most recent as of today.