如何在 R 中为蒙特卡洛创建更有效的模拟循环答案

【问题标题】：How to create a more efficient simulation loop for Monte Carlo in R如何在 R 中为蒙特卡洛创建更有效的模拟循环
【发布时间】：2012-02-19 00:23:57
【问题描述】：

本练习的目的是创建营养摄入值的人口分布。早期数据中有重复测量，这些已被删除，因此每一行在数据框中都是唯一的人。

我有这段代码，在使用少量数据框行进行测试时效果很好。对于所有 7135 行，它非常慢。我试着给它计时，但是当我的机器上经过的运行时间为 15 小时时，我把它搞砸了。 system.time 结果为 Timing stopped at: 55625.08 2985.39 58673.87。

如果有任何关于加速模拟的 cmets，我将不胜感激：

Male.MC <-c()
for (j in 1:100)            {
for (i in 1:nrow(Male.Distrib))  {
    u2        <- Male.Distrib$stddev_u2[i] * rnorm(1, mean = 0, sd = 1)
    mc_bca    <- Male.Distrib$FixedEff[i] + u2
    temp      <- Lambda.Value*mc_bca+1
    ginv_a    <- temp^(1/Lambda.Value)
    d2ginv_a  <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
    mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var / 2
z <- data.frame(
     RespondentID = Male.Distrib$RespondentID[i], 
     Subgroup     = Male.Distrib$Subgroup[i], 
     mc_amount    = mc_amount,
     IndvWeight   = Male.Distrib$INDWTS[i]/100
     )

Male.MC <- as.data.frame(rbind(Male.MC,z))
    }
}

对于我数据集中的 7135 个观察值中的每一个，都会创建 100 个模拟营养值，然后将其转换回原始测量水平（该模拟使用来自 BoxCox 转换营养值的非线性混合效应模型的结果）。

我不想使用for 循环，因为我读到它们在R 中效率低下，但我对基于apply 的选项了解不够，无法使用它们作为替代方案。 R 正在独立机器上运行，通常这将是运行 Windows 7 变体的标准戴尔型台式机，如果这会影响有关如何更改代码的建议。

更新：为了重现这个进行测试， Lambda.Value=0.4 和 Male.Resid.Var=12.1029420429778 和 Male.Distrib$stddev_u2 是所有观察值的恒定值。

str(Male.Distrib)是

'data.frame':   7135 obs. of  14 variables:
 $ RndmEff     : num  1.34 -5.86 -3.65 2.7 3.53 ...
 $ RespondentID: num  9966 9967 9970 9972 9974 ...
 $ Subgroup    : Ord.factor w/ 6 levels "3"<"4"<"5"<"6"<..: 4 3 2 4 1 4 2 5 1 2 ...
 $ RespondentID: int  9966 9967 9970 9972 9974 9976 9978 9979 9982 9993 ...
 $ Replicates  : num  41067 2322 17434 21723 375 ...
 $ IntakeAmt   : num  33.45 2.53 9.58 43.34 55.66 ...
 $ RACE        : int  2 3 2 2 3 2 2 2 2 1 ...
 $ INDWTS      : num  41067 2322 17434 21723 375 ...
 $ TOTWTS      : num  1.21e+08 1.21e+08 1.21e+08 1.21e+08 1.21e+08 ...
 $ GRPWTS      : num  41657878 22715139 10520535 41657878 10791729 ...
 $ NUMSUBJECTS : int  1466 1100 1424 1466 1061 1466 1424 1252 1061 1424 ...
 $ TOTSUBJECTS : int  7135 7135 7135 7135 7135 7135 7135 7135 7135 7135 ...
 $ FixedEff    : num  6.09 6.76 7.08 6.09 6.18 ...
 $ stddev_u2   : num  2.65 2.65 2.65 2.65 2.65 ...

head(Male.Distrib)是

    RndmEff RespondentID Subgroup RespondentID Replicates IntakeAmt RACE INDWTS    TOTWTS   GRPWTS NUMSUBJECTS TOTSUBJECTS  FixedEff stddev_u2
1  1.343753         9966        6         9966      41067 33.449808    2  41067 120622201 41657878        1466        7135  6.089918  2.645938
2 -5.856516         9967        5         9967       2322  2.533528    3   2322 120622201 22715139        1100        7135  6.755664  2.645938
3 -3.648339         9970        4         9970      17434  9.575439    2  17434 120622201 10520535        1424        7135  7.079757  2.645938
4  2.697533         9972        6         9972      21723 43.340180    2  21723 120622201 41657878        1466        7135  6.089918  2.645938
5  3.531878         9974        3         9974        375 55.660607    3    375 120622201 10791729        1061        7135  6.176319  2.645938
6  6.627767         9976        6         9976      48889 91.480049    2  48889 120622201 41657878        1466        7135  6.089918  2.645938

更新 2：导致 NaN 结果的函数行是

d2ginv_a  <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))

感谢大家的帮助和cmets，也感谢大家的回复速度。

更新：@Ben Bolker 是正确的，即负的 temp 值导致了 NaN 问题。我通过一些测试错过了这一点（在注释掉函数以便只返回 temp 值并调用我的结果数据框 Test 之后）。此代码重现了NaN 问题：

> min(Test)
[1] -2.103819
> min(Test)^(1/Lambda.Value)
[1] NaN

但是将值作为值放入然后运行相同的（？）计算会给我一个结果，所以我在手动计算时错过了这个：

> -2.103819^(1/Lambda.Value) 
[1] -6.419792

