# Estimating the number of times each of four pairs of dice was thrown

## PROBLEM

Suppose I have four dice: red (R), green (G), blue (B) and yellow (Y). I throw these dice several times. Each throw uses two of the four dice, and only the following pairs are allowed:

• RB
• RY
• GB
• GY

After each throw I note which dies rolled a 'one'. Across the four allowed pairs there are nine possible results:

1. Both R and B rolled a 'one'
2. Both R and Y rolled a 'one'
3. Both G and B rolled a 'one'
4. Both G and Y rolled a 'one'
5. Only R rolled a 'one'
6. Only G rolled a 'one'
7. Only B rolled a 'one'
8. Only Y rolled a 'one'
9. Neither of the dice rolled a 'one'

As I throw the dice, I make a tally of how many times I got each of results 1 through 8. I do not make a tally of how many times I got result 9.

At the end of this process I give you my tallies, so that you know the number of times I got each of the first eight results. You also know the pairs that were allowed and the probability of rolling a 'one' (which is the same for all dice). However, you know neither the total number of throws I made, nor the number of times I got result 9.

Using this information, your problem is to estimate the number of times I threw each of the four allowed pairs.

Some things to note:

• Order does not matter; both dice are thrown together.
• Multiple allowed pairings can have the same result; for example, both RB and RY can have result 5.

## QUESTIONS

• Is the problem solvable?
• How would you solve this problem?
• What are some names for this type of problem that would be useful search terms if I want to look for more information?
-

This reply addresses the first and third bullet points: Is the problem solvable? Can we frame it in a conventional way to allow for searching or application of conventional methods?

To address the latter question it helps to start generalizing the situation while retaining any special features that might be useful for a solution. Let's begin with the raw data. The experiment is a sequence of trials. In each trial (a) a pair of dice is rolled and (b) we record the outcome of each die (one or not-one) and its color. This can be represented by sixteen values: four possibilities for each of the four pairs of dice. We know the probabilities associated with each set of four outcomes: $1/6^2$ for two ones, $(5/6)^2$ for two non-ones, and $(1/6)(5/6)$ for each of the remaining two outcomes.

To summarize the data, suppose the red and blue dice were thrown $n_{rb}$ times, the red and yellow dice $n_{ry}$ times, etc. This means we have observed the outcomes of $n_{rb}$ independent throws of the red and blue dice, etc. The sum of those red-blue throws is therefore the outcome of a multinomial distribution with count parameter $n_{rb}$ and probabilities $(1/6, 5/36, 5/36, 25/36)$. Similarly the sum of the red-yellow throws is an independent outcome of a multinomial distribution with count parameter $n_{ry}$ and the same probabilities; the green-blue throws have count parameter $n_{gb}$, and the green-yellow throws have count parameter $n_{gy}$. These four distributions, in this order, collectively describe a $16$-variate distribution.

A visualization can help, so let's consider an example. Here are some raw data with descriptive headers:

                      Red-blue      Red-yellow    Green-blue    Green-yellow
00 01 10 11   00 01 10 11   00 01 10 11   00 01 10 11
Red Green Blue Yellow k0 k1 k2 k3   k4 k5 k6 k7   k8 k9 kA kB   kC kD kE kF
*          *            1
*                 *                1
*    *                                              1
*           *                                                  1


The variables are named k0 through kF: evidently they are a dummy coding for the $16$ possible outcomes. The outcomes are schematically shown in the second line: "00" means both dice were non-one, "01" means only the second die (as named on the first line) showed a one, etc. Redundantly, stars indicate which two dice were thrown: I have shown them only to illustrate what's going on. Thus, this dataset describes four trials: in the first, red was non-1 and blue was 1; in the second, red and blue were both non-1; in the third, green and blue were both 1; and in the fourth, green was 1 and yellow was non-1.

A sufficient statistic for this experiment would be the sum of all the data rows (taking blanks to be zeros): this counts each of the 16 kinds of outcomes.

We do not observe these raw data, though: they are condensed for us. Specifically,

1. The count of cases where both R and B showed 1 is the sum in the k3 column.
2. The count of cases where both R and Y showed 1 is the sum in the k7 column..
3. The count of cases where both G and B showed 1 is the sum in the kB column.
4. The count of cases where both G and YB showed 1 is the sum in the kF column.
5. The count of cases where only R showed 1 is the sum of the k2 and k6 columns.
6. The count of cases where only G showed 1 is the sum of the kA and kE columns.
7. The count of cases where only B showed 1 is the sum of the k1 and k9 columns.
8. The count of cases where only Y showed 1 is the sum of the k5 and kd columns.
9. The number of cases where neither die showed 1 is the sum of the k0, k4, k8, and kC columns (but this is not revealed to us).

Writing $\mathbb{k}$ for the column matrix of k's (in the order shown), this information is conveniently written as a linear transformation $\mathbb{A k}$ where the matrix $\mathbb{A}$ is

$$\left( \begin{array}{cccccccccccccccc} 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{array} \right)$$

It is now immediate--your computer will happily tell you this--that the last four rows are redundant (the sum of rows 5 and 6 equals the sum of rows 7 and 8). We might just as well drop the last row from $\mathbb{A}$: it will then be of full rank.

Here, then, is the abstract statement of the problem:

Given one realization $\mathbb{k}$ from a family of multivariate distributions parameterized by natural numbers $\mathbb{n} = (n_{rb}, n_{ry}, n_{gb}, n_{gy})$ and given an observation of $\mathbb{A k}$, estimate the parameters (and obtain standard errors or confidence limits on those estimates).

