В данном задании необходимо проанализировать поведение схожих критериев — t- и z- критерии (обновыборочные) — для задачи проверки гипотезы о равенстве нулю матожидания выборки. Для каждого из критериев требуется построить графики зависимости достигаемых уровней значимости и оценок мощностей от параметров и показать, в каких областях изменения параметров предпочтительнее использовать тот или иной критерий.

Исходные параметры эксперимента выглядят следующим образом:

\[ X^{n},\ X \sim N(\mu, \sigma);\\ H_0: среднее\ X\ равно\ нулю,\\ H_1: среднее\ X\ не\ равно\ нулю:\\ \mu = 0:0.01:2,\ \sigma = 1,\ n = 5:1:100 \] Инициализируем параметры для экспериментов

nsample <- seq(5, 100, by=1)
mu     <- seq(0,   2, by=0.01)
sd  <- 1
N       <- length(nsample)
M       <- length(mu)
grid    <- expand.grid(n=nsample, mu=mu)

Опишем функцию для вычисления p-value z-критерия:

z.test <- function(x, y = NULL, alterative = "two.sided", mu = 0, sd = 1, paired = FALSE) {
  result <- list()
  if (is.null(y)) {
    z = (mean(x) - mu)/(sd/sqrt(length(x)))
    if (alterative == "two.sided") {
      result$p.value <- 2*(1 - pnorm(abs(z), 0, 1))
    } else if (alterative == "greater") {
      result$p.value <- 1 - pnorm(z, 0, 1)
    } else if (alterative == "less") {
      result$p.value <- pnorm(z, 0, 1)
    }
    result$z <- z
  }
  return(result)
}

Проведем одиночный эксперимент с вычислением нобходимых значений

N_exps  <- 1
PV_Z    <- rep(0, N * M)
PV_T   <- rep(0, N * M)
Pow_Z   <- rep(0, N * M)
Pow_T  <- rep(0, N * M)

for (iter in 1:N_exps) {
    X1 <- apply(grid, 1, function(data) rnorm(data['n'], data['mu'], sd = 1))
    # splitList <- split(X1, nsample)
    
    TMP   <- sapply(X1, function(X) z.test(X, sd=sd)$p.value)
    PV_Z  <- PV_Z + TMP
    
    TMP    <- sapply(X1, function(X) t.test(X)$p.value)
    PV_T  <- PV_T + TMP
}
PV_Z   <- matrix(PV_Z   / N_exps, nrow=N, ncol=M)
PV_T  <- matrix(PV_T  / N_exps, nrow=N, ncol=M)

Рассмотрим полученные показатели уровня значимости и мощности

#library(fields)

par(mfrow=c(1,2))
image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), PV_Z, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024), 
           main="Z Test p-values", xlab=expression(n), ylab=expression(mu))

image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), PV_T, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024), 
           main="T Test p-values",xlab=expression(n), ylab=expression(mu))

Однократная генерация пары выборок не позволяет точно оценить границы области, где нулевая гипотеза отклоняется, поэтому необходимо усреднение по большому числу экспериментов. Возьмём 1000 повторений:

N_exps  <- 1000
PV_Z    <- rep(0, N * M)
PV_T   <- rep(0, N * M)
Pow_Z   <- rep(0, N * M)
Pow_T  <- rep(0, N * M)

for (iter in 1:N_exps) {
    X1 <- apply(grid, 1, function(data) rnorm(data['n'], data['mu'], sd = 1))
    # splitList <- split(X1, nsample)
    
    TMP   <- sapply(X1, function(X) z.test(X, sd=sd)$p.value)
    PV_Z  <- PV_Z + TMP
    Pow_Z <- Pow_Z + (TMP <= 0.05)
    
    TMP    <- sapply(X1, function(X) t.test(X)$p.value)
    PV_T  <- PV_T + TMP
    Pow_T <- Pow_T + (TMP <= 0.05)
}
PV_Z   <- matrix(PV_Z   / N_exps, nrow=N, ncol=M)
PV_T  <- matrix(PV_T  / N_exps, nrow=N, ncol=M)
Pow_Z  <- matrix(Pow_Z  / N_exps, nrow=N, ncol=M)
Pow_T <- matrix(Pow_T / N_exps, nrow=N, ncol=M)

Посмотрим сначала на средние достигаемые уровни значимости и мощности критериев при \(n=50\), когда предположения обоих выполняются:

