【问题标题】:simulate a linear model 100 times模拟线性模型 100 次
【发布时间】:2013-01-28 00:36:23
【问题描述】:

我需要模拟 n=100 次线性模型,但在 R 命令中迷失了方向。

我还在学习统计和 R 的基础知识,我对这个练习有点困惑:

我需要使用 OLS 将基本线性模型复制 100 次并收集 N 个估计值,以便执行一致性和效率测试。我试图通过这种方式解决问题:

a <- 3
B <- 0.5
C <- -0.7

for (i in 1:100){
   x1[i] <- rnorm(200, mean=0, sd=1)
   x2[i] <- rnorm(200, mean=0, sd=1) 
   e[i] <- rnorm(200, mean=0, sd=1)
   y1[i] <- a+(B*x1[i])+(C*x2[i])+e[i] 

   y<- lm(y1[i]~x1[i]+x2[i]))
   results <-data.frame(coef(y))
}

但是 R 一直告诉我有错误。有人可以帮我解决这个问题吗?

【问题讨论】:

  • 欢迎来到简历,卡尔!因为您的问题实际上源于您如何在R 中分配和寻址数组,所以它更适合 SO,我们将在其中迁移它。但是,通过searching CV——或者只是阅读足够多的问题——你会在这里找到很多个正是这种模拟的例子,带有完整的代码和输出,你可以立即开始使用。跨度>
  • “CV”代表什么?
  • CrossValidated: stats.stackexchange.com ,用于统计的 StackExchange 站点

标签: r


【解决方案1】:

类似:

a <- 3
B <- 0.5
C <- -0.7

results <- matrix(nrow=100,ncol=3)
for (i in 1:100){
   x1 <- rnorm(200, mean=0, sd=1)
   x2 <- rnorm(200, mean=0, sd=1) 
   e <- rnorm(200, mean=0, sd=1)
   y1 <- a+B*x1+C*x2+e 

   y<- lm(y1~x1+x2)
   results[i,] <- coef(y)
}

这假设您只需要保存每次运行的系数。更优雅的解决方案是:

simfun <- function(a=3,B=0.5,C=-0.7,n=200,x1.sd=1,x2.sd=1,e.sd=1) {
   x1 <- rnorm(n, mean=0, sd=x1.sd)
   x2 <- rnorm(n, mean=0, sd=x2.sd) 
   e <-  rnorm(n, mean=0, sd=e.sd)
   y1 <- a+B*x1+C*x2+e 
   data.frame(x1,x2,y1)
}

statfun <- function(d) {
    coef(lm(y1~x1+x2,data=d))
}

library(plyr)
raply(100,statfun(simfun()))

【讨论】:

  • 您可能会摆脱 plyr 并只使用 replicate - 但每个人都有自己的。
  • 是的。我喜欢raply 方便的内置进度条,但这里既没有必要也没有显示。
【解决方案2】:
    a <- 3
    B <- 0.5
    C <- -0.7
    sims <- 100

    #initialize a data frame to collect results
    df <- data.frame(matrix(ncol = 3, nrow = sims))
    colnames(df) <- c('a', 'B' , 'C')

    for(i in 1:sims){
    ##vectors each  200 long 
    x1 <- rnorm(200)
    x2 <- rnorm(200)
    e <- rnorm(200)

    y <- a + B*x1 + C*x2 +e
    #collect results for each itter
    df[i,] <- data.frame(t(lm(y ~x1 + x2)$coeff))
    }

    #results
    df

【讨论】:

  • 与本的解决方案相同——我可能应该删除它——一旦我弄清楚如何
猜你喜欢
  • 2023-03-08
  • 1970-01-01
  • 2018-03-17
  • 1970-01-01
  • 1970-01-01
  • 2016-11-28
  • 1970-01-01
  • 1970-01-01
  • 2017-01-08
相关资源
最近更新 更多