0%

泊松分布及其应用

说明
1. 译自Aatish Bhatia的博文What does randomness look like?
2. 译自Aatish Bhatia的博文Are mass shootings really random events? A look at the US numbers.
3. 部分借鉴自马同学高等数学中的文章“如何理解泊松分布?”

一、泊松分布的推导

(一)将问题转变为二项分布

楼下有家馒头店,每天早上六点到十点营业,生意挺好,就是发愁一个事情,应该准备多少个馒头才能既不浪费又能充分供应?老板统计了一周每日卖出的馒头(为了方便计算和讲解,缩小了数据):

均值为:\[\overline {X}=\cfrac {3+7+4+6+5}{5}=5\]

按道理讲均值是不错的选择,但是如果每天准备5个馒头的话,从统计表来看,至少有两天不够卖,40%的时间不够卖:

老板尝试把营业时间抽象为一根线段,把这段时间用T来表示:

然后把周一的三个馒头按照销售时间放在线段上:

把T均分为四个时间段:

此时,在每一个时间段上,要不卖出了(一个)馒头,要不没有卖出:

在每个时间段,就有点像抛硬币,要不是正面(卖出),要不是反面(没有卖出):T内卖出3个馒头的概率,就和抛了4次硬币(4个时间段),其中3次正面(卖出3个)的概率一样了。

这样的概率通过二项分布来计算就是:

\[\begin{pmatrix} 4 \\ 3 \end{pmatrix}p^{3}\left( 1-p\right) ^{1}\]

但是,如果把周二的七个馒头放在线段上,分成四段就不够了:

从图中看,每个时间段,有卖出3个的,有卖出2个的,有卖出1个的,就不再是单纯的“卖出、没卖出”了。不能套用二项分布了。

解决这个问题也很简单,把T分为20个时间段,那么每个时间段就又变为了抛硬币:

这样,T 内卖出7个馒头的概率就是(相当于抛了20次硬币,出现7次正面): \[\begin{pmatrix} 20 \\ 7 \end{pmatrix}p^{7}\left( 1-p\right) ^{13}\]

为了保证在一个时间段内只会发生“卖出、没卖出”,干脆把时间切成n份: \[\begin{pmatrix} n \\ 7 \end{pmatrix}p^{7}\left( 1-p\right) ^{n-7}\]

越细越好,用极限来表示: \[\lim _{n\rightarrow \infty }\begin{pmatrix} n \\ 7 \end{pmatrix}p^{7}\left( 1-p\right) ^{n-7}\]

更抽象一点,T时刻内卖出k个馒头的概率为: \[\lim _{n\rightarrow \infty }\begin{pmatrix} n \\ k \end{pmatrix}p^{k}\left( 1-p\right) ^{n-k}\]

(二)计算概率\(p\)

二项分布的期望为: \[E(X)=np=\mu\] 那么:\[p=\cfrac{\mu}{n}\]

(三)泊松分布的推导

有了\(p=\cfrac{\mu}{n}\)了之后,就有:\[\lim _{n\rightarrow \infty }\begin{pmatrix} n \\ k \end{pmatrix}p^{k}\left( 1-p\right) ^{n-k}=\lim _{n\rightarrow \infty }\begin{pmatrix} n \\ k \end{pmatrix}{(\cfrac{\mu}{n})}^{k}\left( 1-\cfrac{\mu}{n}\right) ^{n-k}\]

我们来算一下这个极限: \[\lim _{n\rightarrow \infty }\begin{pmatrix} n \\ k \end{pmatrix}{(\cfrac{\mu}{n})}^{k}\left( 1-\cfrac{\mu}{n}\right) ^{n-k}=\lim _{n\rightarrow \infty }\cfrac {n\left( n-1\right) \left( n-2\right) \ldots \left( n-k+1\right) }{k!}{(\cfrac{\mu}{n})}^{k}\left( 1-\cfrac{\mu}{n}\right) ^{n-k}\] \[=\lim _{n\rightarrow \infty }\cfrac {\mu ^{k}}{k!}\cfrac {n}{n}\cdot \cfrac {n-1}{n}\ldots\cfrac {n-k+1}{n}\left( 1- \cfrac {\mu }{n}\right) ^{-k}\left( 1-\cfrac {\mu }{n}\right) ^{n}\]

\[\lim _{n\rightarrow \infty }\left( 1-\cfrac {\mu }{n}\right) ^{n}=e^{-\mu}\]

\[\lim _{n\rightarrow \infty }\cfrac {n}{n}\cdot \cfrac {n-1}{n}\ldots\cfrac {n-k+1}{n}\left( 1-\cfrac {\mu }{n}\right) ^{-k}=\lim _{n\rightarrow \infty }1\cdot(1-\cfrac {1 }{n})\cdot(1-\cfrac {2 }{n})\ldots(1-\cfrac {k-1 }{n})=1\]

