如何知道数据是否遵循 R 中的泊松分布?

机器算法验证 r 自习 泊松分布 泊松过程
2022-01-19 03:37:46

我是一名本科生,我的概率课有一个项目。基本上,我有一个关于影响我国家多年的飓风的数据集。

在我的概率书(Probability and Statistics with R)中有一个(不完整的)示例说明如何检查数据是否遵循泊松分布,他们开始尝试证明遵循这三个标准:(来自我的书,页面120(标准)第 122-123 页示例)

1- 非重叠间隔中的结果数量是独立的。换言之,时间间隔 (0,t] 内的结果数量与时间间隔 (t, t+h], h > 0) 内的结果数量无关

2- 在足够短的时间间隔内出现两个或多个结果的概率几乎为零。换句话说,如果 h 足够小,那么与在相同时间间隔内获得一个或零个结果的概率相比,在间隔 (t,t+h] 内获得两个或多个结果的概率可以忽略不计。

3- 在足够短的间隔或小区域中恰好出现一个结果的概率与间隔或区域的长度成正比。换句话说,在长度为 h 的区间中出现一个结果的概率是 lambda*h。

但是标准 3 被“作为练习”保留。

A-有人可以告诉我是否有更“简单”的方法来查看我的数据集是否遵循泊松分布?

B- 有人可以用某种类型的例子来解释我的标准 1 和 3(如果是 R 的话,太棒了)?

谢谢!

注意:对不起,很长的帖子。另外,我必须转换数据,以便我有一个像这样的表:

  number of hurricanes       | 0 | 1 | 2  etc.
  -----------------------------------------
  total years that have      |   |   |
  that number of hurricanes  |   |   |
2个回答

分布与泊松分布略有不同的方式有无数种。您无法确定一组数据从泊松分布中提取的。您可以做的是寻找与泊松应该看到的不一致,但缺乏明显的不一致并不能使其成为泊松。

但是,您通过检查这三个标准所谈论的不是通过统计方法(即通过查看数据)检查数据是否来自泊松分布,而是通过评估生成数据的过程是否满足泊松过程的条件;如果条件全部成立或几乎成立(这是对数据生成过程的考虑),您可能会得到来自或非常接近泊松过程的东西,这反过来又是一种从接近泊松过程的东西中获取数据的方式泊松分布。

但是条件在几个方面并不成立……离真实最远的是数字 3。在此基础上没有特别的理由来断言泊松过程,尽管违规可能不会那么严重,以至于结果数据相距甚远从泊松。

所以我们回到了来自检查数据本身的统计论点。

如开头所述,您可以做的是检查数据是否与泊松的基础分布明显不一致,但这并不能告诉您它们是从泊松中提取的(您已经可以确信它们是不是)。

您可以通过拟合优度测试进行此检查。

提到的卡方就是这样一种,但我自己不建议对这种情况进行卡方检验**;对于许多人认为最相关的偏差(例如,相对于具有相同均值的泊松的色散变化、偏度变化等),它的功效很低。如果你的目标是拥有强大的力量来对抗这些效果,你就不会那样做。

卡方的主要价值在于简单,具有教学价值;除此之外,它通常不像拟合优度测试那样具有竞争力。

** 在稍后的编辑中添加:现在很明显这是家庭作业,您需要进行卡方检验以检查数据与泊松不矛盾的机会上升了很多。请参阅我的示例卡方拟合优度检验,在第一个泊松图下方完成


人们经常出于错误的原因进行这些测试(例如,因为他们想说'因此可以对假设数据是泊松的数据做一些其他统计事情')。真正的问题是“这会发生多大的错误?” ...而拟合测试的优劣对这个问题并没有太大帮助。通常,该问题的答案充其量是独立于(/几乎独立于)样本量的——在某些情况下,其后果往往会随着样本量而消失……而拟合优度测试对小样本(违反假设的风险通常最大)。

如果您必须测试泊松分布,则有一些合理的选择。一种方法是根据 AD 统计量做类似于 Anderson-Darling 检验的事情,但使用零下的模拟分布(以解决离散分布的孪生问题并且您必须估计参数)。

一个更简单的替代方法可能是拟合优度的平滑测试 - 这些是通过使用一系列多项式对数据进行建模而为单个分布设计的测试集合,这些多项式与零值中的概率函数正交。低阶(即有趣的)备选方案通过测试基数之上的多项式的系数是否不同于零来测试,这些通常可以通过从测试中省略最低阶项来处理参数估计。泊松有这样的测试。如果您需要,我可以挖掘参考。

您也可以在泊松图中使用相关性(或者,更像是 Shapiro-Francia 检验,可能是) - 例如 vs(参见 Hoaglin,1980 年)——作为检验统计量。n(1r2)log(xk)+log(k!)k

这是该计算(和绘图)的示例,在 R 中完成:

y=rpois(100,5)
n=length(y)
(x=table(y))
y
 0  1  2  3  4  5  6  7  8  9 10 
 1  2  7 15 19 25 14  7  5  1  4 

k=as.numeric(names(x))
plot(k,c(log(x)+lfactorial(k)))

泊松图,显示接近一条直线

这是我建议可用于泊松拟合优度检验的统计数据:

n*(1-cor(k,log(x)+lfactorial(k))^2)
[1] 1.0599

当然,要计算 p 值,您还需要模拟空值下的检验统计量的分布(我还没有讨论如何处理值范围内的零计数)。这应该会产生一个相当强大的测试。还有许多其他替代测试。

这是一个从几何分布 (p=.3) 对大小为 50 的样本进行泊松图的示例:

