0

Say I have 10 coins from the same mint, I flip them each 50 times, now I want to estimate bias of the mint as well as the individual bias of all the coins.

The way I want to do this is like this:

# Generate a list of 10 arrays with 50 flips in each
test = [bernoulli.rvs(0.5, size=50) for x in range(10)]

with pm.Model() as test_model:
    k = pm.Gamma('k', 0.01, 0.01) + 2
    w = pm.Beta('w', 1, 1)

thetas = pm.Beta('thetas', w * (k - 2) + 1, (1 - w) * (k - 2) + 1, shape = len(test))
y = pm.Bernoulli('y', thetas, observed=test)

But this doesn't work, because now it seems like pymc expects 50 coins with 10 flips. I can hack around this issue in this instance. But, I'm both a beginner at python and pymc(3) so I want to learn why it behaves like this and what a proper simulation of this situation should look like.

Jon Sjöberg
  • 159
  • 1
  • 7

1 Answers1

1

If you are new to Python may be you are not familiar with the concept of broadcasting, that is used when working with NumPy arrays, and is also useful for defining PyMC3 models. Broadcasting enables us to operate arithmetically with arrays of difference size under certain circumstances.

For your particular example the problem is that according to the broadcasting rules the shape of the data vector and the shape of the thetas vector are not compatible. The easiest solution for your problem is to transpose the data vector (make rows columns and columns rows). Notice also that using SciPy you can create you mock data without using a list comprehension, you just need to pass the proper shape.

test = bernoulli.rvs(0.5, size=(50, 10))


with pm.Model() as test_model:
    k = pm.Gamma('k', 0.01, 0.01) + 2
    w = pm.Beta('w', 1, 1)
    thetas = pm.Beta('thetas', w * (k - 2) + 1, (1 - w) * (k - 2) + 1, shape = test.shape[1])
    y = pm.Bernoulli('y', thetas, observed=test)
aloctavodia
  • 2,040
  • 21
  • 28