如何在R中使用MonteCarlo()正确运行蒙特卡罗模拟

Question

我需要通过多次迭代实验来运行蒙特卡罗模拟。MonteCarlo包的功能似乎很好地解决了这个问题。该实验包括为每次迭代创建4个矩阵(a、b、c和d)。

矩阵a：

random.sample <- sample(1:100, 1000, replace = TRUE) 
a = matrix(random.sample, nrow = 10, ncol = 10, byrow = TRUE)        

矩阵b：

b <- apply(a[, 1:10], 2, SMA, n=3)

矩阵c：

c <- pmax(a, b)

概率序列：

probability.sequence <- seq(0.1, 0.90, 0.1)

矩阵d：

d <- apply(c, 2, quantile, probs = probability.sequence, na.rm = TRUE)

我想迭代矩阵a到d，给定的次数，n，每一次迭代应该产生一个新的矩阵d，其中包含相应矩阵c的百分位数计算，我想把迭代过程的结果汇总起来，从而从所有矩阵d得到90%百分位数的均值、方差、置信上限和置信下限。

(近似)期望产出：

str(monte.carlo.aggregate.results)
num [1, 1:4] 65 15 85 45 
monte.carlo.aggregate.results         
     mean variance Upper CL Lower CL
[1,]   65       15       85       45

我读过?MonteCarlo，并且很难将矩阵a到d的创建封装在一个函数中，以便与包含在MonteCarlo包中的MonteCarlo()函数一起使用。是否有人建议将此过程设置为与MonteCarlo()一起使用，或者是否有涉及自定义for循环的替代解决方案？

Answer 1

如果您不坚持使用MonteCarlo包(这对我来说是新的，看起来很有趣)，那么您可以将代码打包到一个函数中。

首先，创建一个只需90%行的随机矩阵.然后，在指定了数量之后，您可以执行迭代并使用lapply()将它们打包到一个列表中。最后，只需使用(四舍五入)统计数据创建另一个基准表。

就像这样：

set.seed(54897)  # setting seed for sake of reproducibility
probability.sequence <- seq(0.1, 0.90, 0.1)

# iteration function
mx <- function(X, n) {
  random.sample <- sample(1:1e2, 1e3, replace = TRUE) 
  a = matrix(random.sample, nrow = 10, ncol = 10, byrow = TRUE)  
  b <- apply(a[ ,1:10], 2, SMA, n = 3)
  c <- pmax(a, b)
  d <- apply(c, 2, quantile, probs = probability.sequence,  na.rm = TRUE)
  return(d[9, ])  # returns only the desired 90 % row
}

# amount of iterations
amount <- 1e3 

# iteration
mx.lst <- lapply(1:amount, mx)

# matrix for statistics
r = matrix(NA, nrow = 1, ncol = 4, 
           dimnames = list(NULL, c("mean", "variance", "Upper CL", "Lower CL")))
r[, 1] <- (mean(unlist(lapply(mx.lst, mean))))  # mean
r[, 2]  <- sqrt(var(unlist(lapply(mx.lst, mean))))  # variance
r[, 3]  <- r[, 1] - 1.96 * sqrt(r[, 2])  # lower CI 95 %
r[, 4]  <- r[, 1] + 1.96 * sqrt(r[, 2])  # upper CI 95 %

# output
(monte.carlo.aggregate.results <- round(r, 0))
#      mean variance Upper CL Lower CL
# [1,]   82        3       78       86

str(monte.carlo.aggregate.results)
# num 1, 1:4] 82 3 78 86
# - attr(*, "dimnames")=List of 2
# ..$ : NULL
# ..$ : chr [1:4] "mean" "variance" "Upper CL" "Lower CL"

# speed test
amount <- 1e4
system.time(mx.lst <- lapply(1:amount, mx))
#     user      system     elapsed 
#    48.21        0.00       48.31 

注:请检查自己的信心区间使用的公式是否符合您的需要。