所以: \[\lim _{n\rightarrow \infty }\begin{pmatrix} n \\ k \end{pmatrix}{(\cfrac{\mu}{n})}^{k}\left( 1-\cfrac{\mu}{n}\right) ^{n-k}=\dfrac {\mu ^{k}}{k!}e^{-\mu }\]

上面就是泊松分布的概率密度函数,也就是说,在T时间内卖出k个馒头的概率为: \[P(X=k)=\dfrac {\mu ^{k}}{k!}e^{-\mu }\]

一般来说,我们会换一个符号,让\(\mu = \lambda\),所以: \[P(X=k)=\dfrac {\lambda ^{k}}{k!}e^{-\lambda }\]

这就是教科书中的泊松分布的概率密度函数。泊松分布是单位时间内独立事件发生次数的概率分布

(四)解决馒头店的问题

老板依然蹙眉,不知道\(\mu\)啊?没关系,刚才不是计算了样本均值:\[\overline X=5 \] 可以用它来近似:\[ \overline X \approx \mu\] 于是: \[P(X=k)=\dfrac {5 ^{k}}{k!}e^{-5 }\]

画出概率密度函数的曲线就是:

可以看到,如果每天准备8个馒头的话,那么足够卖的概率就是把前8个的概率加起来:

(五)二项分布与泊松分布的关系

鉴于二项分布与泊松分布的关系,可以很自然的得到一个推论,当二项分布的\(p\) 很小的时候,两者比较接近。

(六)Poisson Process的定义:

假定一个事件在一段时间内随机发生,且符合以下条件:

  1. 将该时间段无限分隔成若干个小的时间段,在这个接近于零的小时间段里,该事件发生一次的概率与这个极小时间段的长度成正比。
  2. 在每一个极小时间段内,该事件发生两次及以上的概率恒等于零。
  3. 该事件在不同的小时间段里,发生与否相互独立。

则该事件称为poisson process。

二、泊松分布的前世今生:随机性是什么样的?

1944年6月13日,诺曼底(Normandy)登陆后的一个星期,一阵嗡嗡声呼啸着划破饱受战火的伦敦上空,这种响声来自当时德国发明的战争武器:V-1飞行炸弹。作为巡航导弹的前身,V-1是一种自我推进、由陀螺仪导航、通过简单脉冲式喷气发动机以每秒50次的频率吸进空气点燃燃料提供动力的飞行炸弹。由于高频率的喷气使得这种炸弹发出独特的声音,因此它有了一个绰号:“嗡嗡炸弹(buzzbombs)”。

从1944年的6月到10月,德国从法国海岸和荷兰总共发射了9521枚这种嗡嗡炸弹,其中有2419枚击中了伦敦的目标。英国人担忧这种无人驾驶战机的精确性。它们只是随机地飞过城市,还是会击中既定的目标?德国人真的发明了一种能自我导航且命中率高的炸弹吗?

幸运的是,英国人小心谨慎地统计了二战期间落到伦敦的几乎所有这种炸弹的地点和时间。根据这些数据,他们可以统计得出到底这些炸弹是随机落到伦敦还是被瞄准发射的,这是一个事关现实后果的数学问题。

想象一下,此时此刻,你供职于英国情报处,你被要求解决这个问题。某人给了你一张布满密密麻麻的点的纸,而你的任务就是判断这些点是否随机。让我们具体一点。现在有来自于Steven Pinker的书《The Better Angels of our Nature》中的两个图案,一个是随机生成的,另一个是模拟自然生成的,你能分辩出来吗?

以下是Pinker的解释:

“左边的图案,有块状的、条状的、孔状的、丝状的(又或者是动物形状、人物裸体像、甚至是圣玛利亚——全凭你的想象),是随机生成的阵列,如星星;而右边的图案,看起来像是杂乱无章的,是由相互间保持一定距离的粒子位置生成的阵列,如萤火虫。”

没错,就是萤火虫。右边图案的点记录的是新西兰ceiling of the Waitomo cave中萤火虫的位置。这些萤火虫并非随机排列的,他们正在相互推挤着抢食物,防止族群抢夺自己的既得利益。

试着均匀地往地面撒沙子,结果会像右边的图案。你本能地避开了已经撒过沙子的地方。随机撒沙子的过程并不带有主观偏见,沙子只会简单地落在该落地方,成堆成组等等,就像你闭着眼睛撒沙子一样。关键的差别在于,随机性并不像均一性。真正的随机性也会产生群组,正如我们在夜空所看到的星座一样。

还有另一个例子。想象一下,一位教授要求她的学生掷100次硬币。一个学生用心地去做,并把结果记下来。另一个学生有点偷懒,不做实验而是编造假的掷硬币结果。你能分辨出哪个学生偷懒吗?

