Where is the bomb: How to estimate the probability, given row and column totals?

Question

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:

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!

Answer 1

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:
https://faculty.math.illinois.edu/~mlavrov/seminar/2018-erdos.pdf
https://arxiv.org/pdf/1701.07101.pdf
https://www.tandfonline.com/doi/abs/10.1198/016214504000001303

Answer 2

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.

Example

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.

Independent

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:

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:

Negative Correlation

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

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

Answer 3

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.