4
$\begingroup$

Suppose you have a die, but you do not know the probabilities of the individual sides. How can you determine the individual probabilities by rolling the dice repeatedly and independently?

$\endgroup$
3
  • 3
    $\begingroup$ So, you have samples from a probability distribution. What is the simplest thing you could do with them? $\endgroup$
    – gunes
    Jan 4, 2021 at 21:24
  • $\begingroup$ Read about the "frequentist" definition of probability and you'll be on your way to figuring out the probabilities $\endgroup$ Jan 4, 2021 at 21:35
  • 1
    $\begingroup$ With enough data from an unfair die, a chi-squared goodness-of-fit test will show its unfairness. (My answer seems to have survived Google translation to German, except for messing up the math notation, as usual--and once 'Broetchen' for rolls of the die and once 'Sterben' for die. Clearly a Google work in progress.) $\endgroup$
    – BruceET
    Jan 5, 2021 at 7:06

3 Answers 3

4
$\begingroup$

Throw it many times and draw the histogram. The more times you throw it the more accurately (law of the large numbers). The resulting shape will approximate the probability mass function describing the random variable dice-side.

This is at the core of Monte carlo methods which use randomness to approximate a stochastic phenomenon.

$\endgroup$
2
$\begingroup$

Consider a $k$-sided die and suppose we observe the counts $\mathbf{n} \equiv (n_1,...,n_k)$ with $n = \sum n_i$ total rolls. Suppose that the sequence of die-rolls is exchangeable, so that the counts of the outcome come from a multinomial distribution with probability vector $\mathbf{p} \equiv (p_1,...,p_k)$. Under these conditions we have $\mathbf{n}/n \rightarrow \mathbf{p}$ (from the law of large numbers) and so the way we determine the probabilities is to roll the die a large number of times and take the sample proportions as estimates of the probabilities.

If you would like to do this with some formal statistical modelling, this inference problem is handled quite well by using Bayesian analysis with a Dirichlet prior for the true probabilities for your die, which leads to the Multinomial-Dirichlet model. Since there is no prior information suggesting bias in the die, we can stipulate a prior that is exchangeable with $n_0 > 0$ "pseudo data-points" in the prior. The model is:

$$\begin{align} \mathbf{n}|\mathbf{p} &\sim \text{Multinomial}(n, \mathbf{p}), \\[12pt] \mathbf{p} &\sim \text{Dirichlet} \Big( \frac{n_0}{k} \cdot \mathbf{1} \Big). \\[6pt] \end{align}$$

Given an observed vector of counts $\mathbf{n}$ we then have the posterior distribution:

$$\mathbf{p}|\mathbf{n} \sim \text{Dirichlet} \Big( \frac{n_0}{k} \cdot \mathbf{1} + \mathbf{n} \Big).$$

For a reasonable choice of the prior strength $n_0$, and a substantial sample size, this posterior ought to give you a reasonable inference for the probability vector $\mathbf{p}$. The Dirichlet distribution is programming in R in various packages (see e.g., here), so it is simple to obtain the density values for any input and generate the posterior density as a function.


Example: Suppose we have a six-sided die with true underlying probabilities:

$$\mathbf{p} = (0.12, \ 0.11, \ 0.13, \ 0.12, \ 0.15, \ 0.37).$$

Here we will conduct a simulation where we roll this die $n=400$ times and derive the resulting posterior distribution, which we program as the function posterior. (Here we use a prior with $n_0 = 10$ pseudo-data points.)

#Generate simulated data
set.seed(1)
n      <- 400
PROB   <- c(0.12, 0.11, 0.13, 0.12, 0.15, 0.37)
k      <- length(PROB)
DATA   <- sample.int(k, size = n, replace = TRUE, prob = PROB)
COUNTS <- table(DATA)

#Determine posterior parameter
n0     <- 10
PAR    <- n0*rep(1, k)/k + COUNTS

#Determine the posterior distribution
posterior <- function(x) { DirichletReg::ddirichlet(x, alpha = PAR) }
$\endgroup$
0
$\begingroup$

Let's look at a formal test of the fairness of a die:

Suppose that the die is biased so that probabilities of faces $1$ through $6$ are $(1,2,2,2,2,3)/12.$ If you roll the die $n = 600$ times, do you have a good chance of detecting the bias?

set.seed(121)
faces = sample(1:6, 600, rep=T, p=c(1,2,2,2,2,3))
freq = tabulate(faces);  freq
[1]  33 111  87 109 109 151

Here is a histogram of the results of 600 rolls of this die.

enter image description here

We do have fewer $1$s and more $6$s than we would expect from a fair die. Does a chi-squared 'goodness-of-fit' test detect the unevenness of the results.

We are testing the null hypothesis that all faces are equally likely (probability $1/6$ each) against the alternative that they are not. The expected counts for each face are $E_i = 100.$ The chi-squared test statistic is $Q = \sum_{i=1}^6 \frac{(X_i - E_i)^2}{E_i},$ where $X_i$ are the observed frequencies. [This is called a chi-squared 'goodness-of-fit' statistic. Notice that large values of $Q$ indicate a bad fit to fairness.]

If the null hypothesis is true, then $Q \stackrel{aprx}{\sim}\mathsf{Chisq}(\nu=6-1=5).$ We would reject the null hypothesis at the 5% level if $Q \ge 11.07.$

qchisq(.95, 5)
[1] 11.0705

In R, the test is shown below: $Q = 75.42 > 11.07$ (X-squared in output). The P-value is the tiny probability of getting a more extreme result if the die were fair.

chisq.test(freq)

        Chi-squared test for given probabilities

data:  freq
X-squared = 75.42, df = 5, p-value = 7.603e-15

This one simulated rolling of the die 600 times is not just a one-time random accident. The following simulation uses 100,000 such tests to show that we will almost always detect such an unfair die with data from $n=600$ rolls.

set.seed(2021)
prob = c(1,2,2,2,2,3);  n = 600
pv = replicate(10^5, 
      chisq.test(tabulate(sample(1:6,n,rep=T,p=prob)) )$p.val)
mean(pv <= .05)
[1] 0.99999

However, $n = 60$ rolls would not be sufficient for vetting such a die. With only 60 rolls of the die, we would detect unfairness only about 35% of the time. [But $n = 300$ roll would suffice for most purposes: 98% rejection. (Not shown.)]

set.seed(121)
prob = c(1,2,2,2,2,3);  n = 60
pv = replicate(10^5, 
      chisq.test(tabulate(sample(1:6,n,rep=T,p=prob)) )$p.val)
mean(pv <= .05)
[1] 0.34642
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.