par(mfrow=c(1,1))
plot(mu, PV_Z[46,],   col="red", type="l", xlab=expression(mu), ylab="Average p-value", main="", ylim=c(0,1))
lines(mu, PV_T[46,], col="blue")
legend("topright", c("Z", "T"), lty=c(1,1), col=c("red", "blue"))
grid()

par(mfrow=c(1,1))
plot(mu, Pow_Z[1,],   col="red", type="l", xlab=expression(mu), ylab="Estimated power", main="", 
     ylim=c(0,1))
lines(mu, Pow_T[1,], col="blue")
legend("bottomright", c("Z criteria", "T criteria"), lty=c(1,1), col=c("red", "blue"))
grid()

Z-критерии равномерно мощнее, и разница в мощностях может достигать порядка 0.2. Оба критерия параметрические, но Z-критерий настроен на известную нам дисперсию, поэтому более точен и при удалении \(\mu\) от истинного значения нуля его оценка мощности возрастает быстрее, к тому же достигает значения единицы уже на разницы в 2 стандартных отклонения (\(\mu=2\)). При рассматриваемом объёме выборок n=50 мощность Z-критерия достигает 60%, если \(mu=1\), т.е. отстоит от нуля на одно стандартное отклоение.
Средние достигаемые уровни значимости и оценки мощности критериев при различных значениях параметров \(n\) и \(\mu\):

par(mfrow=c(1,2))
image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), PV_Z, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024), 
           main="Z Test p-values", xlab=expression(n), ylab=expression(mu))

image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), PV_T, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024), 
           main="T Test p-values", xlab=expression(n), ylab=expression(mu))

par(mfrow=c(1,2))
image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), Pow_Z, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024),
           main="Z Test power", xlab=expression(n), ylab=expression(mu))

image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), Pow_T, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024),
           main="T Test power", xlab=expression(n), ylab=expression(mu))

par(mfrow=c(1,2))
image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), PV_Z - PV_T, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024),
           main="Average p-value difference (Z - T)", xlab=expression(n), ylab=expression(mu))

image.plot(matrix(grid$n, nrow=N, ncol=M), matrix(grid$mu, nrow=N, ncol=M), Pow_Z-Pow_T, 
           col = colorRampPalette(c("blue", "cyan", "yellow", "red"))(1024),
           main="Power difference (Z - T)", xlab=expression(n), ylab=expression(mu))

Видно, что разница в критериях убывает при увеличении размеров выборки, как для мощности, так и для достигаемого уровня значимости. Имеет значение отклонение от нуля матожидания, которое более всего влияет на уровне 0.6. Но различия между критериями перестают быт значимыми как при близких к ную значениях \(\mu\), так и больших различиях. Чтобы оценить корректность критериев, посмотрим отдельно, как часто гипотеза о равенстве нулю матожидания выборки отвергается при \(\mu=0\) (то есть, на частоту ошибок первого рода) при различных \(n\).

N_exps <- 10000
T1_Z  <- rep(0, N)
T1_T  <- rep(0, N)

for (iter in 1:N_exps) {
    X1 <- lapply(nsample,function(n) rnorm(n, 0, sd = 1))
    
    TMP   <- sapply(X1, function(X) z.test(X, sd=sd)$p.value)
    T1_Z  <- T1_Z + (TMP <= 0.05)
    
    TMP    <- sapply(X1, function(X) t.test(X)$p.value)
    T1_T  <- T1_T + (TMP <= 0.05)
}
T1_Z  <- T1_Z  / N_exps
T1_T <- T1_T / N_exps
par(mfrow=c(1,1))
plot(nsample, T1_Z,   col="red", type="l", xlab=expression(n), ylab="Type I error frequency", main="", 
     ylim=c(min(T1_Z, T1_T), max(T1_Z, T1_T)))
lines(nsample, T1_T, col="blue")
legend("topright", c("Z", "T"), lty=c(1,1), col=c("red", "blue"))

Оба критерия корректны, так как при том, что верна нулевая гипотеза, ошибки первого рода мало колеблются вокруг вокруг уровня значимости \(\alpha = 0.05\).

Заключение

В данной модельной задаче все предполоения нормальности выполнены для обоих параметричеких Z-,T-критериев, но так как есть дополнительная информация о распределнии дисперсии выборки, она влияет на мощности Z-критерия, делая его мощнее и более надежным, кгда речь заходит о небольших отклонениях от нулевого матожидания, в 1 или 1.5 стандартных отклонения.