1 2 3

泊松分布的仿真及可视化

本文介绍泊松分布和泊松函数的定义,并通过 Python random 库对泊松分布进行仿真,带你触摸复杂表象下的简单本质。

GitHub项目地址:python-tips/poisson

1. 从泊松函数讲起

泊松分布 表示在给定时间段内发生给定数量的事件的概率。这个定义比较抽象。举个具体的例子,假设你每小时接到电话的概率是固定的,比如每小时 0.05 个,那么你在接下来 1 小时内接到电话个数的概率,就服从泊松分布:

1 小时内接到 0 个电话,对应一个概率值 $P_0$

1 小时内接到 1 个电话,对应一个概率值 $P_1$

… …

1 小时内接到 n 个电话,也对应一个概率值$P_n$

这些概率值组成一个概率分布列,它们的值 $(n, P_n)$ 在二维坐标下连成一条曲线。这条曲线所在的函数就是泊松分布的概率密度函数。其公式及图像如下:

$$\boxed{P(k | t, \lambda)=\frac{(\lambda t)^{k}}{k !} \exp (-\lambda t)}$$

从公式中,我们可以看出:只要确定了 $\lambda$$t$,该式就退化成了概率 $P$ 关于事件发生次数 $k$ 的函数。 类似地,如果我们确定了 $\lambda$$k$,则该式退化成概率$P$ 关于时间范围 $t$ 的函数。

“确定哪些参数,让函数最终退化成哪些参数的函数”,这个选择和我们的研究目的有关。如果你对不同 $k$ 如何影响 $P$ 值感兴趣,那么就应该确定参数 $\lambda$$t$。如果对 $t$$P$ 之间的关系感兴趣,那么就应该确定参数 $\lambda$$k$

$\lambda$, $k$, $t$ 的定义:

  • $\lambda$: 单位时间内,事件发生的频率。
  • $k$: 事件发生次数。
  • $t$: 观测事件发生次数的时间范围。

$\lambda$, $k$, $t$ 三者的关系:

泊松分布衡量的是多长时间内,某事件发生多少次的概率。这里 $t$ 指代的是多长时间$k$ 指代的是某事件发生多少次$\lambda$ 则类似该事件的一个固有属性$\lambda$ 越大,可以简单理解为该事件在一段时间内发生的概率越大。

Note: $t$ 的量纲必须和 $\lambda$的量纲相对应。如果 $\lambda$ 的单位是 次/小时,则 $t$ 应该取 小时 为单位;如果 $\lambda$ 的单位是 次/分钟,则 $t$ 也应该取 分钟 为单位。

2. 引入一个实例

接下来我们用一个直观的实例,解释泊松分布在实际问题中是如何运用的。

Q: 一个淘宝客服,平均每 10 分钟接听 1 个电话,那么他连续 1 小时没接到任何电话的概率是多少?

A:

$P\left(k=0 | t=1, \lambda=6\right)=\exp \left(-6\right)=0.0024$

其中,$t = 1 小时,\lambda = 6 次/小时$

我们发现,这个概率是极低的。这说明该客服也许正在上夜班,或者电话断线了,总之他很可能处在一种非正常的状况下。否则在正常情况下,连续 1 小时没接到任何电话的概率仅有千分之二,是极其罕见的。

3. 泊松分布的本质

尽管泊松分布的函数形式看起来很复杂,但它本质上其实很简单。泊松函数的本质,也就是它的基本假设,可以追溯到一个简单的公式:

$$\boxed{P = \lambda \Delta t}$$

这个公式看起来太过简单,以至于你可能不相信它能推导出上文那个复杂的泊松分布函数。如果你想了解推导过程,可以看我之前写的博客 排队论在网络性能分析中的应用,里面有详细证明。本文的主题不在与于,就不展开讲了。

4. 建模仿真

这个简单公式,就是泊松分布的核心假设。也就是说,我们只需要用这个公式,就能对泊松分布做仿真了。

这里还有个技巧。因为在假设中,$\Delta t$ 代表无限小的时间间隔,但在编程中,我们没法写出无限小这种东西(其实可以,求别杠)。因此只能量力而行。比如,如果 $t$ 代表一个小时,而 $\lambda$ 表示 20 次每小时,则我们的时间片只要切到比 20 这个数字大一到两个数量级即可。20 是的数量级是十,大两个数量级是千。这里我们把一小时时间切割为一千个时间片,并将每个时间片内事件发生的概率设为 201000 ,就保证了一小时内,事件发生次数的期望是 20 了。也就在千分之一小时的精度内,满足了 $P = \lambda \Delta t$ 这个公式。

OK,闲言少叙,我们这就开始建模仿真。

# -*- coding: utf-8 -*-

"""泊松仿真
   author: [email protected]
   usage: python poisson.py [RATE] {TIME}
"""

