Estimate a cross table from a series of row- and columntotals

Consider the following problem:

Known:
 * \( G^r_g \): number of objects of g-type \(g\) in zone \(r\), aka row-totals.
 * \( H^r_h \): number of objects of h-type \(h\) in zone \(r\), aka column-totals.

Requested:
 * \(Q^r_{gh}\), an estimated number of objects with g-type \(g\) and h-type \(h\) in zone \(r\).
 * estimated number of objects per h-type as a function (preferably a linear combination) of object numbers per h-type (for estimating a h-type distribution given a g-type distribution of new objects).

We assume that:
 * For each zone \(r\), the \(F\) and \(G\) count all objects, thus: \( \forall r: \sum\limits_g G^r_g = \sum\limits_h H^r_h \)
 * \( P^i_h \) is equal for all \(i\) in the same \(r\) and with the same \(g\), thus depends only on \( P(h|rg) \).
 * \(Q^r_{gh}\) has a \(G^r_g\) repeated categorical distribution per row \(g\) and zone \(r\); thus \( E[ Q^r_{gh} ] = P(h|rg) \cdot G^r_g \)
 * \( Q^r_{gh} := f_g \cdot f_h \cdot P_{gh} \) such that \( \sum\limits_{h} Q^r_{gh} = G^r_g \) and \( \sum\limits_{g} Q^r_{gh} = H^r_h \), to be determined by Iterative proportional fitting

Thus:
 * \( \sum\limits_{h} P(h|rg) = 1 \), \( f_g = G^r_g / \sum\limits_{h} f_h \cdot P_{gh} \), and \( f_h = H^r_h / \sum\limits_{g} f_g \cdot P_{gh} \).
 * \( P_{gh} \) is to be determined by regression of \( H^r_h \) by \( G^r_g \), thus, written in matrix notation: \( H = G \times P + \epsilon \) with \(\epsilon\) a \(r \times h \) matrix of independent stochasts with zero expectation.
 * it follows that: \( P := (G^T \times G)^{-1} \times (G^T \times H) \)

ToDo:
 * Consider Alternative: \( Q_{bjh} \) is determined by discrete allocation with the given constraints and suitability \( P_{gh} \)
 * Consider Alternative: \( P := (G^T \times H) \times (H^T \times H)^{-1}\) which follows from regression of \( G^r_g \) by \( H^r_h \).
 * Consider Alternative: \( P := (G^T \times H) \)
 * Consider Effect of possible heteroscedasticity of \(\epsilon\) with \( H \); i.e. assume \( var(\epsilon) \sim H\), as each element of \(H\) is assumed to be the sum over \(j\) of the results of \( G^r_g \) trials of categorical distributions with conditional probabilities \( P_{g|h} \).