This question is inspired by a mini-game from Pokemon Soulsilver:

Imagine there are 15 bombs hidden on this 5x6 area (EDIT: maximum 1 bomb/cell):


Now, how would you estimate the probability to find a bomb on a specific field, given the row/column totals?

If you look at column 5 (total bombs = 5), then you might think: Within this column the chance to find a bomb in row 2 is double the chance to find one in row 1.

This (wrong) assumption of direct proportionality, which basically can be described as drawing standard independence-test operations (like in Chi-Square) into the wrong context, would lead to the following estimations:


As you can see, direct proportionality leads to probability estimates over 100%, and even before that, would be wrong.

So I performed a computational simulation of all possible permutations which led to 276 unique possibilities of placing 15 bombs. (given row and column totals)

Here is the average over the 276 solutions: Computational solution

This is the correct solution, but due to the exponential computational work, I would like to find an estimation method.

My question is now: Is there an established statistical method to estimating this? I was wondering if this was a known problem, how it is called and if there are papers/websites you could recommend!

  • 1
    $\begingroup$ Fast and easy approach: For higher number of rows & columns, you could conduct a Monte Carlo simulation, where you would check the random subsample of the possible configurations that is lower then then total number of possibilities. It'd give you an approximate solution. $\endgroup$
    – Tim
    Sep 23, 2019 at 14:24
  • 1
    $\begingroup$ I dont understand your computational solution. What are the numbers in cells? They certainly don't add up to 100%, it's not PMF. They also dont look like CDF, the right/bottom cell is not 100% $\endgroup$
    – Aksakal
    Sep 23, 2019 at 18:12
  • 2
    $\begingroup$ @Aksakal These are the marginal probabilities that any given cell contains a bomb. The numbers add to 15, the number of total bombs on the board. $\endgroup$
    – Danica
    Sep 23, 2019 at 19:21
  • 2
    $\begingroup$ If you're assuming the two margins are independent its relatively straightforward to sample from the distribution of tables conditional on the margins (via Patefield's algorithm). This is implemented in the standard distribution of R in r2dtable (and also used by chisq.test and fisher.test in some circumstances). $\endgroup$
    – Glen_b
    Sep 24, 2019 at 4:08
  • 2
    $\begingroup$ @Glen_b But in the Patefield algorithm the number of events per cell is not limited to one. $\endgroup$ Sep 24, 2019 at 14:30

3 Answers 3


The solution space (valid bomb configurations) can be viewed as the set of bipartite graphs with given degree sequence. (The grid is the biadjacency matrix.) Generating a uniform distribution on that space can be approached using Markov Chain Monte Carlo (MCMC) methods: every solution can be obtained from any other using a sequence of "switches," which in your puzzle formulation look like:

$$ \begin{pmatrix} x & - \\ - & x \end{pmatrix} \to \begin{pmatrix} - & x \\ x & - \end{pmatrix} $$

It's been proven that this has a fast mixing property. So, starting with any valid configuration and setting a MCMC running for a while, you should end up with an approximation of the uniform distribution on solutions, which you can average pointwise for the probabilities you're looking for.

I'm only vaguely familiar with these approaches and their computational aspects, but at least this way you avoid enumerating any of the non-solutions.

A start to the literature on the topic:

  • $\begingroup$ That's an amazing idea! I think I get it! I mix through any known solution for a defined amount of iterations (that I expect to find in the papers) and afterwards average over the unique solutions, hoping most of them are found. Thanks so much! $\endgroup$
    – KaPy3141
    Sep 25, 2019 at 9:22
  • 2
    $\begingroup$ MCMC is exactly the way to go and I also found this: arxiv.org/pdf/1904.03836.pdf $\endgroup$
    – KaPy3141
    Sep 25, 2019 at 10:33
  • $\begingroup$ @KaPy3141 For the above row and column sums, my implementation of the rectangle loop algorithm (in the arxiv preprint) only visits 276 unique states even if I run the algorithm for as many as $10^6$ iterations. $\endgroup$ Sep 25, 2019 at 13:52
  • $\begingroup$ Which suggests that enumeration as suggested by @Aksakal may be more efficient. $\endgroup$ Sep 25, 2019 at 13:52
  • $\begingroup$ @JarleTufto, but OP says there are only 276 unique (valid) states; you've found them all! $\endgroup$ Sep 25, 2019 at 13:54

There's no unique solution

I don't think that true discrete probability distribution can be recovered, unless you make some additional assumptions. Your situation is basically a problem of recovering the joint distribution from marginals. It is sometimes solved by using copulas in the industry, for example financial risk management, but usually for continuous distributions.

Presence, Independent, AS 205

In the presence problem no more than one bomb is allowed in a cell. Again, for the special case of independence, there is relatively efficient computational solution.

If you know FORTRAN, you can use this code that implements AS 205 Algorithm: Ian Saunders, Algorithm AS 205: Enumeration of R x C Tables with Repeated Row Totals, Applied Statistics, Volume 33, Number 3, 1984, pages 340-352. It's related to Panefield's algo that @Glen_B referred to.

