-1

我正在做一个模拟研究，我写了下面的 R 代码。有没有办法在不使用两个循环的情况下编写这段代码for，或者让它更高效（运行得更快）？

S = 10000
n = 100
v = c(5,10,50,100)
beta0.mle = matrix(NA,S,length(v))  #creating 4 S by n NA matrix 
beta1.mle = matrix(NA,S,length(v))
beta0.lse = matrix(NA,S,length(v))
beta1.lse = matrix(NA,S,length(v))
for (j in 1:length(v)){
  for (i in 1:S){
    set.seed(i)
    beta0 = 50
    beta1 = 10
    x = rnorm(n)
    e.t = rt(n,v[j])
    y.t = e.t + beta0 + beta1*x
    func1 = function(betas){
      beta0 = betas[1]
      beta1 = betas[2]
      sum = sum(log(1+1/v[j]*(y.t-beta0-beta1*x)^2))
      return((v[j]+1)/2*sum)
    }
    beta0.mle[i,j] = nlm(func1,c(1,1),iterlim = 1000)$estimate[1]
    beta1.mle[i,j] = nlm(func1,c(1,1),iterlim = 1000)$estimate[2]
    beta0.lse[i,j] = lm(y.t~x)$coef[1]
    beta1.lse[i,j] = lm(y.t~x)$coef[2]
  }
}

func1第二个循环内的函数for用于nlm函数（当错误分布为 t 时查找 mle）。我想parallel在 R 中使用包，但我没有找到任何有用的功能。

4

1 回答 1

2

让任何东西在 R 中运行得更快的关键是用for矢量化函数（例如apply系列）替换循环。nlm此外，对于任何编程语言，您应该寻找使用相同参数多次调用昂贵函数（例如 ）的地方，并查看可以存储结果的位置，而不是每次都重新计算。

在这里，我像您一样通过定义参数开始。同样beta0，beta1自始至终50，10我也将在这里定义这些。

S <- 10000
n <- 100
v <- c(5,10,50,100)
beta0 <- 50
beta1 <- 10

接下来我们将func1在循环外定义以避免每次都重新定义。 func1现在有两个额外的参数，v因此y.t可以使用新值调用它。

func1 <- function(betas, v, y.t, x){
  beta0 <- betas[1]
  beta1 <- betas[2]
  sum <- sum(log(1+1/v*(y.t-beta0-beta1*x)^2))
  return((v+1)/2*sum)
}

现在我们实际上做了真正的工作。我们使用嵌套的应用语句，而不是嵌套循环。外部lapply将为 的每个值创建一个列表，v内部vapply将为您想要为每个值获取的四个值（beta0.mle、beta1.mle、beta0.sle、beta1.lse）创建一个矩阵S。

values <- lapply(v, function(j) vapply(1:S, function(s) {
  # This should look familiar, it is taken from your code
  set.seed(s)
  x <- rnorm(n)
  e.t <- rt(n,j)
  y.t <- e.t + beta0 + beta1*x
  # Rather than running `nlm` and `lm` twice, we run it once and store the results
  nlmmod <- nlm(func1,c(1,1), j, y.t, x, iterlim = 1000)
  lmmod <- lm(y.t~x)
  # now we return the four values of interest
  c(beta0.mle = nlmmod$estimate[1],
    beta1.mle = nlmmod$estimate[2],
    beta0.lse = lmmod$coef[1],
    beta1.lse = lmmod$coef[2])
}, numeric(4)) # this tells `vapply` what to expect out of the function
)

最后，我们可以将所有内容重新组织成四个矩阵。

beta0.mle <- vapply(values, function(x) x["beta0.mle", ], numeric(S))
beta1.mle <- vapply(values, function(x) x["beta1.mle", ], numeric(S))
beta0.lse <- vapply(values, function(x) x["beta0.lse.(Intercept)", ], numeric(S))
beta1.lse <- vapply(values, function(x) x["beta1.lse.x", ], numeric(S))

S最后一点，根据您使用索引设置种子的原因，可以重新组织它以更快地运行。如果知道用什么种子来生成你x的种子很重要，rnorm那么这可能是我能做的最好的事情。但是，如果您这样做只是为了确保您的所有值v都在相同的值上进行测试，x那么我们可以进行更多的重组，这可能会提高使用replicate.

于 2016-11-04T16:49:31.687 回答