上帝的骰子-泊松分布

飞机失事，被雷电击中，被马踢死，枪击案的发生，细胞老化时DNA突变的数目，等交通灯时在你前面的车辆数目，在急诊室中排在你前面的病人数目，某镇上患有白血病的病人数目，某一年中的出生、死亡、结婚、离婚、自杀、他杀的数目，狗身上跳蚤的数目，等等这些小概率事件的发生有什么规律吗? 是上帝掷骰子吗?

爱因斯坦说，我确信上帝是不玩骰子的。所以看似随机的事件也必然有其发生规律。

法国数学家Siméon Denis Poisson于1837年在一本著作中提出了泊松分布的概念，虽然他也不是第一发现人，但后世依然沿用他的名字来命名。泊松分布时一种离散型的分布形态，假设事件是独立发生的且单位时间/空间内事件发生的概率不变，即事件独立且随机发生，则随机变量X服从泊松分布。

60年后，普法战争爆发，一个1898年俄罗斯出生的波兰经济学家Ladislaus Bortkiewicz统计分析了20年间普鲁士士兵被马踢死的数据，试图从中找到规律，发现这是泊松分布的一个完美应用场景, 他在自己的著作"小数定律"中引入了这个经典案例(好像泊松也写过一本"大数定律")

下表记录了从1875年到1894年这20年间，14个军团每年被马踢死的士兵人数 20年14个军团共计280个单位，平均每年每个军团被马踢死0.7人，用python跑一下泊松分布的概率质量函数PMF计算出踢死x个人的概率: 然后用这个概率乘以280得到预期的踢死x个人的次数，再和实际的统计值对比一下，相当吻合:

我想再挑一个数据集验证下，上Kaggle随机找到了土耳其107年的地震记录

数据量比较大，每年大大小小的地震都记录在案了，小地震相对比较频繁，于是选择大于等于里氏6级以上的地震次数进行统计 与前面方式类似，土耳其107年间地震地震97次，平均0.906542次，均值代入PMF计算出占比，然后对比期望值和实际值
基本趋势符合，发生0次地震的次数稍微少了点，发生1次地震的次数稍微多了些，观察原始数据发现1985年到2005年这段数据相当反常，居然连续20年1次地震都没有，不知道是数据集有缺失还是什么其它原因，否则预测趋势将会更贴近实际值。

最后我还抽取了一款耳机在单个门店的销售数据来验证，大多数日子中销量都是0，最多的一天卖了4个出去，日均销量为0.6个。在针对节假日和促销对销量的影响清洗完数据后统计出来的结果也是符合泊松分布的。

下面解释一下前面提到的Poisson概率质量函数和计算过程，PMF公式如下:

e是自然对数，是泊松分布的均值，同样也是泊松分布的方差，因为泊松分布是某种特殊的二项分布，所以这个公式是可以从二项分布公式推导出来的(一个被广泛接受的经验法则是当 n≥20且p≤0.05，可以用泊松分布来估计二项分布值。)。

二项分布公式
当n趋向无穷大而p趋向0的时候，就会变成泊松分布，此时且保持恒定

再结合一个例子说明公式的用法，用于制造核武器和核反应堆的钚239是钚的同位素，一纳克钚239每秒钟发生2.3次放射性衰变，衰变次数符合泊松分布，那么2秒内有3次放射性衰变的概率有多大呢?

原子数量极其庞大，而每个原子在两秒内发生放射性衰变的概率都微乎其微，本质上就是一个n极大而p极小的二项分布。

因为单位时间是两秒，
代入泊松分布的概率密度函数

所以16.3%就是两秒钟内正好发生三次放射性衰变的概率。

如果画图出来可以发现分布呈现右偏，泊松分布的右偏程度取决于的取值，当分布将接近对称，近似正态分布，而接近0的时候，右偏度会非常强。

如果我们要计算每两秒不多于三次放射性衰变的概率，只需要分别计算0,1,2,3次的发生概率再求和就可以了。

所以除了事件要独立且随机发生外，还要发生概率比较低，也就是稀有事件独立且随机发生服从泊松分布，结合下面几个例子来巩固一下:

Case 1: 如果在雨天拿一张很大的纸让雨点落在纸上，会产生下图的独立且随机分布散点，如果随机置放一个正方形，其所覆盖的雨点数是典型的泊松分布。 Case 2: 群居的野生动物分布在广袤的大地上，如果以一百平方米面积的正方形去框定，大概率事件是覆盖0个动物，但是如果一旦框中了则可能数量很集中，在单位空间内的发生概率是不一样的，所以下图不是泊松分布。

Case 3: 每年飞机坠毁次数和空难人数，两者的曲线看起来很相似，但前者符合泊松分布而后者不是，如果你正在一架飞机上，预知到前面有一架飞机会坠毁，但你不会担心自己所在的这架飞机会发生意外，因为这种小概率事件是独立发生的。但是如果你预知到自己隔壁的两个乘客会因为空难而丧生，那对于你来说就是世界末日了。

空难人数要么在某些年份为0，要么就是好几百人，这也会导致值过大，当值大于20的时候泊松分布已经开始向正态分布转化了，当值大于50的时候已经近似正态分布, 下图是=20时一万个样本的分布。

n=10000, lambda=20

from scipy.stats import poisson
import seaborn as sb

data_binom = poisson.rvs(mu=20, size=10000)
ax = sb.distplot(data_binom,
                  kde=True,
                  color='red',
                  hist_kws={"linewidth": 25,'alpha':1})
ax.set(xlabel='Poisson', ylabel='Frequency')

以后有空再写泊松分布与报童模型的结合在单周期库存需求预测中的应用。

References:
Introduction to the Poisson Distribution
What does randomness look like?