第一个学生的数据有群组:一连长达8次的硬币反面。这可能看起来出人意料,但这的确是随机掷硬币得到的结果(我知道这是事实——因为我掷了100次得到那样的数据!)。第二个学生的数据缺乏群组性,这有点可疑。实际上,在100次的掷硬币过程中,没有一连4次或以上的正面或反面,这种情况大概只有0.1%的概率会发生,表明第二个学生在捏造数据(这确实是我捏造的)。

如何判断出一组数字模式是否随机可能看起来有点像晦涩难懂的数学游戏,但这却让我们离事实更接近。关于随机波动性的研究起源于19世纪法国犯罪统计学。当时法国正在迅速城市化,城市的人口密度猛增,犯罪与贫穷成为紧迫的社会问题。

在1825年,法国开始收集犯罪审判的统计数据。随之而来的或许是统计分析应用于社会问题研究的第一个例子。比利时数学家Adolphe Quetelet是当时社会科学的早期先驱者。他具有争议性的目标是把天文学中的概率思想应用于研究治理人类的法律。

Michael Maltz 认为:

在寻找犯罪统计与天文学观察中相同的规律性的过程中,他认为,正如星星有一个具体的位置(不考虑定位测量方法的差异性),我们的社会同样存在着一个犯罪率水平。他构造了“平常人”和“道德人”的概念,并断定平常人具有一种统计意义上的“持续的犯罪倾向”。这使得“社会物理学家”能够计算出随着时间推移的轨迹,“能够揭示出简单的运行规律和预测未来”。(Gigerenzer et al, 1989)

Quetelet注意到犯罪率随时间推移而缓慢下降,并推断法国市民的犯罪倾向必定存在下降趋势。他所用的数据存在一些问题,而他的方法中的关键漏洞被才华横溢的法国博学家和科学家西莫恩•德尼•泊松(Siméon-Denis Poisson)发现了。

泊松的想法既巧妙且现代,用今天的语言来说,他认为Quetelet的数据缺乏一个“模型”。Quetelet没有考虑到陪审员是如何得出审判结果的。根据泊松的说法,陪审员是容易犯错的,我们观察到的数据是定罪率,但我们想知道的是被告人有罪的概率,这两个量并不一样,但可以关联起来。结果是,当把以上过程考虑进去,定罪率会表现出固有的某种程度的波动性,这在法国的犯罪数据中可以看到。

1837年,泊松把研究成果发表在“Research on the Probability of Judgments in Criminal and Civil Matters”。在该论文中他提出了我们现在所说的“泊松分布”公式。文中讲述了大量随机事件发生具体次数的概率(如大多数的法国陪审员做出错误判决的概率)。例如,我们假设平均一年会有45人被雷电击中。把这个数据和人口数量运用到泊松公式中,可以得到一年中有10人、50人或者100人被雷电击中的概率。假设条件是雷击是相互独立且罕见的事件,并且会在任意时间等可能发生。换句话说,泊松公式能够告诉你只因偶然性而导致罕见事件发生的概率。

泊松公式的首次应用来自于一个不太相关的地方。往后过60年,穿过普法战争,来到1898年的普鲁士。一位具有波兰血统的普鲁士统计学家Ladislaus Bortkiewicz正尝试弄明白为什么在某些年份普军中会有异常多的士兵被马踢死。在一个军团中,有时一年中会有4个士兵被马踢死。这只是巧合吗?

发生一次被马踢死的事件是罕见的(并且相互独立的,除非马匹被秘密安排了这样做)。Bortkiewicz意识到他可以运用泊松公式计算出死亡事件预计发生的次数。中间列是预测结果,右侧列是实际数据。

Number of Deaths by Horse Kick in a year Predicted Instances (Poisson) Observed Instances
0 108.67 109
1 66.29 65
2 20.22 22
3 4.11 3
4 0.63 1
5 0.08 0
6 0.01 0

可以看出两组数据如此的吻合吧?如果被马踢死这种事件是纯随机过程,那么零散(sporadic)的死亡次数群组性是可以预测得出的。随机性伴随着群组性(Randomness comes with clusters)。

我决定自己亲自试一下。我在寻找因罕见事件导致死亡的公开可用的数据集的过程中,发现了这份把全世界鲨鱼袭击人类事件列成表格的国际鲨鱼袭击文档(International Shark Attack File)。以下是南非鲨鱼袭击数据。

Year Number of Shark Attacks in South Africa
2000 4
2001 3
2002 3
2003 2
2004 5
2005 4
2006 4
2007 2
2008 0
2009 6
2010 7
2011 5

