在 R 脚本中确定 n 个样本大小以在 alpha = 0.01 处拒绝 null
Determining the n sample size to reject null at alpha =0.01 in R script
创建一个包含 100 个观测值的模拟数据集,其中 x 是均值为 0 且标准差为 1 的随机正态变量,y = 0.1 + 2 * X + e,其中 epsilon 也是均值为 0 的随机正态误差和标准差 1.
set.seed(1)
# simulate a data set of 100 observations
x <- rnorm(100)
y.1 <- 0.1 + 2*x + rnorm(100)
现在提取前 5 个观测值。
y1.FirstFive <- (y.1[1:5]) # extract first 5 observations from y
x.FirstFive <- (x[1:5]) # extract first 5 observations from x
y1.FirstFive # extracted 5 observations from y1
[1] -1.7732743 0.5094025 -2.4821789 3.4485904 0.1044309
x.FirstFive # extracted 5 observations from x
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
假设您根据前五个观察值计算出的样本的均值和标准差不会改变,那么您需要额外观察的最少总数是多少才能得出总体的真实均值是在 p = 0.01 置信水平下与 0 不同?
alpha <- 0.01
mu <- 0
for (i in 5:2000) {
# Recalculate the standard error and CI
stand_err <- Sd_y1 / sqrt(i)
ci <- sample_mean_y1 + c(qt(alpha/2, i-1), qt(1-alpha/2, i-1))*stand_err
if (ci[2] < mu)
break # condition met, exit loop
}
i
[1] 2000
在这里,我写了一个循环,将 n 从最初的 n=5 迭代增加到 n=2000,使用 pt 找到 p 值(给定固定的 y-bar 和 sd),并在 p < 0.01 时停止.但是我一直得到错误的输出。这样,输出始终是我给出的最大范围的数字(这里是 2000),而不是给我特定的最小 n 样本,以便拒绝 p= 处 mu_y = 0 的空值0.01 级。关于如何修复代码有什么建议吗?
- 附加信息:y1.FirstFive 的标准差 = 2.3 和 y1.FirstFive 的平均值 = -0.04
假设:
Sd_y1 = sd(y1.FirstFive)
sample_mean_y1 = mean(y1.FirstFive)
sample_mean_y1
[1] -0.03860587
正如@jblood94 所指出的,您需要更大的样本量。
你不需要 for 循环,你的大部分函数都是矢量化的,所以像这样:
n = 5:30000
stand_err = Sd_y1 / sqrt(n)
ub = sample_mean_y1 + qt(1-alpha/2, n-1)*stand_err
n[min(which(ub<0))]
[1] 23889
因为n>2000
set.seed(1)
x <- rnorm(100)
y.1 <- 0.1 + 2*x + rnorm(100)
Sd_y1 <- sd(y.1[1:5])
sample_mean_y1 <- mean(y.1[1:5])
alpha <- 0.01
sgn <- 2*(sample_mean_y1 > 0) - 1
f <- function(n) qt(alpha/2, n - 1)*Sd_y1 + sgn*sample_mean_y1*sqrt(n)
upper <- 2
while (f(upper) < 0) upper <- upper*2
(n <- ceiling(uniroot(f, lower = upper/2, upper = upper, tol = 0.5)$root))
#> [1] 23889
创建一个包含 100 个观测值的模拟数据集,其中 x 是均值为 0 且标准差为 1 的随机正态变量,y = 0.1 + 2 * X + e,其中 epsilon 也是均值为 0 的随机正态误差和标准差 1.
set.seed(1)
# simulate a data set of 100 observations
x <- rnorm(100)
y.1 <- 0.1 + 2*x + rnorm(100)
现在提取前 5 个观测值。
y1.FirstFive <- (y.1[1:5]) # extract first 5 observations from y
x.FirstFive <- (x[1:5]) # extract first 5 observations from x
y1.FirstFive # extracted 5 observations from y1
[1] -1.7732743 0.5094025 -2.4821789 3.4485904 0.1044309
x.FirstFive # extracted 5 observations from x
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
假设您根据前五个观察值计算出的样本的均值和标准差不会改变,那么您需要额外观察的最少总数是多少才能得出总体的真实均值是在 p = 0.01 置信水平下与 0 不同?
alpha <- 0.01
mu <- 0
for (i in 5:2000) {
# Recalculate the standard error and CI
stand_err <- Sd_y1 / sqrt(i)
ci <- sample_mean_y1 + c(qt(alpha/2, i-1), qt(1-alpha/2, i-1))*stand_err
if (ci[2] < mu)
break # condition met, exit loop
}
i
[1] 2000
在这里,我写了一个循环,将 n 从最初的 n=5 迭代增加到 n=2000,使用 pt 找到 p 值(给定固定的 y-bar 和 sd),并在 p < 0.01 时停止.但是我一直得到错误的输出。这样,输出始终是我给出的最大范围的数字(这里是 2000),而不是给我特定的最小 n 样本,以便拒绝 p= 处 mu_y = 0 的空值0.01 级。关于如何修复代码有什么建议吗?
- 附加信息:y1.FirstFive 的标准差 = 2.3 和 y1.FirstFive 的平均值 = -0.04
假设:
Sd_y1 = sd(y1.FirstFive)
sample_mean_y1 = mean(y1.FirstFive)
sample_mean_y1
[1] -0.03860587
正如@jblood94 所指出的,您需要更大的样本量。
你不需要 for 循环,你的大部分函数都是矢量化的,所以像这样:
n = 5:30000
stand_err = Sd_y1 / sqrt(n)
ub = sample_mean_y1 + qt(1-alpha/2, n-1)*stand_err
n[min(which(ub<0))]
[1] 23889
因为n>2000
set.seed(1)
x <- rnorm(100)
y.1 <- 0.1 + 2*x + rnorm(100)
Sd_y1 <- sd(y.1[1:5])
sample_mean_y1 <- mean(y.1[1:5])
alpha <- 0.01
sgn <- 2*(sample_mean_y1 > 0) - 1
f <- function(n) qt(alpha/2, n - 1)*Sd_y1 + sgn*sample_mean_y1*sqrt(n)
upper <- 2
while (f(upper) < 0) upper <- upper*2
(n <- ceiling(uniroot(f, lower = upper/2, upper = upper, tol = 0.5)$root))
#> [1] 23889