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,
- The count of cases where both R and B showed 1 is the sum in the
k3
column.
- The count of cases where both R and Y showed 1 is the sum in the
k7
column..
- The count of cases where both G and B showed 1 is the sum in the
kB
column.
- The count of cases where both G and YB showed 1 is the sum in the
kF
column.
- The count of cases where only R showed 1 is the sum of the
k2
and k6
columns.
- The count of cases where only G showed 1 is the sum of the
kA
and kE
columns.
- The count of cases where only B showed 1 is the sum of the
k1
and k9
columns.
- The count of cases where only Y showed 1 is the sum of the
k5
and kd
columns.
- 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.