rs = matrix(0,10,100)
xx = rnorm(10*10)
dim(xx) = c(10,10)
colnames(xx) = paste('x', 1:10, sep='')
ff = 2*xx[,1]+2*xx[,2]
yy = ff + rnorm(10, 0, 1)
for (k in 1:100)
{
x = rnorm(100*10)
dim(x) = c(100,10)
colnames(x) = paste('x', 1:10, sep='')
y = 2*x[,1]+2*x[,2] + rnorm(100, 0, 1)
mm = coef(lm(y~as.matrix(x)))
pyy = mm[1] + xx %*% mm[-1]
rs[,k] = pyy[]
}
rss = bais2 = var = 0
for ( k in 1:10)
{
rss = rss+mean((rs[,k] - yy[k])^2)/10
bias2 = (mean(rs[,k]) - ff[k])^2/10
var = var + var(rs[,k])/10
}rss
10.4668988127494
bias2
0.115522433888335
var
5.86777707464027