如何将此 for 循环的结果保存为向量而不是单个值？

Question

我无法以我想要的方式保存 for 循环的结果。

我当前正在运行的循环如下所示：

# Setup objects
n = 100
R = (1:1000)
P = seq(-.9, .9, .1)
betahat_OLS = rep(NA, 1000)
Bhat_OLS = rep(NA, 19)

# Calculate betahat_OLS for each p in P and each r in R
for (p in P) {
  for (r in R) {
    # Simulate data
    v = rnorm(n)
    e = rnorm(n)
    z = rnorm(n)
    u = p*v+e
    x = z+v
    y = 0*x+u
    #Calculate betahat_OLS
    betahat_OLS[r] = sum(x*y)/sum(x^2)
  }
  #Calculate Bhat_OLS
  Bhat_OLS = sum(betahat_OLS)/1000-0
}

# Make a scatterplot with p on the x-axis and Bhat_OLS on the y-axis
plot(P, Bhat_OLS)

循环似乎工作正常，除了因为我想最终得到 19 个 Bhat_OLS 值并且目前只得到 1 个值。我想在P 中为p 的每个值设置一个Bhat_OLS 值，以便我可以针对p 绘制Bhat_OLS。

Answer 1

您可以将结果写入具有两列的数据框中，其中包含P 和Bhat_OLS。

# Setup objects
n = 100
R = (1:1000)
P = seq(-.9, .9, .1)
betahat_OLS = rep(NA, 1000)
Bhat_OLS = rep(NA, 19)

# initialize result data frame
results <- data.frame(matrix(ncol = 2, nrow = 0, 
                      dimnames = list(NULL, c("P", "Bhat_OLS"))))

# Calculate betahat_OLS for each p in P and each r in R
for (p in P) {
    for (r in R) {
        # Simulate data
        v = rnorm(n)
        e = rnorm(n)
        z = rnorm(n)
        u = p*v+e
        x = z+v
        y = 0*x+u
        #Calculate betahat_OLS
        betahat_OLS[r] = sum(x*y)/sum(x^2)
    }
    #Calculate Bhat_OLS
    Bhat_OLS = sum(betahat_OLS)/1000-0
    
    # insert P and Bhat_OLS into results
    results[nrow(results) + 1,] = c(p, Bhat_OLS)
}

# Make a scatterplot with p on the x-axis and Bhat_OLS on the y-axis
plot(results$P, results$Bhat_OLS)

Answer 2

循环遍历概率的事实使索引变得困难。你可以循环遍历seq(P) 和子集P[i]。另外，最后你需要Bhat_OLS[i]。然后就可以了。

# Setup objects
n <- 100
R <- (1:1000)
P <- seq(-.9, .9, .1)
betahat_OLS <- rep(NA, length(R))
Bhat_OLS <- rep(NA, length(P))

set.seed(42)  ## for sake of reproducibility

# Calculate betahat_OLS for each p in P and each r in R
for (i in seq(P)) {
  for (r in R) {
    # Simulate data
    v <- rnorm(n)
    e <- rnorm(n)
    z <- rnorm(n)
    u <- P[i]*v + e
    x <- z + v
    y <- 0*x + u
    #Calculate betahat_OLS
    betahat_OLS[r] <- sum(x*y)/sum(x^2)
  }
  #Calculate Bhat_OLS
  Bhat_OLS[i] <- sum(betahat_OLS)/1000 - 0
}

# Make a scatterplot with p on the x-axis and Bhat_OLS on the y-axis
plot(P, Bhat_OLS, xlim=c(-1, 1))

替代解决方案vapply

以更 R-ish 的方式（现在它更 C-ish），您可以在函数 sim() 中定义模拟，并将 vapply 用于外部循环。 （其实也适用于内循环，不过我已经测试过了，这样更快。）

sim <- \(p, n=100, R=1:1000) {
  r <- rep(NA, max(R))
  for (i in R) {
    v <- rnorm(n)
    e <- rnorm(n)
    z <- rnorm(n)
    u <- p*v + e
    x <- z + v
    y <- 0*x + u
    r[i] <- sum(x*y)/sum(x^2)
  }
  return(sum(r/1000 - 0))
}

set.seed(42)
Bhat_OLS1 <- vapply(seq(-.9, .9, .1), \(p) sim(p), 0)

stopifnot(all.equal(Bhat_OLS, Bhat_OLS1))

注意：

R.version.string
# [1] "R version 4.1.2 (2021-11-01)"