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?
 A: 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!
