这个练习的目的是创建营养素摄入值的人口分布。早期的数据中有重复的测量数据,已被删除,所以数据框中的每一行都是一个独特的人。
我有这段代码,在少量数据框行测试时运行得非常好。但对于全部 7135 行,速度非常慢。我尝试计时,但在我的机器上经过的时间达到了15小时,我不得不强制退出。system.time 的结果为 Timing stopped at: 55625.08 2985.39 58673.87。
我希望能够得到任何关于如何加快模拟速度的建议。
针对我的数据集中的7135个观察值,每个观察值都会创建100个模拟营养价值,然后将其反向转换为原始测量水平(该模拟使用来自非线性混合效应模型在BoxCox转换营养价值上的结果)。
我不想使用for循环,因为我读到它们在R中效率低下,但我不了解基于apply的选项足以用作替代。 R正在独立机器上运行,通常这将是运行Windows 7变体的标准Dell式台式机,如果这影响更改代码的建议。
更新:为了进行测试, Lambda.Value = 0.4,Male.Resid.Var = 12.1029420429778,而且Male.Distrib$stddev_u2是所有观察值上的恒定值。
str(Male.Distrib)是:
更新2:导致
我有这段代码,在少量数据框行测试时运行得非常好。但对于全部 7135 行,速度非常慢。我尝试计时,但在我的机器上经过的时间达到了15小时,我不得不强制退出。system.time 的结果为 Timing stopped at: 55625.08 2985.39 58673.87。
我希望能够得到任何关于如何加快模拟速度的建议。
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变体的标准Dell式台式机,如果这影响更改代码的建议。
更新:为了进行测试, 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)
is
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))
感谢大家的帮助和评论,也感谢回复的速度。
更新:@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
我现在有可用的代码,使用了向量化技术,速度非常快。以防其他人遇到同样的问题,我在下面发布可用的代码。我已经添加了一个最小值来防止计算中出现<0的问题。感谢所有帮助过我的人,还有咖啡。我尝试将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
和矩阵+数组数学运算来将其简化为一行代码,以便一次性处理所有观测数据。不过,你能否提供一个小的可重现示例,这样我们就可以给出更具体的建议? - John Colbyrbind()
函数来增加对象的成本很高。您可以在开始时创建一个空的数据框(例如,用虚拟变量填充它),然后在循环中填充它,这样做会更好。 - Sacha Epskampboot
函数来完成工作,这样做可能比你现在所做的更具统计学意义。 - IRTFMmax(0,...)
的功能,您可以选择使用max(0,...,na.rm=TRUE)
或分别测试(1-Lambda.Value)
和temp
组件。 - Ben Bolker