# Likelihood of two dice. One fair the other unfair

first of all sorry about the title. I couldn't think anything better.

Let me describe the problem first and ask later. Imagine I have 1500 observations of the results of two dice being played (I don't know how many times each die was played). E.g.:

   Sample1
1  187
2  168
3  164
4  320*
5  187
6  174


Seeing that there is something odd with these dice (Number 4 seems to be over-represented), imagine I'm able physically to measure the result bias from each die, E.g.

   Dice 1(fair)  Dice 2(unfair)
[1] 0.1666667    0.04
[2] 0.1666667    0.04
[3] 0.1666667    0.04
[4] 0.1666667    0.80*
[5] 0.1666667    0.04
[6] 0.1666667    0.04


My goal is to measure the proportion of the 1500 observations each die contributed to. One way of doing this is by adding a parameter P (proportion) and perform a non-negative least square regression. However this is limited since it does not consider the observed results a distribution (right?). Another solution may be to use Maximum likelihood estimation. And here is the question. What is the likelihood function that would model this problem? p(O|b,P). Where b is the bias signature and P are the proportion of throws from each die.

The idea here is to have a generic solution. Let's say I have more dice, 10 for example, each on with a specific bias (signature). Is it possible to recover how many times each die was played?

One idea, trying to use likelihood. Dice throws can be modelled by multinomial distributions. For the fair dice, let $\pi=1/6$ be the common probability, for the unfair dice, let the eye probabilities (assumed known) be $\phi_i, i=1,\dotsc,6$ (for instance $\phi_1=0.04$). Let $N=n_1+n_2$, $N=1500$ known, $n_1$ the number of throws with fair dice, $n_2=N-n_1$ number of throws with unfair dice, both unknown. Then the data is $N_i=\text{total number of throws showing$i$eyes}$, which is observed, with $N=N_1+N_2+N_3+N_4+N_5+N_6$. To be able to write down the likelihood, we need the unobserved (potential, or missing) data $N_{1i}, N_{2i}$ which sums to $N_i$. Note that the likelihood here is similar to the one in Estimating parameters for a binomial so some of the points noted in my answer there is relevant here too.
Using the multinomial distribution the likelihood function is $$L=\binom{n_1}{N_{11}\dots N_{16}} \prod_{i=1}^6 \pi^{N_{1i}} \times \binom{n_2}{N_{21}\dots N_{26}} \prod_{i=1}^6 \phi_i^{N_{2i}}$$ so the log likelihood can be written $$\ell = \log(N-n_2)! - \sum_i \log(N_i-N_{2i})! + \log n_2! -\sum_i \log N_{2i}! + N\log \pi +\sum_i N_{2i} \log(\phi_i/\pi)$$ which should be maximized in the unknown parameters (which are integers!) $N_{21}, \dotsc, N_{26}$ (and their sum, which is $n_2$).
To maximize this expression is probably difficult, especially when taking into account the integer restrictions. Another approach is solving the relaxed problem, forgetting that the $N_{2i}$'s are integers. Then we can start by first rewriting the likelihood using that $n!=\Gamma(n+1)$ and https://en.wikipedia.org/wiki/Polygamma_function the polygamma function $\psi(z)$ the first derivative of $\log \Gamma(z)$. $$\ell =\log \Gamma(N-\sum_i N_{2i} +1) - \sum_i \log \Gamma (N_i-N_{2i}+1) + \log\Gamma(\sum_i N_{2i}+1) - \sum_i \log \Gamma(N_{2i}+1) +\sum_i N_{2i} \log(\phi_i/\pi)$$ Maximizing that symbolically doesn't seem easy, so probably numerical optimization is better. I leave that for now!