'How to Simulate a Biased 6-sided Dice using Pymc3?

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?



Solution 1:[1]

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!

Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source
Solution 1 Luxspes