I think the easiest option is to move the correlation out of the output distributions and into the structure of the model. For example, this paper (implemented in PyMC here) uses latent structure to model conditionally independent Poisson variables, which end up being correlated due to the hierarchical model structure .
I think the other alternative would be to work with copula models. @jonsedar has been doing a lot of work on these lately, I think he might be able to chime in with some general remarks on the subject?