I want to generate random monthly (m) temperature (T) and Precipitation (P) data considering that both variables are intercorrelated (rTP[m]) The tricky thing is that my random variables that have specific quantitative properties: temperatures are normally distributed, while precipitations follow a log-normal distribution and should be log-transformed

mvrnorm of the package MASS could be used.

mP=c(3.9,3.7,3.9,4.1,4.5,4.7,4.8,4.8,4.4,4.1,4.2,3.9) #log-transformed
sdP=c(0.7,0.8,0.7,0.6,0.4,0.4,0.4,0.5,0.6,1,0.8,0.6)  #log-transformed

for (m in 1:12)
out=mvrnorm(500, mu = c(mT[m],mP[m]), Sigma = matrix(c(sdT[m]*sdT[m],covTP[m],covTP[m],sdP[m]*sdP[m]), ncol = 2), empirical = TRUE)

In this case I generate two random time-series that are inter-correlated which is great.

However, simulated precipitations (simP) are on average higher than the observed one (mP)

plot(exp(mP)-1, type="l", lwd=2, ylim=c(0,250)); points(simP, type="l", lwd=2, lty=2)

I could use rlnorm or rlnorm.rplus to consider than precipitations are log-transformed, but then I have troubles with temperatures that are normally distributed.

My question is: How can I create random sampling for variables that have specific quantitative properties (log-normal and normal distributions)?

  • $\begingroup$ Are you asking how to produce samples from distributions whose multivariate moments will match the moments of the data? $\endgroup$ – whuber Sep 10 '14 at 15:09
  • $\begingroup$ Yes, but within this multivariate data one variable is log-transformed and not the other one... $\endgroup$ – user52463 Sep 10 '14 at 15:22
  • 1
    $\begingroup$ Yes, but the lognormality creates no problems in terms of generating correlated data. In this question, though, you observe that "simulated (values) ... are on average higher than the observed (ones)," which suggests that your concern is not actually about creating correlated variables but rather focuses on how to generate random variables that have specific quantitative properties. If you could edit your question to clarify this point you might increase your chances of getting some really on-target answers. $\endgroup$ – whuber Sep 10 '14 at 15:41

The the problem is in calculating the mean of the original distribution in the plot statement. The correct way is

plot(exp(mP + sdP^2/2), type="l", lwd=2, ylim=c(0,250)); points(simP, type="l", lwd=2, lty=2)

Mean and standard deviation do not commute with logarithm or exponent. In general

log(mean(mydata)) != mean(log(mydata))

and similarly also for sd, exp in any combination of the above You should also make sure that mP and sdP are the mean and sd of the transformed data and not the logarithms of the mean and sd of the raw data. ie you want

mP <- mean(log(P))
stP <- sd(log(P))

You can check the formula for the mean and other statistics of the log normal on the wikipedia

| cite | improve this answer | |
  • $\begingroup$ Thanks! I'm a little bit ashamed, it was a simple mathematical issue... $\endgroup$ – user52463 Sep 11 '14 at 7:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.