机器算法验证 - 球是随机抽取的吗（与其颜色中存在的球数量无关）？ - 吾爱随笔录

球是随机抽取的吗（与其颜色中存在的球数量无关）？

机器算法验证采样斯皮尔曼罗超几何分布

2022-04-08 08:55:09

我们有一个大骨灰盒，里面装着 $N_{Tot}$ 球。球是 $r$ 不同的颜色。球数 $i^{th}$ 颜色（采样前）是 $N_i$ . 约翰抽样 $x$ 这个瓮中的球总数（无需更换）。每种颜色的球数是 $x_1, x_2, ..., x_{r-1}, x_r$ .

我试图回答的问题：约翰是否独立于球的颜色中存在的球数量对球进行采样？

我的想法

我首先想到了一个简单的回归 $r$ 数据点，其中x变量是该颜色的球y数，坐标是该颜色采样的球数。但这行不通。一个明显的原因是每个点都不是一个独立的观察，因为从一种颜色中采样第一个球意味着不以另一种颜色采样第一个球。

$x_i$ 应该遵循带参数的超几何分布 $N_{tot}-N_i$ （失败）， $N_i$ （成功）和 $x$ （抽出的球总数），其中 $N_{tot} = \sum_{i=1}^r N_i$ 是球的总数。我可以比较 $x_i$ 到超几何分布的中位数，以判断每种颜色的瓮在样本中的比例过高还是过低。我们称这个二进制变量 $B$

然后，我可以执行 Spearman 相关性 $B$ 以及每种颜色的球数 $N_1, N_2, ..., N_{r-1}, N_r$ . 如果 Spearman 相关性的 p 值足够低，我可以得出结论，这些球不是随机抽取的，与球的颜色无关。

这似乎是一个很好的方法？

这会是低功率测试还是完全错误？

1个回答

因为您提出的测试没有理论基础，也没有考虑计数之间的相关性，所以我们不能很好地利用我们的时间来评估它。相反，让我们开发一些可行的测试。

当所有计数都相对较大时，可以保证似然比检验运行良好。卡方检验也是如此，只要样品中只有一小部分瓮中的球被取出。（通常的经验法则是当样本超过总数的 10% 时要小心。）在其他情况下，从零分布进行模拟是有效的。

观察计数的可能性 $x_1,\ldots, x_r$ 在一个瓮中的样本中 $N_1,\ldots, N_r$ 每种颜色的球是

L (x; N) = (\binom{N_{1}}{x_{1}}) (\binom{N_{2}}{x_{2}}) \dots (\binom{N_{r}}{x_{r}}) / (\binom{N_{1} + \dots + N_{r}}{x_{1} + \dots + N_{r}}) .

$\mathcal{L}(x;N) = \binom{N_1}{x_1}\binom{N_2}{x_2}\cdots \binom{N_r}{x_r}\ /\ \binom{N_1+\cdots+N_r}{x_1+\cdots+N_r}.$

这最好用它的偏差来表达 $D=-2\log\mathcal L$ 因为渐近地，偏差有一个 $\chi^2$ 分布与 $r-1$ 自由程度

通过采样几千次（使用计算机），您可以估计 $D.$ 观测值的 p 值是由观测偏差确定的右尾面积。

例如，考虑一个带有 $r=7$ 数量上的颜色 $(N_i) = (34,45,41,35,49,47,51,42).$ 我计算了一万个大小样本的零偏差 $167$ （每个样本中有一半的球），如左图所示。我还计算了另外一万个样本中的偏差，其中颜色的选择概率从 20%（对于颜色 $1$ ) 降至仅 6%（用于颜色 $7$ ）。后一种偏差往往很大。实际上，其中大约 96% 超过了零偏差的第 95 个百分位。从这个意义上说，这个测试的力量，当在 $alpha=100-95\% = 5\%$ 水平，为 96%。

我们可以和平常一样玩同样的游戏 $\chi^2$ 测试。如果该检验直接适用，则其 p 值的零分布将是均匀的。但是，零分布是偏斜的，大多数 p 值都非常大：请参见下图中的左图。但是，在检查这个模拟后，我们可以声明卡方检验结果在 $\alpha$ 每当其 p 值小于 $\alpha$ 此零分布的分位数。该分位数接近 $0.4$ （而不是预期的 $0.05$ )，用阴影条显示。

右图同样显示了备择假设下卡方检验 p 值的分布。它的“标称功效”是报告的（标称）p 值小于 $\alpha,$ 显示为红色条的区域。这只有 48.2%，远低于似然比 (LR) 测试所达到的 96%。然而，当我们使用 $\alpha$ 零分布的分位数作为我们的阈值，现在零假设在 95% 的模拟数据中被拒绝。95% 的功效与 LR 测试基本相同。

