比较理论和经验 alpha

Question

我编写了以下代码来将理论 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 <- 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(n, z)

Answer 1

单独使用

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