# function simsample
simsample = function(mean = 10, stdev = 1, maxsize = 20, reps = 10, ...){
index = 0
x = rep(1:maxsize, times = reps)
y = c(1:maxsize*reps)
for (i in 1:reps){
for (j in 1:maxsize){
index = index + 1
y[index] = mean(rnorm(j, mean, stdev))
}
}
plot(x, y, pch = 19, col = "blue", cex = 0.5, ylim = c(mean - 3 * stdev, mean + 3 * stdev), ylab = "mean", xlab = "sample size", ...)
abline(h = mean, col = "green", lwd = 2)
abline(h = mean + qnorm(0.25) * stdev, lty = 2, col = "red", lwd = 2)
abline(h = mean - qnorm(0.25) * stdev, lty = 2, col = "red", lwd = 2)
out = list("samplesize" = x, "average" = y)
invisible(out)
}