Извлечение методом Монте-Карло значения моделирования для хи-квадрат тест

Question

Я пытаюсь построить график распределения нуль в хи-квадрат test.In R это возможно сделать метод Монте-Карло, чтобы получить эмпирическое значение р, используя код:

chisq.test(d,simulate.p.value=TRUE,B=10000) 

Но он не возвращает граф распределения. Есть ли способ заставить R вернуть имитированные значения для теста?

Answer 1

Если посмотреть на определение функции для chisq.test (вокруг линии 56 из capture.output(chisq.test)), вы попадете в раздел моделирования:

if (simulate.p.value && all(sr > 0) && all(sc > 0)) { 
     setMETH() 
     tmp <- .Call(C_chisq_sim, sr, sc, B, E) 
     STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE)) 
     PARAMETER <- NA 
     PVAL <- (1 + sum(tmp >= almost.1 * STATISTIC))/(B + 
      1) 
    } 

Это вызывает функцию C. Сначала приготовьте несколько фиктивных данных

## Some data 
x <- as.table(rbind(c(762, 327, 468), c(484, 239, 477))) 
dimnames(x) <- list(gender = c("F", "M"), 
       party = c("Democrat","Independent", "Republican")) 

Затем возьмите немного, что вам нужно

sr <- rowSums(x) 
sc <- colSums(x) 
n <- sum(x) 
E <- outer(sr, sc, "*")/n 
v <- function(r, c, n) c * r * (n - r) * (n - c)/n^3 
V <- outer(sr, sc, v, n) 
dimnames(E) <- dimnames(x) 
B = 2000 
tmp <- .Call(stats:::C_chisq_sim, sr, sc, B, E) 
STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE)) 
almost.1 <- 1 - 64 * .Machine$double.eps             
PVAL <- (1 + sum(tmp >= almost.1 * STATISTIC))/(B + 1) 

Переменная tmp содержит вывод, который вы хотите. Переменная PVAL согласует выход

chisq.test(x, simulate.p.value = T, B=2000)$p.value 

---

Примечание Я использовал :::, так как функция C_chisq_sim не экспортируется из статистики.