Matrix multinomial distribution

I am trying to fit a model to explain counts of occurrences which occur within two category classes: blood type and (a kind of) HLA type. We have data on the number of individuals in each bin of each category class. So for example, we may have n_{i,j} individuals with blood type i and HLA type j, and we know we have say N_i individuals of blood type i and Z_j individuals of HLA type j in the population.

I was wondering, is there a matrix version of the multinomial distribution? Has anyone worked with such a thing in the past? What I’m hoping is that such a distribution has (the usual) multinomial distributions as its marginals.

Note: I’ve found this paper which discusses the multivariate Dirichlet-multinomial distribution but haven’t yet found anything for the matrix multinomial.

Consider utilizing some basic probability theory. For example, if the categories are independent then P[BT = A \, \text{AND}\, HLA = B] = P[BT = A] * P[HLA = B]. If there are correlations then you’d have to take the overlaps into account.

If you can compute the latent probabilities of each joint event then you can just flatten everything and model all of the joint events with a standard multinomial.

Ok not sure I made myself clear. Forgive me if this seems trivial to some but it’s not to me.

The problem I can’t figure out is how to count all the possible combinations of possible counts such that the marginal counts are fixed. This seems pretty similar to a Latin Square problem (where you have nxn matrices and each row / column contains 1,2,…,n), whose number of combinations are a famously unsolved problem.

Happy to be shown to be wrong here though!