How to column bind monte-carlo simulation results for log-normal distribution?

Question

I want perform log-normal monte-carlo simulations for 12 iterations with location and shape of distribution varying for each iteration. As an output of the following code I want each simulation result in a different columns as a data frame. For instance, if I want to perform 10,000 simulations for each set of location and shape, the resultant data frame will have 12 columns and 10000 rows. I wrote the code in two parts where I first defined the function for lapply and then tried to perform simulation with results from each iteration being stored into a different column. Part 1 code:

myfxn <- function(i,s,location,shape){
cvsq=1
  i <-seq(5,60,5)
  s <-i*cvsq
  location <- log(i^2 / sqrt(s^2 + i^2))
  shape <- sqrt(log(1 + (s^2 / i^2)))
}

Part 2 code:

DF1 <- do.call(cbind, lapply(seq(5,60,5), function(myfxn) setNames(data.frame(rlnorm(n=10000, location, shape)), i)))
head(DF1)

The issue is I keep getting errors with the defined function.

Answer 1

You seem to be confused about your function declaration. Your function looks as though it should take 4 parameters (i, s, location and shape). However, inside the function these inputs are never used - they are all created from scratch according to the value of cvsq. They are therefore redundant as inputs to the function.

Inside your lapply, you are passing the values seq(5,60,5) to your function, but it isn't clear what these are for if the variables i, s, location and shape are all calculated anew every time the function is called. What is it you want your function to do with these numbers?

My guess is that you want to feed a single number into your function and have it alter these variables. If so, the only variable it makes sense to have being passed as a single parameter is cvsq, since altering cvsq alters the value of all the other variables. If this is wrong, then it's just not clear what you're trying to do.

If I'm on the right track, your code would look like this:

myfxn <- function(cvsq)
{
  i <-seq(5,60,5)
  s <-i*cvsq
  location <- log(i^2 / sqrt(s^2 + i^2))
  shape <- sqrt(log(1 + (s^2 / i^2)))
  rlnorm(n = 10000, location, shape)
}

my_cvsq <- seq(5, 60, 5)
df <- as.data.frame(do.call(cbind, lapply(my_cvsq, myfxn)))

head(df)
#>           V1        V2         V3        V4         V5        V6           V7
#> 1  0.1786204 1.2587228 0.04588056 1.1925041 0.01985825 0.6924621 0.0008467456
#> 2 28.6147331 0.1243106 0.08306054 0.5057983 4.65819135 3.0851280 0.0429381886
#> 3  0.6666314 0.5546972 2.83435439 0.6512173 0.04994503 6.4365528 0.3779170170
#> 4  7.7163909 4.5954509 0.26637679 0.4085584 3.31589074 0.4894957 1.5195320720
#> 5  1.0931151 0.0575365 0.80689457 0.5033203 1.35200408 0.1427171 2.4058576253
#> 6  0.2573164 7.6190362 1.01222386 1.8598690 5.91525940 0.5600217 0.7858050404
#>            V8         V9        V10         V11       V12
#> 1 1.263350450 0.09485695 0.02615021 0.331885430 2.9109114
#> 2 0.004408080 0.14120821 0.05990997 0.053248650 0.1013357
#> 3 0.888602537 1.32448399 1.17625454 0.766153864 0.0196438
#> 4 5.254249758 0.01337194 0.08463637 2.996341976 0.8942556
#> 5 0.004375628 0.05052838 0.51024824 0.533356862 1.6166927
#> 6 0.985499853 0.09939285 0.86272478 0.001223073 1.0652496

And we can see what your 12 samples look like by plotting their density:

d <- density(df[,1], from = 0, to = 10)
plot(d$x, d$y, type = "l", xlim = c(0, 10), ylim = c(0, 0.7))
for(i in 2:12){
  d <- density(df[,i], from = 0, to = 10)
  lines(d$x, d$y, col = i)
}

^{Created on 2020-07-11 by the reprex package (v0.3.0)}