2
$\begingroup$

My scenario is that I have two six-sided dice (D1 and D2), either of which may be fair or loaded (biased). I have samples of combined roll data (i.e. D1 + D2).

I would like to view the posterior distribution of rolls for each individual die using Markov Chain Monte Carlo (MCMC) analysis. Specifically, I would like to see if either one or both dice could potentially be loaded or unfair

What I'm unsure of is how to incorporate my combined roll data in this analysis, so that I could view the posterior distribution over the rolls of D1 and D2, and perhaps see that they are not so uniform afterall.

Can anyone point me in the right direction? Or, if perhaps I'm going about this the wrong way, I would appreciate some pointers.

Edit: I provided some beginning code below. I approximate the distribution of the combined dice to a Normal distribution after having read this .

import numpy as np
import pymc as pm
from pymc.Matplot import plot as mcplot

pm.numpy.random.seed(0)

#==============================================================================
# instantiate sample data
#==============================================================================
count_data = {2: 2, 3: 6, 4: 7, 5: 14, 6: 8, 7: 15, 8: 6, 9: 16, 10: 7, 11: 3,
              12: 1}

dice_data = list()
for index, count in count_data.iteritems():
    for i in xrange(count):
        dice_data.append(index)

dice_data = np.array(dice_data)

#==============================================================================
# create variables and model
#==============================================================================
dice1 = pm.DiscreteUniform('dice1', lower=1, upper=6)
dice2 = pm.DiscreteUniform('dice2', lower=1, upper=6)
std = pm.Uniform("std", 0, 100, trace=False)


@pm.deterministic
def mean(dice1=dice1, dice2=dice2):
    return dice1 + dice2


@pm.deterministic
def prec(U=std):
    return 1.0 / (U) ** 2

combined_roll = pm.Normal('combined_roll', mean, prec, value=dice_data, observed=True)

model = pm.Model([dice1, dice2, combined_roll, prec, std])
mcmc = pm.MCMC(model)
mcmc.sample(iter=30000, burn=10000, thin=20)

#==============================================================================
# plot results
#==============================================================================
dice1_trace = mcmc.trace('dice1')[:]
dice2_trace = mcmc.trace('dice2')[:]

mcplot(dice1_trace, 'dice1')
mcplot(dice2_trace, 'dice2')
$\endgroup$
  • $\begingroup$ Is your objective to "perform MCMC" or is it to learn something about the dice? In the latter case, what exactly do you want to learn about them? In the former case, you need to tell us more specifically what model you have in mind: "latent variables" can cover a lot of different models. $\endgroup$ – whuber Jan 25 '15 at 21:46
  • $\begingroup$ I would like to learn about the posterior distributions of both individual dice. Specifically, I would like to see if either one or both dice could potentially be loaded or unfair. In my model, I am assigning a Uniform prior distribution for both dice, although personally I believe either one or both die are biased. Using PyMC to perform inference is not a requirement but rather a learning exercise for myself. $\endgroup$ – trianta2 Jan 25 '15 at 22:03
  • 4
    $\begingroup$ If one of the dice is unfair, it would be impossible to figure out which one, as we could swap the two and get the same result. Therefore, I think you have to consider the joint posterior, as the individual posteriors won't be informative. $\endgroup$ – Sven Jan 25 '15 at 22:23
  • $\begingroup$ @Sven could you perhaps point me to some documentation which goes over this? $\endgroup$ – trianta2 Jan 25 '15 at 22:49
  • 2
    $\begingroup$ My remark has nothing to do with pyMC. All I'm saying is that you could always swap the two dice and the sum will be the same. Hence we are never able to identify the biased dice if there is one. But you know the distribution of outcomes when both are fair, so you can check whether the data supports that. If not, you know that something is wrong but you'll not be able to figure out what is wrong in all cases. $\endgroup$ – Sven Jan 26 '15 at 0:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.