ran.sd <- c(20,5)
ran.var <- ran.sd^2
ran.size <- c(10,10)
ran.mean <- 100
ran <- data.frame(matrix(rnorm(prod(ran.size),0,ran.sd[2]),
   ran.size[2],ran.size[1]))
names(ran) <- 1:(ran.size[1])
for (i in names(ran))
   ran[[i]] <- ran[[i]] + rnorm(1,ran.mean,ran.sd[1])
ran.df <- c(ran.size[1]-1,ran.size[1]*(ran.size[2]-1))

ran.mean <- mean(unlist(ran))
ran.class <- apply(ran,2,mean)
ran.element <- apply(ran,2,function(x){x-mean(x)})
ran.var <- (ran.df[2]*ran.var[2] + 
   ran.df[1]*(ran.var[2]+ran.size[2]*ran.var[1])) / sum(ran.df)