我無法按照我想要的方式保存 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 個值的事實。我想為in 的Bhat_OLS每個值設定一個值，以便我可以針對.pPBhat_OLSp
uj5u.com熱心網友回復：
您可以將結果寫入包含兩列的資料框中，其中包含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)
uj5u.com熱心網友回復：
您遍歷概率的事實使得索引變得困難。您可以seq(P)改為回圈使用和 subset 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)"
轉載請註明出處，本文鏈接：https://www.uj5u.com/gongcheng/386699.html