我现在有（我认为）使用矢量化的工作代码，而且速度非常快。以防万一其他人有这个问题，我在下面发布工作代码。我必须添加一个最小值以防止计算出现 rnorm 结果放入数据框，这确实减慢了速度，以这种方式创建它们然后使用cbind 非常快。 Male.Distrib 是我的 7135 个观察值的完整数据框，但此代码应该适用于我之前发布的缩减版本（未经测试）。

Min_bca <- ((.5*min(Male.AddSugar$IntakeAmt))^Lambda.Value-1)/Lambda.Value
Test <- Male.Distrib[rep(seq.int(1,nrow(Male.Distrib)), 100), 1:ncol(Male.Distrib)]
RnormOutput <- rnorm(nrow(Test),0,1)
Male.Final <- cbind(Test,RnormOutput)
Male.Final$mc_bca    <- Male.Final$FixedEff + (Male.Final$stddev_u2 *     Male.Final$RnormOutput)
Male.Final$temp      <- ifelse(Lambda.Value*Male.Final$mc_bca+1 > Lambda.Value*Min_bca+1,
                           Lambda.Value*Male.Final$mc_bca+1, Lambda.Value*Min_bca+1)
Male.Final$ginv_a    <- Male.Final$temp^(1/Lambda.Value)
Male.Final$d2ginv_a  <- ifelse(0 > (1-Lambda.Value)*Male.Final$temp^(1/Lambda.Value-2),
                           0, (1-Lambda.Value)*Male.Final$temp^(1/Lambda.Value-2))
Male.Final$mc_amount <- Male.Final$ginv_a + Male.Final$d2ginv_a * Male.Resid.Var / 2

当天的课程：

如果您尝试执行我之前尝试的操作，分布函数似乎不会在循环中重新采样
你不能像我尝试的那样使用max()，因为它返回列中的最大值，而我想要两个值中的最大值。 ifelse 语句是要替换的语句。

【问题讨论】：

我怀疑这将简化为使用replicate 和矩阵+数组数学一次完成所有观察的单行。不过，您能否发布一个可重现的小示例，以便我们为您提供更具体的建议？
使用rbind() 增长对象非常昂贵。您可以更好地在开始时创建一个空数据框（例如，用虚拟变量填充它）并将其填充到循环中。
除了@SachaEpskamp 所说的，不需要内部循环。您使用的所有功能都是矢量化的；充分利用这一点。
我同意@JohnColby。我认为“启动”包中的复制或boot 都可以完成工作，并且可能比您现在所做的更具统计有效性。
根据您希望max(0,...) 执行的操作，您可以使用max(0,...,na.rm=TRUE) 或单独测试(1-Lambda.Value) 和temp 组件。

标签： r loops simulation montecarlo performance

【解决方案1】：

这是解决两个最大速度问题的方法：

我们不是循环观察观察结果 (i)，而是一次计算所有观察结果。
我们使用 replicate，而不是循环遍历 MC 复制 (j)，这是用于此目的的简化 apply。

首先，我们加载数据集并为您正在执行的操作定义一个函数。

Male.Distrib = read.table('MaleDistrib.txt', check.names=F)

getMC <- function(df, Lambda.Value=0.4, Male.Resid.Var=12.1029420429778) {
  u2        <- df$stddev_u2 * rnorm(nrow(df), mean = 0, sd = 1)
  mc_bca    <- df$FixedEff + u2
  temp      <- Lambda.Value*mc_bca+1
  ginv_a    <- temp^(1/Lambda.Value)
  d2ginv_a  <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
  mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var / 2
  mc_amount
}

然后我们复制了很多次。

> replicate(10, getMC(Male.Distrib))
         [,1]      [,2]     [,3]     [,4]      [,5]     [,6]     [,7]     [,8]     [,9]    [,10]
[1,] 36.72374 44.491777 55.19637 23.53442 23.260609 49.56022 31.90657 25.26383 25.31197 20.58857
[2,] 29.56115 18.593496 57.84550 22.01581 22.906528 22.15470 29.38923 51.38825 13.45865 21.47531
[3,] 61.27075 10.140378 75.64172 28.10286  9.652907 49.25729 23.82104 31.77349 16.24840 78.02267
[4,] 49.42798 22.326136 33.87446 14.00084 25.107143 25.75241 30.20490 33.14770 62.86563 27.33652
[5,] 53.45546  9.673162 22.66676 38.76392 30.786100 23.42267 28.40211 35.95015 43.75506 58.83676
[6,] 34.72440 23.786004 63.57919  8.08238 12.636745 34.11844 14.88339 21.93766 44.53451 51.12331

然后您可以重新格式化、添加 ID 等，但这是主要计算部分的想法。祝你好运！

【讨论】：

谢谢约翰，这看起来确实是要走的路，不过我得到了NaN 每次复制的结果，但我不知道为什么。它适用于测试数据，然后在我在完整数据帧上运行时失败。
请注意，用replicate 替换外循环是装饰性的——那里没有速度增益。速度增益来自避免 rbind 和元素操作。
唯一可以明显给出NaN结果的操作是将负数（temp）提高到分数幂（1/Lambda.Value，1/Lambda.Value-2）。从Male.Distrib 发布summary 结果？
糟糕，您确实发布了str（与summary 一样好）。您的标准差是 2.65，因此您可以预期 u2 会定期下降到 -5 或 -6，这可能会使 mc_bca
@BenBolker 不错的提示！我自己从来没有调查过。