Как итеративно вычислить p-значения для t-теста - PullRequest
1 голос
/ 18 апреля 2019

а) Создайте 50 значений из X ~ N (μX = 25, σX = 4) и 50 значений из Y ~ N (μY = 25, σY = 4). Используйте t-критерий для проверки на равенство средних.

в) Повторите часть (а) 2500 раз и сохраните значение p для каждого из 2500 тестов. Каждое повторение должно генерировать новый образец для x и новый образец для y. НЕ ПЕЧАТЬ р-значения. НЕ используйте петлю.

Я решил для Части A на одной rnorm выборке, но я не понимаю, с чего начать, чтобы получить 2500 различных случайных выборок x и 2500 различных случайных выборок y, чтобы получить 2500 различных значений p.

Я также не знаю, как написать свой код, чтобы мой профессор получал те же ответы, что и я. Я попытался установить начальное значение, но это только делает так, чтобы значения p были одинаковыми, используя мой код выше.

# Part A

set.seed(1081)
x = rnorm(50,25,4)
y = rnorm(50,25,4)

t.test(x,y)

#Part B
#The p-value is 0.3752.
#We do not reject the null hypothesis.

#Part C

x1 = sample(x, 2500, replace = T)
y1 = sample(y, 2500, replace = T)
pval = sample(t.test(x1,y1)$p.value, 2500, replace = T)

Ответы [ 3 ]

1 голос
/ 19 апреля 2019

другая возможность заключается в использовании replicate:

Обратите внимание, что вы должны установить случайное начальное число вне функции.

myfun <- function(){
  x <- rnorm(50, 25, 4)
  y <- rnorm(50, 25, 4)

  return(t.test(x, y)$p.value)
}


set.seed(1)
p_vals <- replicate(2500, myfun())
1 голос
/ 19 апреля 2019

Другой подход заключается в следующем:

    library(MASS)       #load MASS library

    s <- 4*diag(2500)   #create the variance matrix for the simulation
    set.seed(123)        # seed to replicate results

    x <- mvrnorm( 50, m= rep(25,times=2500), Sigma=s)  #draw 50 values, 25000 times 

    y <- mvrnorm( 50, m = rep(25, times=2500), Sigma=s) #draw 50 values, 2500 times

    diff <- x - y

    test <- apply(diff,2,t.test) #do the t.tests

    names(test) #some of the results you can print

Если у вас есть вопросы по поводу кода, вы можете задать их мне.

0 голосов
/ 19 апреля 2019

Еще одна возможность:

set.seed(1081)
n <- 50
times <- 2500
x <- data.frame(matrix(rnorm(n*times, mean=25, sd=4), nrow=n))
y <- data.frame(matrix(rnorm(n*times, mean=25, sd=4), nrow=n))
pvals <- mapply(FUN = function(x,y) t.test(x,y)$p.value, x, y)
mean(pvals < .05)  # should be ~= .05

Зацикливание одновременно над двумя списками в R (комментарий Джого)

Но если мы возьмем ", каждое повторение должно генерировать новоеОбразцы "буквально, ответ @ Cettt может быть то, что хотел.

...