这里有三个问题:
模拟可能并不完全像看起来那样。
该分析不评估均值和方差之间的关系是否为线性关系。
如何检查数据是否出现泊松?
#3 的许多答案已经出现在我们网站的其他地方,因此我将重点关注前两个问题。
模拟的作用
代码行
pos<-rpois(100,5)
...
samp<-sample(pos, 10)
从均值为的泊松分布中抽取iid 值,然后随机获得个值(无需替换)。迭代第二个命令是重采样分析的一种形式。它并没有真正测试当重复从相同的基础泊松分布中提取数据时会发生什么;相反,它探索了在第一行中获得值的分布,这些值在模拟期间是固定的。100510100
为了研究泊松分布,我们需要从中抽取重复的独立样本:
set.seed(17) # Create reproducible results
mu <- 3 # Mean
n <- 10 # Sample size
n.trials <- 10000 # Number of trials
sim <- replicate(n.trials, {x <- rpois(n, mu); c(mean(x), var(x))})
此时sim是一个包含两行的矩阵:一个用于均值,下一个用于方差。每列总结一个单独的模拟(ergo,本例中有列)。10000
分析的作用
这是一个紧凑的形式:
xy <- apply(sim, 1, function(x) quantile(x, probs=seq(0.001, 0.999, length.out=99)))
这将创建两行输出的百分位数数组,从而追踪它们的分布。一个情节很有启发性:
plot(xy, xlab="mean", ylab="variance",col=hsv(0, 0, 0.25, .4))

与任何 QQ 图一样,该图系统地将一个分布的百分位数(模拟泊松值样本的方差)与另一种分布的百分位数(模拟泊松值样本的平均值)进行比较。1010
曲线并不是真的笔直,我们可以通过使用包含二次项的 x 值回归 y 值来发现:
z <- xy[,1]^2 # Quadratic term
fit <- lm(xy[,2] ~ xy[,1] + z) # Regression of y-values against x and x^2
#
# Overplot the fit in red.
x0 <- seq(mu - 2 * sqrt(mu/n), mu + 2 * sqrt(mu/n), length.out=99)
u <- cbind(rep(1,length(x0)), x0, x0^2) %*% coefficients(fit)
lines(x0, u, lwd=2, col=hsv(0, .8, .8, 1))

二次项非常显着并且具有很大的系数:
>coefficients(fit)[3]
z
0.8779665
我鼓励感兴趣的读者使用参数(mu、n和n.trials)来看看事情(不要)如何变化。还可以考虑rpois将模拟更改为其他内容(例如rnorm)。
一些理论
次独立抽取的平均值将近似为正态分布。泊松(均值)也不例外:其均值应具有均值和方差的正态分布。10555/10
具有相当大均值(例如,或更大)的泊松分布的 CDF 本身将开始逼近正态分布的 CDF,并且对于更大均值的逼近会有所改善。因此,我们可以希望次独立泊松抽取样本的方差将开始假设次独立抽取样本的方差分布来自可比较的正态分布。(方差具有与卡方分布成比例的分布。)我们可以通过计算这些近似分布的确切百分位数并将它们重叠在模拟值上进行检查:31010
q <- (1:100 - 0.5) / 100 # An array of percentiles for plotting
lines(qnorm(q, mean=mu, sd=sqrt(mu/n)), mu*qchisq(q, df=n-1) / (n-1), col="Blue")

即使对于这个相当小的平均值 ( ) 和适度的样本量 ( ),拟合也非常好。不幸的是,这表明对于许多不同的基础分布,我们应该期待几乎任何中等样本量的相似图,从而表明该图很少或没有告诉我们关于数据是否来自泊松分布的有用信息。510
我们能做些什么呢?
通过指出数据可能偏离该形状的方式,可以改进“检查”分布是否具有特定形状。
例如,泊松分布是负二项分布的特殊例子,可以认为是过度分散的泊松分布。具体来说,如果我们获得一系列泊松值,其中泊松分布的基本均值以近似Gamma分布的方式变化,则收集的数据集由负二项式描述。
在这样的设置中,我们可以拟合负二项式参数(例如,使用最大似然)并测试“非泊松”参数是否与对应于泊松分布的值显着不同。出于篇幅的考虑,并且为了使这个回答保持合理的基本,我不会在这里寻求这个选项。
不太正式,我们可以简单地应用泊松广义线性模型并查看其诊断:
set.seed(17)
y <- rnbinom(100, mu=5, size=1) # NOT Poisson
x <- rpois(100, 5) # Poisson
fit.y <- glm(y ~ 1, family=poisson())
fit.x <- glm(x ~ 1, family=poisson())
plot(fit.y) # Includes a residual Q-Q plot
plot(fit.x) # Includes a residual Q-Q plot
Poisson 数据的残差 QQ 图几乎是线性的,而负二项式数据的残差 QQ 图有一些明显的极端


可以看出,检测这些差异需要适度的样本量,但至少它们确实出现了。