On Thu, 26 Sep 2013, Umberto Triacca wrote: The problem is not trivial, in fact. As long as the (joint) density you're interested in is invariant to linear transformations, such as the multivariate normal, the trick is easy and well-known (exemplified in the following example script): <hansl> matrix Omega = unvech({2; 1; 1; 2; 0; 1}) X = mnormal(100,3) matrix K_om = cholesky(Omega)' matrix K = cholesky(inv(mcov(X))) * K_om matrix Z = X * K_om # draw from density with population variance Omega matrix W = X * K # sample variance = Omega printf "Sample covariance matrix for Z:\n%10.4f\n", mcov(Z) printf "Sample covariance matrix for W:\n%10.4f\n", mcov(W) </hansl> In the more general case, in which your joint density is _not_ invariant to linear (or affine) transformations, a solution does not necessarily exist. For example, I think I remember that a bivariate probability distribution such that marginal distributions are Poisson and covariance is negative does not exist (but I may be wrong). In those cases, you'd have to resort to copulas. ------------------------------------------------------- Riccardo (Jack) Lucchetti Dipartimento di Scienze Economiche e Sociali (DiSES) Università Politecnica delle Marche (formerly known as Università di Ancona) r.lucchetti(a)univpm.it http://www2.econ.univpm.it/servizi/hpp/lucchetti -------------------------------------------------------