We probably should add that the entries of $\mathbb{k}$ are themselves natural numbers: this is a discrete distribution.

In general, the structure of $\mathbb{A}$ induces dependencies among the entries of $\mathbb{A k}$. (Note that the entries of $\mathbb{k}$ themselves already have some slight dependencies arising from the underlying multinomial distributions). This, along with the discrete distribution of $\mathbb{A k}$ and discreteness of the parameter space, are going to create difficulties in developing estimators.

We can at least get started by taking expectations, because it's easy to write the expectation of $\mathbb{k}$ in terms of the $n_{*}$: $E[k_0] = n_{rb}25/36$, $E[k_1] = n_{rb}5/36$, ..., and $E[k_A] = n_{gy}/36$. The linearity of expectation tells us that $E[\mathbb{Ak}] = \mathbb{A}E[\mathbb{k}]$. Working this out gives us a lot of possible method-of-moments estimators (not just one). (One of them was posted as a reply by the OP.) So yes, the problem is solvable. (The generalized problem could have a unique method-of-moments estimator or it might have none at all when $\mathbb{A}$ does not give enough information to identify all the parameters.)

The important questions left to solve are:

• How well can it be solved? Can we find good (e.g., admissible) estimators?

• Can we obtain good confidence intervals or other expressions of uncertainty in the estimates?

We can continue and compute the variance of these observations using standard rules, using second moments of the multinomial distribution. With this in hand one might be tempted to combine the seven observations (the counts in 1-7) using generalized least squares. Or, one might proceed directly to try a maximum likelihood approach (but this would be quite tricky to compute). When the components of $\mathbb{k}$ are expected to be large, normal approximations to the multinomial distribution will work nicely, for then $\mathbb{A k}$ will also be (approximately) multivariate normal and maximum likelihood estimates of the $n_{*}$ might be well behaved.

That's the (limited) extent of my analysis. I wanted to share it at this point to give something for future answers to build on and to show the complexities and pitfalls involved in this apparently simple situation.

-

Let p be the probability of throwing a 1 on any die. Let ni be the probability of the outcome designated by i. The probability of the occurrences 1,2,3 and 4 are all p$^2$ as the dice are rolled independently. The probability of occurrences 5, 6, 7, and 8 are all 2p(1-p) since for example in case 5 a red must roll a 1 and the blue roll a number other than 1 or red rolls 1 and the yellow die rolls a number other than 1.

Hence E(n1)=E(n2)=E(n3)=E(n4)=Np$^2$

and E(n5)=E(n6)=E(n7)=E(n8)=2Np(1-p).

Estimate N by equating the nis with their expectations to get since

n9=N-(n1+n2+n3+n4+n5+n6+n7+n8)

and E(n1+n2+n3+n4)=4Np$^2$

and E(n5+n6+n7+n8) =4(2Np(1-p))

E(n9)=N-4N(p$^2$ +2p(1-p))=

N(1-4(p$^2$+2p-2p$^2$))=N(1-8p+4p$^2$)

This shows that if N were known I could estimate n9 by the nearest integer to N(1-8p+4p$^2$). Or on the other hand if I knew n9 I could estimate N by the nearest integer to n9/(1-8p+4p$^2$).

But if I did not know N and I also do not know n9 then for any solution I could multiply N and n9 by any integer (say 2 or 5) and get another answer. So without additional information I cannot find a sensible and unique estimate for n9.

-
Very good response, but I'm still wondering if there might be way to estimate N. A conceptual example: suppose p is close to one. Suppose also that I got response 1 one time, and responses 2-8 zero times. Having this result would be very unlikely if N were very large, so we can conclude that N is probably not very large. Is there a way to estimate the N that maximizes the probability of the results n1-n8? – ericstalbot Jun 27 '12 at 18:15
@est I think that if you knew the proportion of times each pair is tossed you could maximize the likelihood of the observed outcomes given N and compute a maximum likelihood estimate for N. But you haven't specified the frequencies with which each of the 4 pairs is tossed and so you can't uniquely define the likelihood. – Michael Chernick Jun 27 '12 at 18:25
$N_9$ is easy to estimate; you know p and the total number of observed rolls N, so, knowing p, you can calculate the probability that a roll is not observed: $(1−p)^2$ (neither of the two dice comes up 1.) From that, you get $N_9 = N(1−p)^2/(1−(1−p)^2)$. – jbowman Jun 27 '12 at 19:18
@jbowman The Op specified initially that neither n9 nor N are known. Otherwise the problem has solutions. Also he is specifying occurrences for specific pairs of dice but does not specify whether or not each pair is rolled the same number of times. He has too many unknowns. – Michael Chernick Jun 27 '12 at 19:27
There aren't too many unknowns. If you ignore the counts 5-8 you have a simple standard problem. You have more information than this. – Douglas Zare Jun 27 '12 at 20:08
Let $p$ be the probability of any die rolling a one. Then the probability of outcomes 1, 2, 3 and 4 are all $p^2$. Let $n_1$, $n_2$, $n_3$ and $n_4$ be number of observations for outcomes 1 through 4. Let $n_{rb}$, $n_{ry}$, $n_{gb}$ and $n_{gy}$ be the values we want to estimate, which are the actual number of times each of the dice pairs was thrown. Then the best estimates are:
$$n_{rb} = \frac{n_1}{p^2}$$ $$n_{ry} = \frac{n_2}{p^2}$$ $$n_{gb} = \frac{n_3}{p^2}$$ $$n_{gy} = \frac{n_4}{p^2}$$