几何变量的泊松图; 它在开始时向下倾斜,然后有一个扭结并以接近线性的方式向上

如您所见,它显示出明显的“扭结”,表明非线性


泊松图的参考是:

David C. Hoaglin (1980),
“泊松图”,
美国统计学家
卷。34,第 3 期(8 月),第 146-149 页

Hoaglin, D. 和 J. Tukey (1985),
“9. 检查离散分布的形状”,
探索数据表、趋势和形状
(Hoaglin、Moseller 和 Tukey 编辑)
John Wiley & Sons

第二个参考包含对小计数图的调整;你可能想要合并它(但我没有提到手)。


进行卡方拟合优度检验的示例:

除了执行卡方拟合优度之外,通常期望在很多课程中完成它的方式(尽管不是我这样做的方式):

1:从你的数据开始,(我将把它作为我在上面'y'中随机生成的数据,生成计数表:

(x=table(y))
y
 0  1  2  3  4  5  6  7  8  9 10 
 1  2  7 15 19 25 14  7  5  1  4 

2:计算每个单元格中的期望值,假设 ML 拟合泊松:

 (expec=dpois(0:10,lambda=mean(y))*length(y))
 [1]  0.7907054  3.8270142  9.2613743 14.9416838 18.0794374 17.5008954 14.1173890  9.7611661
 [9]  5.9055055  3.1758496  1.5371112

3:注意端类小;这使得卡方分布不太适合作为测试统计量分布的近似值(一个常见规则是您希望期望值至少为 5,尽管许多论文表明该规则具有不必要的限制性;我会接受接近,但一般方法可以适应更严格的规则)。折叠相邻类别,以便最小预期值至少不会低于 5 (在 10 多个类别中,一个预期倒数接近 1 个的类别并不算太糟糕,两个非常接近)。另请注意,我们尚未考虑超出“10”的概率,因此我们还需要将其纳入:

expec[1]=sum(expec[1:2])
expec[2:8]=expec[3:9]
expec[9]=length(y)-sum(expec[1:8])
expec=expec[1:9]
expec
sum(expec) # now adds to n

4:同样,在观察到的情况下折叠类别:

(obs=table(y))
obs[1]=sum(obs[1:2])
obs[2:8]=obs[3:9]
obs[9]=sum(obs[10:11])
obs=obs[1:9]

5:将(可选)连同对卡方的贡献和 Pearson 残差(贡献的有符号平方根)一起放入表中,这些在尝试查看哪里时很有用它不太适合:(OiEi)2/Ei

print(cbind(obs,expec,PearsonRes=(obs-expec)/sqrt(expec),ContribToChisq=(obs-expec)^2/expec),d=4)
  obs  expec PearsonRes ContribToChisq
0   3  4.618   -0.75282      0.5667335
1   7  9.261   -0.74308      0.5521657
2  15 14.942    0.01509      0.0002276
3  19 18.079    0.21650      0.0468729
4  25 17.501    1.79258      3.2133538
5  14 14.117   -0.03124      0.0009761
6   7  9.761   -0.88377      0.7810581
7   5  5.906   -0.37262      0.1388434
8   5  5.815   -0.33791      0.1141816

6:计算,预期总数与观察到的总数匹配损失 1df,估计参数损失 1:X2=i(EiOi)2/Ei

(chisq = sum((obs-expec)^2/expec))
[1] 5.414413
(df = length(obs)-1-1) # lose an additional df for parameter estimate
[1] 7
(pvalue=pchisq(chisq,df))
[1] 0.3904736

诊断和 p 值都显示不缺乏拟合……这是我们所期望的,因为我们生成的数据实际上是 Poisson。


编辑:这里是 Rick Wicklin 博客的链接,其中讨论了泊松图,并讨论了 SAS 和 Matlab 中的实现

http://blogs.sas.com/content/iml/2012/04/12/the-poissonness-plot-a-goodness-of-fit-diagnostic/


Edit2:如果我做对了,1985 年参考文献中修改后的泊松图将是*:

y=rpois(100,5)
n=length(y)
(x=table(y))
k=as.numeric(names(x))
x=as.vector(x)
x1 = ifelse(x==0,NA,ifelse(x>1,x-.8*x/n-.67,exp(-1)))
plot(k,log(x1)+lfactorial(k))

* 他们实际上也调整了截距,但我在这里没有这样做;它不会影响绘图的外观,但是如果您执行参考中的任何其他内容(例如置信区间),如果您执行的操作与他们的方法完全不同,则必须小心。

(对于上面的例子,外观与第一个泊松图几乎没有变化。)

执行卡方拟合优度检验。如果是计数数据,我们可以使用goodfit()包含在 vcd 包中。注意,如果p值大于0.05,我们不能拒绝h0:这个过程是泊松过程。否则,它不是泊松过程。

# load the vcd package
library(vcd) ## loading vcd package

# generate two processes for test
set.seed(2014);y=rpois(200,5)
set.seed(2014);y=rnorm(100, 5, 0.3) # goodfit asks for non-negative values
# output the results
gf = goodfit(y,type= "poisson",method= "ML")
plot(gf,main="Count data vs Poisson distribution")
summary(gf)

# to automatically get the pvalue
gf.summary = capture.output(summary(gf))[[5]]
pvalue = unlist(strsplit(gf.summary, split = " "))
pvalue = as.numeric(pvalue[length(pvalue)]); pvalue

# to mannualy compute the pvalue
chisq = sum(  (gf$observed-gf$fitted)^2/gf$fitted )

df = length(gf$observed)-1-1
pvalue = pchisq(chisq,df)
pvalue