在数据帧的行上迭代 mvrnorm 的复杂函数

编程问答 2022-04-29

问题描述

我有这样的数据：

data = data.frame(a.coef = c(.14,.15,.16),b.coef = c(.4,.5,.6),a.var = c(0.0937,0.0934,0.0945),b.var = c(0.00453,0.00564,0.00624),ab.cov = c(0.000747,0.000747,0.000747))

我想在数据集的每一行上运行以下代码（来源：http://www.quantpsy.org/medmc/medmc.htm）。

require(MASS)
a = data$a.coef
b = data$b.coef
rep = 10000
conf = 95
pest = c(a,b)
acov <- matrix(c(data$a.var,data$ab.cov,data$b.var),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical = FALSE)
ab <- mcmc[,1] * mcmc[,2]
low = (1 - conf / 100) / 2
upp = ((1 - conf / 100) / 2) + (conf / 100)
LL = quantile(ab,low)
UL = quantile(ab,upp)
LL4 = format(LL,digits = 4)
UL4 = format(UL,digits = 4)

我创建了一个相对简单的函数，它将数据和行号作为输入：

MCMAM <- function(data_input,row_number) {
  data = data_input[row_number,]
  a = data[["a.coef"]]
  b = data[["b.coef"]]
  rep = 10000
  conf = 95
  pest = c(a,b)
  acov <- matrix(c(data[["a.var"]],data[["ab.cov"]],data[["b.var"]]),2)
  require(MASS)
  mcmc <- mvrnorm(rep,empirical = FALSE)
  ab <- mcmc[,2]
  low = (1 - conf / 100) / 2
  upp = ((1 - conf / 100) / 2) + (conf / 100)
  LL = quantile(ab,low)
  UL = quantile(ab,upp)
  return(c(LL,UL))
}

MCMAM(data,1)

    2.5%      97.5% 
-0.1901272  0.3104614

但如果有办法摆脱行规范，让函数逐行运行数据集并将输出保存到数据集中的新列，那就太好了。

我一直在试验 for 循环和 apply 函数，但没有取得任何成功，主要是因为 matrix() 和 mvrnorm() 函数都采用值而不是向量。

解决方法

我们可以使用lapply

do.call(rbind,lapply(seq_len(nrow(data)),MCMAM,data_input = data))

-输出

    2.5%     97.5%
[1,] -0.1832449 0.3098362
[2,] -0.2260856 0.3856575
[3,] -0.2521126 0.4666583

或使用 rowwise

library(dplyr)
library(tidyr)
data %>% 
     rowwise %>% 
     mutate(new = list(MCMAM(cur_data(),1))) %>%
     unnest_wider(new)
# A tibble: 3 x 7
#  a.coef b.coef  a.var   b.var   ab.cov `2.5%` `97.5%`
#   <dbl>  <dbl>  <dbl>   <dbl>    <dbl>  <dbl>   <dbl>
#1   0.14    0.4 0.0937 0.00453 0.000747 -0.185   0.309
#2   0.15    0.5 0.0934 0.00564 0.000747 -0.219   0.396
#3   0.16    0.6 0.0945 0.00624 0.000747 -0.259   0.472
dplyrmassrr