# Maximum likelihood estimator of joint distribution given only marginal counts

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?

• The margins don't really contain information* about the joint distribution (indeed this is the point of copulas). $\:$ * or at least hardly any - obviously the margins do contain at least some information, since the interior counts can't exceed the margins they occur in. Do you have a specific joint distribution in mind? Why did you use the maximum-entropy tag? Are you after a maximum-entropy solution? – Glen_b Nov 29 '14 at 9:15
• I'm not very familiar with copulas. Do they hold for the categorical case as well? What would that mean - that every joint distribution with the same margins would have the same likelihood? (I tagged maximum-entropy because I thought it might be relevant.) – R S Nov 29 '14 at 10:42
• We don't even have a specified distributional model yet, so we aren't really in a position to compute $P(x|\theta)$. There are numerous possibilities here. Copulas exist for the ordered categorical case (if not unique), but my aim in raising it was to give a motivation for why marginals weren't very informative in general. In respect of the categorical-count case, Fisher treated the margins as uninformative about the joint, whence the Fisher-Irwin exact test. If you want maximum entropy, you probably can get a maximum entropy solution, but I don't know that it will be very informative about... – Glen_b Nov 29 '14 at 21:48
• (ctd) ... structure. In either ME or ML cases, I think you'll first need some kind of model, whether it be bivariate multinomial, bivariate hypergeometric, or something with more structure. See this question, where the author puts a reference into an answer. That may be of help. – Glen_b Nov 29 '14 at 21:51
• I meant a general bivariate multinomial distribution. The question speaks about the case where the sums of the distribution are given and we see samples from the joint distribution. Here we have the sums of the sample. I think the problem is well defined in the ML case (the solution may not be unique but I don't know). – R S Nov 30 '14 at 22:29

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.

• Is it not possible to obtain an analytic solution? – Ben Kuhn Dec 11 '14 at 17:29
• Thanks! The paper seems relevant, although it seems to be from the Bayesian perspective. What about the specific case where $\theta$ is in fact the distribution itself, namely $\theta=\{\theta_{x,y}\}$, for all $(x,y)$ pairs? Would suspect there would be an analytic solution in this case? – R S Dec 12 '14 at 16:10
• I would not suspect there to be an analytic solution. I added an example to illustrate this. – Tom Minka Dec 12 '14 at 16:55
• Thanks. Perhaps it is true asymptotically? Then, conditioning on margin totals is the same as conditioning on margin distributions (after normalizing), and the log-likelihood for each unobserved integer table is proportional to the its entropy. Maybe something with AEP then? – R S Dec 13 '14 at 11:07

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.

• The likelihood in this can would be multinomial. Why is that insufficient? – R S Dec 6 '14 at 17:19
• As I understand it, the likelihood is a function of the parameters given the data. Here, you do not have values for each cell, only the marginals, hence you do not have a single function of the parameters that you can compute, let alone maximize. There are in general many cell configurations compatible with the margins, and each would give a different likelihood. – F. Tusell Dec 6 '14 at 17:28
• Yeah, but that's OK. The parameters are $p$, the data are the marginals. I can still compute the probability of the marginals given $p$ - it's the sum over all the probabilities of cell configurations that give the marginals. That's a single function I can maximize. – R S Dec 6 '14 at 17:50

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.$$

• For the first example: What is given is the marginal counts, not the marginal probabilities. In the case you've described, the probability of $S_1=S_2=T_1=T_2=10$ for the left $p$ is the probability of $[[10,0],[0,10]]$ which is $2^{-20}$. For the right $p$, it is $\sum_{0\le a \le 10}{Pr[[a,10-a],[10-a,a]]}$, which is $10\cdot 4^{-20}$. Even if there is no unique solution, it doesn't mean we can't point to some solution. Maximum entropy gives a unique solution, but it might not be maximum likelihood. – R S Dec 8 '14 at 22:21
• You've calculated the probabilities incorrectly; for instance, you forgot to include the binomial coefficients. But you're right in that the two matrices give different joint distributions of marginal counts even though they give the same marginal distribution of marginal counts. (Yikes!) I'll think about this more. – Ben Kuhn Dec 8 '14 at 23:05