那么，这些模拟所证明的是

应用标准卡方检验或 LR 检验是无效的，除非样本是骨灰盒的一小部分并且样本数量相当大。第二个图中的偏态空分布显示出了什么问题。
然而，如果 p 值是使用实际的零分布（通过模拟估计的）而不是使用标准公式（依赖于具有大量人口的渐近大样本）计算的，则可以使用这两种测试。

R执行这些计算并显示这些数字的代码如下。

#
# Multivariate hypergeometric distribution.
# Draw `n` balls without replacement from an urn with length(N) colors, each
# appearing N[i] times.
#
# Returns a vector of counts of the colors.
#
rmhyper <- function(n, N, p) {
  if (missing(p)) p <- rep(1/sum(N), length(N))
  prob <- rep(p, N)
  prob <- prob / sum(prob)
  x <- sample(rep(seq_along(N), N), n, prob=prob)
  tabulate(x, length(N))
}
#
# Returns the likelihood of observing `k` in a hypergeometric draw specified
# by `N`.
#
dmhyper <- function(k, N, log.p=TRUE) {
  q <- sum(lchoose(N, k)) - lchoose(sum(N), sum(k))
  if (!isTRUE(log.p)) q <- exp(q)
  return(q)
}
#
# Perform a chi-squared test.
#
mhyper.test <- function(k, N, ...) {
  chisq.test(k, p=N / sum(N), ...)
}
#
# Simulations.
#
alpha <- 0.05
n.sim <- 1e4
set.seed(17)
N <- rpois(8, 40) + 1          # Determine the population randomly
p <- rgamma(length(N), 10/4)   # Determine the alternative hypothesis randomly
p <- rev(sort(p / sum(p)))
n <- ceiling(sum(N)*0.5)
#
# Display the sampling probabilities.
#
plot(p, type="h", col=rainbow(length(N), .9, .8),
     lwd=3, ylim=c(0, max(p)), xlab="Color", main="Probabilities")
#
# The chi-squared test p-values are ok when the sample is a small fraction
# of the population.
#
p.values.null <- replicate(n.sim, mhyper.test(rmhyper(n,N), N, simulate.p.value=FALSE)$p.value)
p.values.alt <- replicate(n.sim, mhyper.test(rmhyper(n,N, p), N, simulate.p.value=FALSE)$p.value)
power <- round(100*mean(p.values.alt <= alpha), 1)
power.alt <- round(100*mean(p.values.alt <= quantile(p.values.null, alpha)))
k <- round(20 * quantile(p.values.null, alpha)) - 1

b <- seq(0, 1, by=0.05)
par(mfrow=c(1,2))
h <- hist(p.values.null, freq=FALSE, breaks=b, xlab="p",
     col=c("#d08080", rep("#e0e0e0", k), rep("White", 20-k-1)),
     main="Histogram of Chi-square P-values\nNull Distribution")
hist(p.values.alt, freq=FALSE, breaks=b, xlab="p", ylim=c(0, max(h$density)),
     col=c("#d08080", rep("#e0e0e0", k), rep("White", 20-k-1)),
     sub=bquote(paste("Nominal power = ", .(power), "%; Simulation power = ", .(power.alt), "%")),
     main="Histogram of Chi-square P-values\nAlternative Distribution")

q.null <- -2 * replicate(n.sim, dmhyper(rmhyper(n, N), N))
q <- -2 * replicate(n.sim, dmhyper(rmhyper(n, N, p), N))


xlim <- range(c(q, q.null))

h <- hist(q.null, xlim=xlim, freq=FALSE, breaks=50, xlab=expression(-2~~log(p)),
     col="#f0f0f0",
     main="Null distribution of deviance")
abline(v = quantile(q.null, 1-alpha), col="Red", lwd=2)

power <- signif(mean(q >= quantile(q.null, 1-0.05)), 2)
hist(q, xlim=xlim, freq=FALSE, breaks=50, xlab=expression(-2~~log(p)),
     col="#f0f0f0", ylim=c(0, max(h$density)),
     main="Alternative distribution of deviance",
     sub=bquote(paste("Power = ", .(100*power), "%")))
abline(v = quantile(q.null, 1-alpha), col="Red", lwd=2)
par(mfrow=c(1,1))

其它你可能感兴趣的问题

上一篇广义加性模型 (GAM) 的限制下一篇关于生命表中有效样本量的问题