import sys
import random
import collections
import matplotlib.pyplot as plt


class Poisson:
    def __init__(self, rate=sys.argv[1]):
        self.rate = int(rate)

        self.time = 1  # 单位时间
        if len(sys.argv) > 2 and sys.argv[2] != "":
            self.time = int(sys.argv[2])  # 手动指定时间范围

        self.EXP_NUM = 100000  # 实验次数
        self.NUM_LEVEL = 2  # 数量级

    def generator(self, prob):
        """仿真结果生成器"""
        while True:
            if random.random() < prob:
                yield 1
            else:
                yield 0

    def perform_exp(self, rate, time):
        """进行一次实验
        每次实验中,时间分片的数量比rate高两个数量级
        """
        level = len(str(rate))
        shard_num = 10 ** (level + self.NUM_LEVEL)  # 计算时间分片的数量

        gen = self.generator(rate / shard_num)

        cnt = 0
        for _ in range(time * shard_num):
            cnt += next(gen)
        
        return cnt

    def perform_exps(self, exp_num, rate, time):
        """多次实验,得到分布"""
        lst = []
        for _ in range(exp_num):
            lst.append(self.perform_exp(rate, time))

        return sorted(collections.Counter(lst).items(), key=lambda e: e[0])

    def draw(self, sorted_list):
        """画图"""
        s = sum([e[1] for e in sorted_list])
        x = [e[0] for e in sorted_list]
        y = [e[1] / s for e in sorted_list]

        plt.plot(x, y)
        
        plt.xlabel("k")
        plt.ylabel("P(k)")
        plt.show()

    def main(self):
        sorted_list = self.perform_exps(self.EXP_NUM, self.rate, self.time)
        self.draw(sorted_list)

    @staticmethod
    def calculator(rate, t, k):
        """用于计算泊松函数的概率 P(k|t,lambda)
           rate: lambda
           t: t
           k: k
        """
        import math
        return (rate * t) ** k / math.factorial(k) * math.exp(-rate * t)


if __name__ == '__main__':
    p = Poisson()
    p.main()

看完代码,聪明的你一定已经知道这段代码在干什么了

那我们直接来说如何使用代码吧!

如果你想知道 $\lambda = 3, t = 2$ 情况下泊松函数的图像,在命令行执行:

python poisson.py 3 2

如果不设置 $t$ 的值,则 $t$ 的值默认为 1。

下面放一些仿真结果:

  1. $\lambda = 2, t = 1$
  2. $\lambda = 5, t = 1$
  3. $\lambda = 10, t = 1$
  4. $\lambda = 80, t = 1$

5. 三维可视化

泊松分布有三个自变量,一个因变量,需要四维空间才能把它完整地绘制出来。但幸好,通过观察我们发现 $\lambda$$t$ 的关系很近,可以归约成一个变量 $\mu$。我们设 $\mu = \lambda t$,则绘制 $P = f(\mu, k)$ 的图像,就足以表达泊松函数四个变量间的函数关系了。

Note: 绘图代码见 3dplot.py

6. 泊松的应用

泊松分布在生活中也有很多应用,下面我举一个比较贴近生活的例子。

停车问题:公司有个自行车停车场,临近上班就停满了车。停车区域是一个长条形,从这头走到那头的过程中,我们需要找到合适的车位将车停入。但是,车位有宽有窄,我们倾向于选择比较宽的车位。我们需要一种判断方式,来决定何时将车停入是合适的。

为了简化问题,我们做出如下假设:

  • 假设1:车辆在停车区域内均匀分布
  • 假设2:车辆到达的速率是一个常数 $\lambda$
  • 假设3:假设在某个时间点 T,停车场被停满。此刻,我们让时间倒流,由于车辆停放的逆过程是车辆离开。因此对于时间 T - x,车辆离开的过程可视作时间为 x 离开速率为 $\lambda$ 的泊松过程。

现在我们提出某种策略:首先将长条停车场均分为前后两段。前段用于观察车辆的分布,后段用于实际停放车辆。假设我们在前段发现有 5 个勉强能将车停入的车位,则大胆估计前段停车场的车位数量服从 $\lambda t = 5$ 的泊松分布。由于存在假设1,因此我们可以近似认为后段停车场的车位数量也服从 $\lambda t = 5$ 泊松分布。

$\lambda t = 5$ 泊松分布在前文中已经出现过了,它的图像如下:

由上图可知,在后段停车场中,遇不见类似之前 5 个车位那样,能将车勉强停入的车位的概率仅为 0.7% 左右,这是真的是非常低的概率。尽管这是一个近似的结果,但也能帮助我们判断,目前停车场的车位状况是比较宽松的。这时我们就可以有信心找个好车位停车,而不是找个狭窄的位置硬挤。