Simulate I x J contingency tables using correlation coefficient using bivariate normal distribution

Step 1. generate random values from a multivariate normal distribution
Step 2. convert the random values into probabilities
Step 3. convert the probabilities into values from the target distribution

This requires a cautionary note.
The method above is simple in the univariate case (where people would
usually start with a uniform distribution), however, generalizing it
to the multivariate case, is not as simple as Thierry as suggested.

That doesn't necessarily mean that his solution is wrong.
(Because there's more than one way of defining a bivariate Poisson
distribution).