使用 apply for 模拟而不是嵌套 for 循环

Using apply for simulation instead of nested for loops

我在 R 中复制了一个最初在 Stata 中完成的模拟。我使用了 'for' 循环,因为这是我知道如何使它工作的唯一方法。 运行 需要相当长的时间,所以我想使用 'apply' 命令之一来查看它是否更快,但我不知道该怎么做。有人可以帮忙吗?这是代码:

simdiffuse <- function(a, b, c, d) {

  endo <- 1/a        # innovation endogenous effect
  endomacro <- 1/b   # category endogenous effect
  appeal <- c        # innovation's ex ante appeal
  ninnov <- d        # number of innovations in category 

  results <- data.frame(catdensity = rep(0:ninnov, each = 25), t = 1:25, endo = endo, endomacro = endomacro, appeal = appeal, adopt = NA)    

  prop <- rnorm(1000)
  diff <- data.frame(prop)
  diff$adopt <- 0
  diff$adopt[1:5] <- 1

  for (catdensity in 0:ninnov) {
    diff$adopt <- 0
    diff$adopt[1:5] <- 1

    for (t in 1:25) {
      results[results$catdensity == catdensity & results$t == t,]$adopt <- mean(diff$adopt)
      for (obs in 1:nrow(diff)) {
        if(appeal+(mean(diff$adopt)*endo)+(catdensity*endomacro) > rnorm(1, diff[obs,]$prop)) diff[obs,]$adopt <- 1
      }
    }
  }
  return(results)
}

results <- simdiffuse(.2, 20, -3, 60)

您可以使用 data.table 提高函数的速度。但是,您仍然必须使用 for 循环(这不是一件坏事)。

library(data.table)
simdiffuse <- function(a, b, c, d) {

  endo <- 1/a        # innovation endogenous effect
  endomacro <- 1/b   # category endogenous effect
  appeal <- c        # innovation's ex ante appeal
  ninnov <- d        # number of innovations in category 

  results <- data.table(catdensity = rep(0:ninnov, each = 25), t = 1:25, 
                        endo = endo, endomacro = endomacro, appeal = appeal, 
                        adopt = as.numeric(NA))    


  for (cc in 0:ninnov) {
    diff <- data.table(prop = rnorm(1000), adopt = c(rep(1,5), rep(0, 995)))
    for (tt in 1:25) {
      results[catdensity == cc & t == tt, adopt := diff[, mean(adopt)]]
      diff[, rr := rnorm(1, prop), by="prop"]
      diff[appeal + mean(adopt) * endo + cc * endomacro > rr, adopt := 1]
    }
  }
  return(results)
}

results <- simdiffuse(.2, 20, -3, 60)

我没有足够的时间等到你的代码被执行,但这是我得到的 adoptt 的图表: