Maximum likelihood estimator of joint distribution given only marginal counts

Question

Let $p_{x,y}$ be a joint distribution of two categorical variables $X,Y$, with $x,y\in\{1,\ldots,K\}$. Say $n$ samples were drawn from this distribution, but we are only given the marginal counts, namely for $j=1,\ldots,K$:

$$ S_j = \sum_{i=1}^{n}{\delta(X_i=l)}, T_j = \sum_{i=1}^{n}{\delta(Y_i=j)}, $$

What is the maximum likelihood estimator for $p_{x,y}$, given $S_j,T_j$? Is this known? Computationally feasible? Are there any other reasonable approaches to this problem other than ML?

Answer 1

This kind of problem was studied in the paper "Data Augmentation in Multi-way Contingency Tables With Fixed Marginal Totals" by Dobra et al (2006). Let $\theta$ denote the parameters of the model, let $\mathbf{n}$ denote the unobserved integer table of counts for each $(x,y)$ pair, and let $C(S,T)$ be the set of integer tables whose marginal counts equal $(S,T)$. Then the probability of observing the marginal counts $(S,T)$ is: $$ p(S,T | \theta) = \sum_{\mathbf{n} \in C(S,T)} p(\mathbf{n} | \theta) $$ where $p(\mathbf{n} | \theta)$ is the multinomial sampling distribution. This defines the likelihood function for ML, but direct evaluation is infeasible except for small problems. The approach they recommend is MCMC, where you alternately update $\mathbf{n}$ and $\theta$ by sampling from a proposal distribution and accepting the change according to the Metropolis-Hastings acceptance ratio. This could be adapted to find an approximate maximum over $\theta$ using Monte Carlo EM.

A different approach would use variational methods to approximate the sum over $\mathbf{n}$. The marginal constraints can be encoded as a factor graph and inference over $\theta$ could be carried out using Expectation Propagation.

To see why this problem is difficult and does not admit a trivial solution, consider the case $S=(1,2), T=(2,1)$. Taking $S$ as the row sums and $T$ as the column sums, there are two possible tables of counts: $$ \begin{bmatrix} 0 & 1 \\ 2 & 0 \end{bmatrix} \qquad \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix} $$ Therefore the likelihood function is $$ p(S,T|\theta) = 3 p_{12} p_{21}^2 + 6 p_{11} p_{21} p_{22} $$ The MLE for this problem is $$ \hat{p}_{x,y} = \begin{bmatrix} 0 & 1/3 \\ 2/3 & 0 \end{bmatrix} $$ which corresponds to assuming the table on the left. By contrast, the estimate that you would get by assuming independence is $$ q_{x,y} = \begin{bmatrix} 1/3 \\ 2/3 \end{bmatrix} \begin{bmatrix} 2/3 & 1/3 \end{bmatrix} = \begin{bmatrix} 2/9 & 1/9 \\ 4/9 & 2/9 \end{bmatrix} $$ which has a smaller likelihood value.

Answer 2

As has been pointed by @Glen_b, this is insufficiently specified. I do not think you can use maximum likelihood unless you can fully specify the likelihood.

If you were willing to assume independence, then the problem is quite simple (incidentally, I think the solution would be the maximum entropy solution that has been suggested). If you are not willing nor able to impose additional structure in your problem and you still want some kind of approximation to the values of the cells, may be you could use the Fréchet–Hoeffding copula bounds. Without additional assumptions, I do not think you can go any further.

Answer 3

Edit: This answer is based on an incorrect assumption that likelihood of the marginal counts given $p_{x,y}$ is only a function of the marginal probabilities $p_x = \sum_y p_{x,y}$ and $p_y = \sum_x p_{x,y}$. I'm still thinking about it.

Wrong stuff follows:

As mentioned in a comment, the problem with finding "the" maximum-likelihood estimator for $p_{x, y}$ is that it's not unique. For instance, consider the case with binary $X, Y$ and marginals $S_1 = S_2 = T_1 = T_2 = 10$. The two estimators

$$p = \left(\begin{array}{cc} \frac12 & 0 \\ 0 & \frac12\end{array}\right), \qquad p = \left(\begin{array}{cc} \frac14 & \frac14 \\ \frac14 & \frac14\end{array}\right)$$

have the same marginal probabilities $p_x$ and $p_y$ in all cases, and hence have equal likelihoods (both of which maximize the likelihood function, as you can verify).

---

Indeed, no matter what the marginals are (as long as two of them are nonzero in each dimension), the maximum likelihood solution is not unique. I'll prove this for the binary case. Let $p = \left(\begin{array}{cc}a & b \\ c & d\end{array}\right)$ be a maximum-likelihood solution. Without loss of generality suppose $0 < a \le d$. Then $p = \left(\begin{array}{cc}0 & b + a \\ c + a & d - a\end{array}\right)$ has the same marginals and is thus also a maximum-likelihood solution.

---

If you want to additionally apply a maximum-entropy constraint, then you do get a unique solution, which as F. Tussell stated is the solution in which $X, Y$ are independent. You can see this as follows:

The entropy of the distribution is $H(p) = -\sum_{x,y} p_{x,y} \log p_{x,y}$; maximizing subject to $\sum_x p_{x,y} = p_y$ and $\sum_{y} p_{x,y} = p_x$ (equivalently, $\vec g(p) = 0$ where $g_x(p) = \sum_y p_{x,y} - p_x$ and $g_y(p) = \sum_x p_{x,y} - p_y$) using Lagrange multipliers gives the equation:

$$\nabla H(p) = \sum_{ k \in X \cup Y} \lambda_k \nabla g_k(p) $$

All the gradients of each $g_k$ are 1, so coordinate-wise this works out to

$$1 - \log p_{x,y} = \lambda_x + \lambda_y \implies p_{x,y} = e^{1-\lambda_x-\lambda_y}$$

plus the original constraints $\sum_x p_{x,y} = p_y$ and $\sum_{y} p_{x,y} = p_x$. You can verify that this is satisfied when $e^{1/2 - \lambda_x} = p_x$ and $e^{1/2 - \lambda_y} = p_y$, giving $$p_{x,y} = p_xp_y.$$