Sure there is.
Let's assume that we could see $k$ different types of dice with $Y=\{y_1, \dots, y_k\}$ pips. The "classical" D&D dice would yield $Y=\{2, 4, 6, 8, 10, 12, 20, 100\}$, but anything else could be made to work as below. Let $|Y|$ denote the number of different possible dice.
Assume that the minimum observed roll is $m$ and the maximum is $M$. This gives us bounds on the possible number of dice involved: we must have
$$m\leq x\leq M.$$
Now, simulate a large number of rolls from each remaining candidate xdy (by the formula above, there are $(M-m+1)\times |Y|$ different candidates). This large number can be larger than the number of original observations. In each case, create a histogram of total outcomes.
Finally, compare these density histograms to the histogram of the original observations. (We work with the densities, not the frequencies, so we can use different numbers of simulations from the original number of observations.) The earth mover's distance gives you a notion of distance between histograms. The simulation-based histogram with the lowest earth mover's distance to your original observation histogram is most consistent with your observations.
Here is some R code. We simulate 1000 rolls of 3d6 and use the standard D&D dice as candidates $Y$.
set.seed(1)
obs <- rowSums(matrix(ceiling(runif(3000,0,6)),nrow=1000))
candidates.y <- c(4,6,8,10,12,20,100) # standard D&D dice
breaks <- c(seq(.5,max(obs)+.5,by=1),max(obs)*max(candidates.y))
hist.obs <- hist(obs,breaks=breaks,freq=FALSE)
candidates.x <- 1:max(obs)
n.sims <- 1000
library(emdist)
emd.dist <- matrix(NA,nrow=length(candidates.x),ncol=length(candidates.y),
dimnames=list(candidates.x,candidates.y))
for ( ii in seq_along(candidates.x) ) {
for ( jj in seq_along(candidates.y) ) {
set.seed(1) # for replicability
obs.sim <- rowSums(matrix(ceiling(
runif(candidates.x[ii]*n.sims,0,candidates.y[jj])),
nrow=n.sims))
hist.sim <- hist(obs.sim,breaks=breaks,plot=FALSE)
emd.dist[ii,jj] <- emd2d(matrix(hist.obs$density,nrow=1),
matrix(hist.sim$density,nrow=1))
}
}
emd.dist
Output:
4 6 8 10 12 20 100
1 7.9029999 6.8909998 5.8920002 4.8950000 3.88700008 1.9228243 0.39887184
2 5.4299998 3.4460001 1.4499999 0.8147613 1.08018708 0.8099990 0.04417419
3 2.9579999 0.0000000 1.8063331 1.9216927 1.61559486 0.2684011 0.35928789
4 0.6180000 2.4049284 2.2750988 1.3528987 0.71173453 0.2177625 1.00000000
5 1.9768735 2.8317218 1.7545075 0.5033562 0.17459151 0.1226299 1.00000000
6 3.4169219 1.9913402 0.4674550 0.1918677 0.05264558 1.0000000 1.00000000
7 3.3924465 0.9393466 0.2968451 0.1572183 0.21881166 1.0000000 1.00000000
8 2.6134274 0.3687483 0.1226299 1.0000000 1.00000000 1.0000000 1.00000000
9 1.5483801 0.3894975 0.3592879 1.0000000 1.00000000 1.0000000 1.00000000
10 0.3270817 0.2188117 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
11 0.1052912 0.3592879 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
12 0.3592879 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
13 1.0000000 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
14 1.0000000 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
15 1.0000000 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
16 1.0000000 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
17 1.0000000 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
18 1.0000000 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000 1.00000000
And look, the 3d6 in fact does have the simulated histogram with the lowest distance from your observations' histogram. Hurray!
A couple of comments:
- This works better if you have "many" observations. With few observations, you get a lot of noise, and many different simulated histograms are consistent with your noisy data.
- Of course, the larger your candidate set $Y$, the lower your chance will be of getting the true values.
- If the true $y$ is not in $Y$, you will at least get something "close to it".
- The formally-statistically-correct treatment would be a Bayesian one, putting equal prior probabilities on all possible $x$ and $y\in Y$, then calculating or simulating draws and deriving posterior probabilities on pairs $(x,y)$. Yes, this can be done rigorously. Anyone want to have a go?