如何列绑定对数正态分布的蒙特卡洛模拟结果？

Question

我想对 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)))
}

第 2 部分代码：

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

问题是我一直在使用定义的函数时出错。

Answer 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-07-11 创建

如何列绑定对数正态分布的蒙特卡洛模拟结果？

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

r

montecarlo

lapply

do.call