2
我試圖繪製的卡方test.In R中的零分佈是可行的做蒙特卡羅模擬使用代碼來獲得經驗p值:檢索的卡方檢驗蒙特卡羅模擬值
chisq.test(d,simulate.p.value=TRUE,B=10000)
但它不返回分佈圖。有沒有辦法讓R返回測試的模擬值?
我試圖繪製的卡方test.In R中的零分佈是可行的做蒙特卡羅模擬使用代碼來獲得經驗p值:檢索的卡方檢驗蒙特卡羅模擬值
chisq.test(d,simulate.p.value=TRUE,B=10000)
但它不返回分佈圖。有沒有辦法讓R返回測試的模擬值?
如果你看一下函數定義chisq.test
(左右線capture.output(chisq.test)
56)你會來模擬部分:
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
包含所需的輸出。的
chisq.test(x, simulate.p.value = T, B=2000)$p.value
注輸出我使用:::
變量PVAL
比賽,因爲函數C_chisq_sim
不是從統計出口。