R语言 蒙特卡洛运行错误求大神帮忙看看

最后print出问题 而且print怎么又两个b1_se?

> # Set the true values for the intercept and slope
> b1=-7
> b2=25
> b3=300
>
> # READ IN THE DATA FOR THE SEEDS - TO ENSURE RESULTS CAN BE REPLICATED
> SEEDS <- read.csv("SEEDS.csv")
>
>
> # SET YOUR DATASET NUMBER (DIFFERENT FOR EACH STUDENT)
> dataset_number <- 30
>
>
> # Now run a loop which creates the data, runs the regressions and stores the results
> # Note - SEEDS have been added, as discussed in class, to ensure replicability
> for (MC_rep in 1:MC_rep) {
+
+ ### FIRST GENERATE THE X VARIABLES FOR DIFFERENT NUMBERS OF OBSERVATIONS
+ x10 <- (1:30)
+ x11 <- (1:60)
+ x12 <- (1:120)
+ x13 <- (1:250)
+ x14 <- (1:500)
+ x15 <- (1:1000)
+ x20 <- runif(30,min=1000,max=25000)
+ x21 <- runif(60,min=1000,max=25000)
+ x22 <- runif(120,min=1000,max=25000)
+ x23 <- runif(250,min=1000,max=25000)
+ x24 <- runif(500,min=1000,max=25000)
+ x25 <- runif(1000,min=1000,max=25000)
+
+ ### NOW DRAW THE RESIDUALS
+ set.seed(SEEDS[MC_rep,dataset_number]+0); u0 <- runif(rnorm(30))
+ set.seed(SEEDS[MC_rep,dataset_number]+1); u1 <- runif(rnorm(60))
+ set.seed(SEEDS[MC_rep,dataset_number]+2); u2 <- runif(rnorm(120))
+ set.seed(SEEDS[MC_rep,dataset_number]+3); u3 <- runif(rnorm(250))
+ set.seed(SEEDS[MC_rep,dataset_number]+4); u4 <- runif(rnorm(500))
+ set.seed(SEEDS[MC_rep,dataset_number]+5); u5 <- runif(rnorm(1000))
+
+ ### ADD THE RESIDUALS TO THE KNOWN DATA GENERATING PROCESS (y = b1 + b2*x)
+ y0 <- b1 + b2*x10+b3*x20 + u0
+ y1 <- b1 + b2*x11+b3*x21 + u1
+ y2 <- b1 + b2*x12+b3*x22 + u2
+ y3 <- b1 + b2*x13+b3*x23 + u3
+ y4 <- b1 + b2*x14+b3*x24 + u4
+ y5 <- b1 + b2*x15+b3*x25 + u5
+
+ ### ESTIMATRE THE MODELS
+ fit0 <- summary( lm( y0 ~ 1 + x10+x20 ) )
+ fit1 <- summary( lm( y1 ~ 1 + x11+x21 ) )
+ fit2 <- summary( lm( y2 ~ 1 + x12+x22 ) )
+ fit3 <- summary( lm( y3 ~ 1 + x13+x23 ) )
+ fit4 <- summary( lm( y4 ~ 1 + x14+x24 ) )
+ fit5 <- summary( lm( y5 ~ 1 + x15+x25 ) )
+
+ ### STORE THE BIAS E.G. THE DIFFERENCE BETWEEN THE KNOWN VALUE OF A
+ ### COEFFICIENT AND THE ESTIMATED VALUE
+ b1_bias[MC_rep,1] = fit0$coef[1,1]-b1
+ b1_bias[MC_rep,2] = fit1$coef[1,1]-b1
+ b1_bias[MC_rep,3] = fit2$coef[1,1]-b1
+ b1_bias[MC_rep,4] = fit3$coef[1,1]-b1
+ b1_bias[MC_rep,5] = fit4$coef[1,1]-b1
+ b1_bias[MC_rep,6] = fit5$coef[1,1]-b1
+
+ b2_bias[MC_rep,1] = fit0$coef[2,1]-b2
+ b2_bias[MC_rep,2] = fit1$coef[2,1]-b2
+ b2_bias[MC_rep,3] = fit2$coef[2,1]-b2
+ b2_bias[MC_rep,4] = fit3$coef[2,1]-b2
+ b2_bias[MC_rep,5] = fit4$coef[2,1]-b2
+ b2_bias[MC_rep,6] = fit5$coef[2,1]-b2
+
+ b3_bias[MC_rep,1] = fit0$coef[3,1]-b3
+ b3_bias[MC_rep,2] = fit1$coef[3,1]-b3
+ b3_bias[MC_rep,3] = fit2$coef[3,1]-b3
+ b3_bias[MC_rep,4] = fit3$coef[3,1]-b3
+ b3_bias[MC_rep,5] = fit4$coef[3,1]-b3
+ b3_bias[MC_rep,6] = fit5$coef[3,1]-b3
+
+ ### STORE THE STANDARD ERRORS OF THE ESTIMATE
+ b1_se[MC_rep,1]   = fit0$coef[1,2]
+ b1_se[MC_rep,2]   = fit1$coef[1,2]
+ b1_se[MC_rep,3]   = fit2$coef[1,2]
+ b1_se[MC_rep,4]   = fit3$coef[1,2]
+ b1_se[MC_rep,5]   = fit4$coef[1,2]
+ b1_se[MC_rep,6]   = fit5$coef[1,2]
+
+ b2_se[MC_rep,1]   = fit0$coef[2,2]
+ b2_se[MC_rep,2]   = fit1$coef[2,2]
+ b2_se[MC_rep,3]   = fit2$coef[2,2]
+ b2_se[MC_rep,4]   = fit3$coef[2,2]
+ b2_se[MC_rep,5]   = fit4$coef[2,2]
+ b2_se[MC_rep,6]   = fit5$coef[2,2]
+
+ b3_se[MC_rep,1]   = fit0$coef[3,2]
+ b3_se[MC_rep,2]   = fit1$coef[3,2]
+ b3_se[MC_rep,3]   = fit2$coef[3,2]
+ b3_se[MC_rep,4]   = fit3$coef[3,2]
+ b3_se[MC_rep,5]   = fit4$coef[3,2]
+ b3_se[MC_rep,6]   = fit5$coef[3,2]
+
+ }
>
> ### COMBINE THE RESULTS INTO AN EASY-TO-READ TABLE
> results.mat<-rbind(
+ round( abs( colMeans( b1_bias ) ) ,3),
+ round( abs( colMeans( b2_bias ) ) ,3),
+ round( abs( colMeans( b3_bias ) ) ,3),
+ round( abs( colMeans( b1_se ) ) ,3),
+ round( abs( colMeans( b2_se ) ) ,3),
+ round( abs( colMeans( b3_se ) ) ,3),
+
+ )
Error in rbind(round(abs(colMeans(b1_bias)), 3), round(abs(colMeans(b2_bias)),  :
  缺少参数, 也没有缺省值
>
> colnames(results.mat)<-c("n=5","n=10","n=50","n=100","n=200","n=1000")
> rownames(results.mat)<-c("b1.Bias","b2.Bias","b3.Bias","b1.SE","b2.SE","b3.SE")
Error in `rownames<-`(`*tmp*`, value = c("b1.Bias", "b2.Bias", "b3.Bias",  :
  'dimnames'的长度[1]必需与陈列范围相等
>
> ### NOW LOOK AT THE RESULTS
> print(results.mat)
          n=5  n=10  n=50 n=100 n=200 n=1000
b1.Bias 0.052 0.057 0.034 0.013 0.004  0.002
b2.Bias 0.017 0.005 0.002 0.000 0.000  0.000
b1.SE   1.708 1.158 0.496 0.349 0.246  0.110
b1.SE   0.515 0.187 0.017 0.006 0.002  0.000
>

一个MC模拟写成这样子，楼主也是辛苦了。 建议写成函数， 根据不同次数输出结果，再合并比较。 用apply函数吧，数据量比较大的时候尽量少用rbind，cbind。
round函数有两个参数~