数字都比较小,平均值为3.75。但比较2008年和2009年的鲨鱼袭击数据,一年为0,而下一年为6,到了2010年则为7。你可以想象出新闻头条在大声报道“鲨鱼袭击!”。次数这么多,是因为鲨鱼的反攻呢,还是因为偶然性而导致的呢?为了查找原因,我把实际数据与泊松预测结果进行对比。

蓝色的是每年观察到的发生0、1、2、3……次鲨鱼袭击的数据。例如,蓝色长条代表发生4次鲨鱼袭击的年份(2000,2005和2006)。红色点线是泊松分布,代表的是当鲨鱼袭击是纯随机过程时可以预料得到的发生次数。数据吻合得很好:我没有找到任何数据群是超出泊松分布预测的(p=0.87)。这恐怕排除2010年南非鲨鱼袭击数量激增的可能性。再一次证明随机事件并不是随意发生的。

让我们回到嗡嗡炸弹的例子。以下是掉落到不同区域的炸弹个数的形象化表示,由Charles Franklin根据在British Archives in Kew的原始地图重现。

注:澄清说明。上图显示的是掉落到伦敦的炸弹的分布。我现在问的问题是,如果你把城市受到严重攻击的区域放大来看(尤其是上图中的高峰区),那么炸弹是受到更精确的操纵而击中明确的目标吗?这远远不是均匀分布,但它显示出了精确瞄准的迹象了吗?现在你应该猜出如何回答这个问题了。在一篇题为“An Application of the Poisson Distribution”的报告中,一位名叫R. D. Clarke的英国统计学家写道:

在飞行炸弹袭击伦敦的过程中,人们普遍声称炸弹的袭击地点倾向于聚集一起。于是,一项统计测试被用来发掘是否存在支持这种主张的论据。

Clarke把南伦敦中被严重袭击的12km$$12km的区域划分为网格。他总共划分了576个方格子,每个方格子大约是25条城市街区。接着,他算出被0个、1个、2个……炸弹击中的方格子个数。

总共有537个炸弹落到这576个方格子。平均每个方格子被略少于1个炸弹击中。他把这个参数代入泊松公式中去计算出因为偶然性而导致可以预料的炸弹聚集性大小。以下是取自他的论文中相关的表格:

对比两栏数据,你可以看到预测数据与实际数据是如此惊人地接近。有7个方格子被4个炸弹击中——但这是根据偶然性能够预料得到的。在伦敦的大范围之内,炸弹并没有被瞄准。它们只是大量地随机击下来,犹如在玩一个毁灭性的、城市范围内的俄式轮盘游戏。

泊松分布已经渗透到各种各样的应用场景中,其中一些无关紧要,而另一些则性命攸关。如细胞老化时DNA突变的数目,等交通灯时在你前面的车辆数目,在急诊室中排在你前面的病人数目,我博客中错别字的数目,某镇上患有白血病的病人数目,某一年中的出生、死亡、结婚、离婚、自杀、他杀的数目,你的狗身上跳蚤的数目,等等。

从普通平凡的小事到关乎生死的大事,这些维多利亚时代的科学家们教会了我们,随机性发挥着比我们所承认的更大的作用。可悲的是,这个事实并不能在我们生活不如意时提供多少的慰藉。

“So much of life, it seems to me, is determined by pure randomness.(以我所看,生活的很大部分是由随机性主导的)” – Sidney Poitier

三、美国枪击案

资料显示,1982年至2012年,美国共发生62起(大规模)枪击案。其中,2012年发生了7起,是次数最多的一年。

2012年有这么多枪击案,这是巧合,还是表明美国治安恶化了?假定美国枪击案满足"泊松分布"的三个条件:

  • 枪击案是小概率事件。
  • 枪击案是独立的,不会互相影响。
  • 枪击案的发生概率是稳定的。

显然,第三个条件是关键。如果成立,就说明美国的治安没有恶化;如果不成立,就说明枪击案的发生概率不稳定,正在提高,美国治安恶化。

根据资料,1982-2012年枪击案的分布情况如下:

计算得到,平均每年发生2起枪击案,所以λ = 2 。

上图中,蓝色的条形柱是实际的观察值,红色的虚线是理论的预期值。可以看到,观察值与期望值还是相当接近的。

我们用"卡方检验"(Chi-Squared Test),检验观察值与期望值之间是否存在显著差异,得到p-value为0.09,这是什么意思?这意味着,美国的治安没有恶化,枪击案的数量并没有超过随机过程的预期。换句话说,1982-2012年发生的大规模枪击事件与枪击事件是独立事件的假设并不矛盾,大规模枪击事件平均每年发生2次。

然而,p值为0.09并不是特别高,如果随后继续出现与2012年一样的极端年份,很可能会排除大规模枪击事件是随机事件的假设。

\[卡方统计量=\sum \cfrac {\left(A-E\right)^2}{E}\]

E为观察值
A为观察值