Suppose I have a random vector X = (X_1, X_2, ..., X_n) for which X_i ~ Po(\lambda_i), but the components are not independent, so that the joint distribution is not the product of the marginal distributions, and the covariance matrix of X is not diagonal. It seems the joint distribution does not have an explicit form. I would like to obtain some properties from this distribution, as mean, mode, and variance. How can I sample from it? Is there any good approximation to it other than the multivariate normal?