3

How do I simulate a 6-side Dice roll using Pymc3? Also, what is I know that different sides of the dice have different distributions?

ashish trehan
  • 413
  • 1
  • 5
  • 9
  • Please consider accepting an answer: https://stackoverflow.com/help/someone-answers – merv Dec 10 '19 at 17:03

1 Answers1

10

The easiest way to simulate 1000 rolls of a fair 6-sided die in PyMC3 is

import pymc3 as pm

with pm.Model():
    rolls = pm.DiscreteUniform('rolls', lower=1, upper=6)
    trace = pm.sample(1000)
trace['rolls']  # shows you the result of 1000 rolls

Note that this is slower, but equivalent, to just calling np.random.randint(1, 7, size=1000).

For 1000 rolls of an unfair die

probs = np.array([0.1, 0.2, 0.3, 0.2, 0.1, 0.1])

with pm.Model():
    rolls = pm.Multinomial('rolls', n=1000, p=probs, shape=6)
    trace = pm.sample(1)

Which is again equivalent, but slower, than np.random.multinomial(1000, pval=probs).

The situtation in which you would want to use PyMC3 is if you observe, say, 50 rolls of an unfair die, have some prior expectation that it is a fair die, and want to evaluate the posterior of that expectation. Here's an example of that:

observations = np.array([20, 6, 6, 6, 6, 6]) # sums up to 50
with pm.Model():
    probs = pm.Dirichlet('probs', a=np.ones(6))  # flat prior
    rolls = pm.Multinomial('rolls', n=np.sum(observations), p=probs, observed=observations)
    trace = pm.sample(1000)
trace['probs']  # posterior samples of how fair the die are

You can use the built-in traceplot to see how the samples look:

posterior plot

Note that we correctly work out that one of the sides comes up more often than the others!

Luxspes
  • 6,268
  • 2
  • 28
  • 31
colcarroll
  • 3,632
  • 17
  • 25