I've written the following code to compare the theoretical alpha = 0.05 with the empirical one from the buit-in t.test in Rstudio:
set.seed(1)
N <- 1000
n <- 20
k <- 500
poblacion <- rnorm(N, 10, 10) #Sample
mu.pob <- mean(poblacion)
sd.pob <- sd(poblacion)
p <- vector(length=k)
for (i in 1:k) {
muestra <- poblacion[sample(1:N, n)]
p[i] <- t.test(muestra, mu=mu.pob)$p.value
}
a_teo <- 0.05
a_emp <- length(p[p < a_teo])/k
sprintf("alpha_teo = %.3f <-> alpha_emp = %.3f", a_teo, a_emp)
And it works printing both theoretical and empirical values. Now I wanna make it more general, to different values of 'n', so I wrote this:
set.seed(1)
N <- 1000
n <- 20
k <- 500
z <-c()
for (i in n){
poblacion <- rnorm(N, 10, 10)
mu.pob <- mean(poblacion)
sd.pob <- sd(poblacion)
p <- vector(length=k)
for (j in 1:k){
muestra <- poblacion[sample(1:N, length(n))]
p[j] <- t.test(muestra, mu = mu.pob)$p.value
}
a_teo = 0.05
a_emp = length(p[p<a_teo])/k
append(z, a_emp)
print(sprintf("alpha_teo = %.3f <-> alpha_emp = %.3f", a_teo, a_emp))
}
plot(a_teo, z)
Finally, I want to plot in one single graph the empirical (a_emp) values against the 0.05. I've made a 'z' vector to store the empirical values, but I get only one point in the graph. Any ideas?