**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 *-μ*_{i}μ_{j.}. Very neat!

**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.