This algo enumerates all presence tables, i.e. goes through all possible tables where only one bomb is in a field. It also calculates the multiplicity, i.e. multiple tables that look the same, and calculates some probabilities (not those you're interested in). With this algorithm you may be able to run the complete enumeration faster than you did before.

Presence, not independent

The AS 205 algorithm can be applied to a case where the rows and columns are not independent. In this case you'd have to apply different weights to each table generated by the enumeration logic. The weight will depend on the process of placement of bombs.

Counts, independence

The count problem allows more than one bomb placed in a cell, of course. The special case of independent rows and columns of count problem is easy: $P_i^j=P_i\times P^j$ where $P_i$ and $P^j$ are marginals of rows and columns. For instance, row $P_6=3/15=0.2$ and column $P^3=3/15=0.2$, hence the probability that a bomb is in row 6 and column 3 is $P_6^3=0.04$. You actually produced this distribution in your first table.

Counts, Not independent, Discrete Copulas

In order to solve the counts problem where rows and columns are not independent, we could apply discrete copulas. They have issues: they're not unique. It doesn't make them useless though. So, I'd try applying discrete copulas. You can find a good overview of them in Genest, C. and J. Nešlehová (2007). A primer on copulas for count data. Astin Bull. 37(2), 475–515.

Copulas can be especially useful, as they usually allow to explicitly induce dependence, or to estimate it from data when the data is available. I mean the dependence of row and columns when placing bombs. For instance, it could be the case when if the bomb is one the first row, then it is more likely that it will be one the first column too.


Let's apply Kimeldorf and Sampson copula to your data, assuming again that more than one bomb can be placed in a cell. The copula for a dependency parameter $\theta$ is defined as: $$C(u,v)=(u^{-\theta}+u^{-\theta}-1)^{-1/\theta}$$ You can think of $\theta$ as an analog of the correlation coefficient.


Let's start with the case of weak dependence, $\theta=0.000001$, where we have the following probabilities (PMF) and marginal PDFs are shown too on the panels on the right and at the bottom:

enter image description here

You can see how in the column 5 the second row probability does have twice higher probability than the first row. This is not wrong contrary to what you seemed to imply in your question. All probability do add up to 100%, of course, as do the marginals on the panels match the frequencies. For instance, the column 5 in the lower panel shows 1/3 which corresponds to stated 5 bombs out of total 15 as expected.

Positive Correlation

For stronger dependency (positive correlation) with $\theta=10$ we have the following:

enter image description here

Negative Correlation

The same for stronger but negative correlation (dependency) $\theta=-0.2$:

enter image description here

You can see that all probabilities add up to 100%, of course. Also, you can see how dependency impacts the PMF's shape. For positive dependency (correlation) you get the highest PMF concentrated on the diagonal, while for the negative dependency it is off-diagonal

  • $\begingroup$ Thanks so much for your answer and your interesting links to copulas! Unfortunately, I have never used copulas, so it will be hard for me to find a solution that enforces only 1 bomb per cell, but I will definitely try once I have a better understanding! $\endgroup$
    – KaPy3141
    Sep 24, 2019 at 11:02
  • $\begingroup$ @KaPy3141, I added reference to the code that you can use to solve the problem. It's in F90, but relatively straightforward to convert to Python with numpy $\endgroup$
    – Aksakal
    Sep 24, 2019 at 20:24
  • $\begingroup$ How are the copula's with some parameter $\theta $ a solution to the problem? How do you determine $\theta $ and how are you gonna know that it is the answer (for instance, a strange effect in your answer is that rows with the same marginal probability give different cell probabilities). The problem seems like a combinatorial problem to me. $\endgroup$ Sep 25, 2019 at 16:31
  • $\begingroup$ You'd have to fit the parameters to the process. The problem is pure combinatorial if the generating process is consistent with it. $\endgroup$
    – Aksakal
    Sep 25, 2019 at 16:44

Your question does not make this clear, but I'm going to assume that the bombs are initially distributed via simple-random-sampling without replacement over the cells (so a cell cannot contain more than one bomb). The question you have raised is essentially asking for the development of an estimation method for a probability distribution that can be computed exactly (in theory), but which becomes computationally infeasible to compute for large parameter values.

The exact solution exists, but it is computationally intensive

As you point out in your question, it is possible for you to perform a computational search over all possible allocations, to identify the allocations that match the row and column totals. We can proceed formally as follows. Suppose we are dealing with an $n \times m$ grid and we allocate $b$ bombs via simple random sampling without replacement (so each cell cannot contain more than one bomb).

Let $\mathbf{x} = (x_1,...,x_{nm})$ be a vector of indicator variables indicating whether or not a bomb is present in each cell, and let $\mathbf{s} = (r_1, ..., r_n, c_1, ..., c_m)$ denote the corresponding vector of row and column sums. Define the function $S: \mathbf{x} \mapsto \mathbf{s}$, which maps from the allocation vector to the row and column sums.

The goal is to determine the probability of each allocation vector conditional on knowledge of the row and column sums. Under simple-random-sampling we have $\mathbb{P}(\mathbf{x}) \propto 1$, so the conditional probability of interest is:

$$\begin{equation} \begin{aligned} \mathbb{P}(\mathbf{x} | \mathbf{s}) = \frac{\mathbb{P}(\mathbf{x}, \mathbf{s})}{\mathbb{P}(\mathbf{s})} &= \frac{\mathbb{P}(\mathbf{x}) \cdot \mathbb{I}(S(\mathbf{x}) = \mathbf{s})}{\sum_\mathbf{x} \mathbb{P}(\mathbf{x}) \cdot \mathbb{I}(S(\mathbf{x}) = \mathbf{s})} \\[6pt] &= \frac{\mathbb{I}(S(\mathbf{x}) = \mathbf{s})}{\sum_\mathbf{x} \mathbb{I}(S(\mathbf{x}) = \mathbf{s})} \\[6pt] &= \frac{1}{|\mathscr{X}_\mathbf{s}|} \cdot \mathbb{I}(S(\mathbf{x}) = \mathbf{s}) \\[6pt] &= \text{U}(\mathbf{x} | \mathscr{X}_\mathbf{s}), \\[6pt] \end{aligned} \end{equation}$$

where $\mathscr{X}_\mathbf{s} \equiv \{ \mathbf{x} \in \{ 0, 1\}^{nm} | S(\mathbf{x}) = \mathbf{s} \}$ is the set of all allocation vectors compatible with the vector $\mathbf{s}$. This shows that (under simple-random-sampling of the bombs) we have $\mathbf{x} | \mathbf{s} \sim \text{U}(\mathscr{X}_\mathbf{s})$. That is, the conditional distribution of the allocation vector for the bombs is uniform over the set of all allocation vectors compatible with the observed row and column totals. The marginal probability of a bomb in a given cell can then be obtained by marginalising over this joint distribution:

$$\begin{equation} \begin{aligned} \mathbb{P}(x_{ij} = 1 | \mathbf{s}) = \sum_{\mathbf{x}: x_{ij} = 1} \text{U}(\mathbf{x} | \mathscr{X}_\mathbf{s}) = \frac{|\mathscr{X}_{ij} \cap \mathscr{X}_\mathbf{s}|}{|\mathscr{X}_\mathbf{s}|}. \end{aligned} \end{equation}$$

where $\mathscr{X}_{ij} \equiv \{ \mathbf{x} \in \{ 0, 1\}^{nm} | x_{ij} = 1 \}$ is the set of all allocation vectors with a bomb in the cell in the $i$th row and $j$th column. Now, in your particular problem you computed the set $\mathscr{X}_\mathbf{s}$ and found that $|\mathscr{X}_\mathbf{s}| = 276$, so the conditional probability of the allocation vectors is uniform over the set of allocations you computed (assuming you did this correctly). This is ane exact solution to the problem. However, it is computationally intensive to compute the set $\mathscr{X}_\mathbf{s}$ and so the computation of this solution may become infeasible when $n$, $m$ or $b$ become larger.

Searching for good estimation methods

In the case where it is infeasible to compute the set $\mathscr{X}_\mathbf{s}$, you want to be able to estimate the marginal probabilities of a bomb being in any particular cell. I am not aware of any existing research that gives estimation methods for this problem, so this is going to require you to develop some plausible estimators and then test their performance against the exact solution using computer simulations for parameter values that are sufficiently low for this to be feasible.

The naive empirical estimator: The estimator you have proposed and used in your green table is:

$$\hat{\mathbb{P}}(x_{ij} = 1 | \mathbf{s}) = \frac{r_i}{b} \cdot \frac{c_j}{b} \cdot b = \frac{r_i \cdot c_j}{b}.$$

This estimation method treats the rows and columns as independent, and estimates the probability of a bomb in a particular row/column by the relative frequencies in the row and column sums. It is simple to establish that this estimator sums to $b$ over all the cells, as you would want. Unfortunately, it has the major drawback that it can yield an estimated probability above one in some cases. That is a bad property for an estimator.

  • $\begingroup$ Thanks so much for your in-depth answer! Actually, in my green chart, there already are values up to 133%. It's good to know that there is no popular method for this problem and it's acceptable to experiment for oneself! My most accurate estimator is similar to the "green" approach, but instead of allocating the bombs proportional to P(row)/sum(P(rows)) * P(c)/sum(P(cols)), I use an imaginary P(r)/(1-P(r))/sum(rows) and afterwards bring the product back: P(real) = P(imag)/(1+P(imag). This forces P < 1. Now I guess, I just need to computationally enforce the (slightly violated) row/column sums. $\endgroup$
    – KaPy3141
    Sep 24, 2019 at 10:57
  • $\begingroup$ @KaPy3141 you might use the value that a specific bomb is in a cell (which does not have the problem of being above 1) and then describe the problem as a draw of 15 bombs out of that distribution with the condition that each cell has only values 0 or 1 (drawing without replacement). This will provide you with a probability that does not exceed 1. $\endgroup$ Sep 25, 2019 at 16:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.