sampleSize <- 10;

numIter <- 1000;
ubVar <- rep(0,numIter)

for (i in 1:numIter) {
    x <- rnorm(sampleSize);
    ubVar[i] <- var(x)
}

bVar <- ubVar * ((sampleSize-1)/sampleSize)

var <- c(ubVar,bVar);
label <- as.factor(c(rep("sigmaPrime",numIter),rep("sigmaHat",numIter)))
boxplot(var~label)
abline(1,0,col=1)
points(1,mean(bVar),pch=16,col=1)
points(2,mean(ubVar),pch=16,col=2)