我编写了以下代码来将理论 alpha = 0.05 与 Rstudio 中内置 t.test 的经验值进行比较:
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)
它可以打印理论值和经验值。现在我想让它更通用,针对不同的“n”值,所以我写了这个:
set.seed(1)
N <- 1000
n <- c(2, 10, 15, 20)
k <- 500
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
sprintf("alpha_teo = %.3f <-> alpha_emp = %.3f", a_teo, a_emp)
}
但我没有得到“打印”。关于哪里出了问题有什么想法吗?
单独使用
sprintf
无法形成 for
循环,您需要将其包裹在 print
中。
> 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
+ print(sprintf("alpha_teo = %.3f <-> alpha_emp = %.3f", a_teo, a_emp))
+ }
[1] "alpha_teo = 0.050 <-> alpha_emp = 0.056"
[1] "alpha_teo = 0.050 <-> alpha_emp = 0.050"
[1] "alpha_teo = 0.050 <-> alpha_emp = 0.064"
[1] "alpha_teo = 0.050 <-> alpha_emp = 0.048"
一种更 R 风格的方法是将逻辑包装在函数中。
> comp_fn <- \(N, n, k, alpha=.05, verbose=FALSE) {
+ poblacion <- rnorm(N, 10, 10)
+ mu.pob <- mean(poblacion)
+ sd.pob <- sd(poblacion)
+ p <- replicate(k, t.test(poblacion[sample(1:N, n)], mu=mu.pob)$p.value)
+ a_emp <- length(p[p < alpha])/k
+ if (verbose) {
+ message(sprintf("alpha_teo = %.3f <-> alpha_emp = %.3f", a_teo, a_emp))
+ }
+ c(a_teo, a_emp)
+ }
>
> set.seed(1)
> comp_fn(1000, 20, 500)
[1] 0.050 0.058
> comp_fn(1000, 20, 500, verbose=TRUE)
alpha_teo = 0.050 <-> alpha_emp = 0.042
[1] 0.050 0.042
要循环不同的参数,
mapply
是你的朋友。
> set.seed(1)
> mapply(comp_fn, 1000, c(2, 10, 15, 20), 500)
[,1] [,2] [,3] [,4]
[1,] 0.050 0.050 0.050 0.050
[2,] 0.058 0.054 0.048 0.046