**P**aulo (from the Instituto de Matemática e Estatística, Universidade de São Paulo, Brazil) has posted an answer to my earlier question both as a comment on the ‘Og and as a solution on StackOverflow (with a much more readable LaTeX output). His solution is based on the observation that the multidimensional log-normal distribution still allows for closed form expressions of both the mean and the variance and that those expressions can further be inverted to impose the pair *(μ,Σ)* on the log-normal vector. In addition, he shows that the only constraint on the covariance matrix is that the covariance *σ _{ij}* is larger than

*-μ*. Very neat!

_{i}μ_{j.}**I**n the meanwhile, I corrected my earlier R code on the gamma model, thanks to David Epstein pointing a mistake in the resolution of the moment equation and I added the constraint on the covariance, already noticed by David in his question. Here is the full code:

sol=function(mu,sigma){ solub=TRUE alpha=rep(0,3) beta=rep(0,2) beta[1]=mu[1]/sigma[1] alpha[1]=mu[1]*beta[1] coef=mu[2]*sigma[1]-mu[1]*sigma[3] if (coef<0){ solub=FALSE}else{ beta[2]=coef/(sigma[1]*sigma[2]-sigma[3]^2) alpha[2]=sigma[3]*beta[2]/sigma[1]^2 alpha[3]=mu[2]*beta[2]-mu[1]*alpha[2] if (alpha[3] } list(solub=solub,alpha=alpha,beta=beta) } mu=runif(2,0,10);sig=c(mu[1]^2/runif(1),mu[2]^2/runif(1));sol(mu,c(sig,runif(1,max(-sqrt(prod(sig)), -mu[1]*mu[2]),sqrt(prod(sig)))))

and I did not get any FALSE outcome when running this code several times.