如何为对数正态分布列绑定蒙特卡罗模拟结果？答案

【问题标题】：How to column bind monte-carlo simulation results for log-normal distribution?如何为对数正态分布列绑定蒙特卡罗模拟结果？
【发布时间】：2020-07-11 20:41:07
【问题描述】：

我想对 12 次迭代执行对数正态蒙特卡罗模拟，每次迭代的位置和分布形状都不同。作为以下代码的输出，我希望每个模拟结果在不同的列中作为数据框。例如，如果我想为每组位置和形状执行 10,000 次模拟，则生成的数据框将有 12 列和 10000 行。我分两部分编写代码，首先为 lapply 定义函数，然后尝试执行模拟，将每次迭代的结果存储到不同的列中。第 1 部分代码：

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)))
}

第二部分代码：

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

问题是我不断收到定义函数的错误。

【问题讨论】：

在第 2 部分代码中，您将 myfxn 作为匿名函数应用
另外，在函数myfxn中，指定的参数是i、s、location、shape，而那些是在函数中创建的对象？

标签： r lapply montecarlo do.call

【解决方案1】：

您似乎对函数声明感到困惑。您的函数看起来好像应该采用 4 个参数（i、s、位置和形状）。但是，在函数内部，这些输入从未使用过——它们都是根据cvsq 的值从头开始创建的。因此，它们作为函数的输入是多余的。

在您的lapply 中，您将值seq(5,60,5) 传递给您的函数，但是如果每次函数执行时都重新计算变量i、s、位置和形状，则不清楚这些值是什么叫。你希望你的函数用这些数字做什么？

我的猜测是你想将一个数字输入你的函数并让它改变这些变量。如果是这样，作为单个参数传递的唯一变量是cvsq，因为更改cvsq 会更改所有其他变量的值。如果这是错误的，那么就不清楚您要做什么。

如果我在正确的轨道上，您的代码将如下所示：

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

我们可以通过绘制它们的密度来查看您的 12 个样本的样子：

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)
}

^{由reprex package (v0.3.0) 于 2020 年 7 月 11 日创建}

【讨论】：