Dear Stan users,

I’ve been considering the so called matrix balancing problem - which consist in estimating a matrix, based on known row and column totals (x and y, respectively) only. Rather than estimating matrix elements directly, I estimated row shares (Lambda), for which there is also some prior information available (formulated in terms of Dirichlet priors). The following code illustrates the 2x2 case, for simplicity.

```
data {
vector[2] y; // column totals
vector[2] x; // row totals
vector[2] alpha1;
vector[2] alpha2;
real s;
}
parameters {
simplex[2] Lambda1; // row shares of the unknown matrix (row 1)
simplex[2] Lambda2; // row shares of the unknown matrix (row 2)
}
model {
Lambda1 ~ dirichlet(alpha1);
Lambda2 ~ dirichlet(alpha2);
y[1] ~ normal(Lambda1[1]*x[1] + Lambda2[1]*x[2], s);
y[2] ~ normal(Lambda1[2]*x[1] + Lambda2[2]*x[2], s);
}
```

That formulation (seemingly) worked fine in practice. However, it includes what I consider a workaround, i.e. normal errors, with arbitrarily chosen standard deviations. These standard deviations were set to very small values, implying that estimated column totals equal their target values (y) closely, though still only approximately. But in fact these errors are not a part of the original estimation problem.

My question is - how to program the above problem in Stan avoiding the addition of errors?

That is, instead of assuming “y[1] ~ normal(Lambda1[1]*x[1] + Lambda2[1]*x[2], s)”, I would rather like to impose a strict identity, such as “y[1] = Lambda1[1]*x[1] + Lambda2[1]*x[2]” etc.

Note that reverting to identity does not rule out the need for estimation - the unknown matrix is still underdetermined. Rather it implies using a peculiar likelihood function, which takes the value of 1 for any combination of Lambdas satisfying (in matrix notation) y=Lambda’*x, and 0 otherwise.

I’ll be grateful for your suggestions

Thanks

Jakub

[edit: escaped Stan code]