1 甜在心馒头店
公司楼下有家馒头店:
每天早上六点到十点营业,生意挺好,就是发愁一个事情,应该准备多少个馒头才能既不浪费又能充分供应?
老板统计了一周每日卖出的馒头(为了方便计算和讲解,缩小了数据):
均值为:
按道理讲均值是不错的选择(参见“如何理解最小二乘法?”),但是如果每天准备5个馒头的话,从统计表来看,至少有两天不够卖, 的时间不够卖:
你“甜在心馒头店”又不是小米,搞什么饥饿营销啊?老板当然也知道这一点,就拿起纸笔来开始思考。
2 老板的思考
老板尝试把营业时间抽象为一根线段,把这段时间用 来表示:
然后把 的三个馒头(“甜在心馒头”,有褶子的馒头)按照销售时间放在线段上:
把 均分为四个时间段:
此时,在每一个时间段上,要不卖出了(一个)馒头,要不没有卖出:
在每个时间段,就有点像抛硬币,要不是正面(卖出),要不是反面(没有卖出):
内卖出3个馒头的概率,就和抛了4次硬币(4个时间段),其中3次正面(卖出3个)的概率一样了。
这样的概率通过二项分布来计算就是:
但是,如果把 的七个馒头放在线段上,分成四段就不够了:
从图中看,每个时间段,有卖出3个的,有卖出2个的,有卖出1个的,就不再是单纯的“卖出、没卖出”了。不能套用二项分布了。
解决这个问题也很简单,把 分为20个时间段,那么每个时间段就又变为了抛硬币:
这样, 内卖出7个馒头的概率就是(相当于抛了20次硬币,出现7次正面):
为了保证在一个时间段内只会发生“卖出、没卖出”,干脆把时间切成 份:
越细越好,用极限来表示:
更抽象一点, 时刻内卖出 个馒头的概率为:
3 的计算
“那么”,老板用笔敲了敲桌子,“只剩下一个问题,概率 怎么求?”
在上面的假设下,问题已经被转为了二项分布。二项分布的期望为:
那么:
4 泊松分布
有了 了之后,就有:
我们来算一下这个极限:
其中:
所以:
上面就是泊松分布的概率密度函数,也就是说,在 时间内卖出 个馒头的概率为:
一般来说,我们会换一个符号,让 ,所以:
这就是教科书中的泊松分布的概率密度函数。
5 馒头店的问题的解决
老板依然蹙眉,不知道 啊?
没关系,刚才不是计算了样本均值:
可以用它来近似:
于是:
画出概率密度函数的曲线就是:
可以看到,如果每天准备8个馒头的话,那么足够卖的概率就是把前8个的概率加起来:
这样 的情况够用,偶尔卖缺货也有助于品牌形象。
老板算出一脑门的汗,“那就这么定了!”
6 二项分布与泊松分布
鉴于二项分布与泊松分布的关系,可以很自然的得到一个推论,当二项分布的 很小的时候,两者比较接近:
7 总结
这个故事告诉我们,要努力学习啊,要不以后馒头都没得卖。
生活中还有很多泊松分布。比如物理中的半衰期,我们只知道物质衰变一半的时间期望是多少,但是因为不确定性原理,我们没有办法知道具体哪个原子会在什么时候衰变?所以可以用泊松分布来计算。
还有比如交通规划等等问题。
这篇文章可以继续扩充:如何理解指数分布?
文章最新版本在(有可能会有后续更新):如何理解泊松分布?
这篇回答节选自我在专栏《机器学习中的数学:概率统计》中的一篇文章,我们来谈一下泊松分布的问题。
欢迎关注我的知乎账号 @石溪 ,将持续发布机器学习数学基础及算法应用等方面的精彩内容。
谈到泊松分布必须要结合着先说说二项分布
1.1.分布列及PMF
我们还是举抛硬币的例子:将一个硬币抛掷 次,每次抛掷出现正面的概率为 ,每次抛掷彼此之间都是相互独立的,随机变量 对应 次抛掷得到的是正面的次数。
这里,随机变量 服从二项分布,二项分布中的核心参数就是上面提到的 和 ,随机变量的分布列可以通过下面这个熟悉的公式计算得到:
下面我们通过依次指定不同的 参数: ,来绘制 图,来观察一下二项随机变量的分布情况:
代码片段:
from scipy.stats import binom import matplotlib.pyplot as plt import seaborn seaborn.set() fig, ax = plt.subplots(3, 1) params = [(10, 0.25), (10, 0.5), (10, 0.8)] x = range(0, 11) for i in range(len(params)): binom_rv = binom(n=params[i][0], p=params[i][1]) ax[i].set_title('n={},p={}'.format(params[i][0], params[i][1])) ax[i].plot(x, binom_rv.pmf(x), 'bo', ms=8) ax[i].vlines(x, 0, binom_rv.pmf(x), colors='b', lw=3) ax[i].set_xlim(0, 10) ax[i].set_ylim(0, 0.35) ax[i].set_xticks(x) ax[i].set_yticks([0, 0.1, 0.2, 0.3]) plt.show()
运行结果:
挺好看的一张图,我们来简要解释一下代码:
第11行:生成服从指定参数 , 的二项分布随机变量。
第12行~第18行: 分别对其进行PMF图绘制,因为是离散型随机变量,因此不建议画成折线图,这种形态更为合适一些。
在这个例子中,我们直接通过scipy中的stats模块得到的二项分布的概率质量函数,也就是反映了不同参数条件下,随机变量 各取值点所对应的取值概率。
1.2.随机变量的采样
我们可以使用 模块中的 方法进行二项随机变量的采样模拟,我们可以指定所要采样的随机变量个数,这里指定重复采样10万次。我们使用三组参数 :分别是 , 和 。
通过上述模拟采样试验可以得到每种实验结果所对应的次数,然后我们通过归一化,可以计算出随机变量每一种取值所对应的频数,并将其作为概率的近似进行绘图观察。
代码片段
from scipy.stats import binom import matplotlib.pyplot as plt import seaborn seaborn.set() fig, ax = plt.subplots(3, 1) params = [(10, 0.25), (10, 0.5), (10, 0.8)] x = range(0, 11) for i in range(len(params)): binom_rv = binom(n=params[i][0], p=params[i][1]) rvs = binom_rv.rvs(size=100000) ax[i].hist(rvs, bins=11, normed=True) ax[i].set_title('n={},p={}'.format(params[i][0], params[i][1])) ax[i].set_xlim(0, 10) ax[i].set_ylim(0, 0.4) ax[i].set_xticks(x) print('rvs{}:{}'.format(i, rvs)) plt.show()
运行结果:
rvs0:[0 4 2 ... 3 2 3] rvs1:[6 6 5 ... 5 7 8] rvs2:[7 8 9 ... 9 7 8]
程序打印的结果是三个数组,这就是我们在不同参数下分别做10万次采样试验的结果数组。
1.3.随机变量的数字特征
服从二项分布的随机变量,他的期望和方差的表示很简单,服从参数为 的二项分布的随机变量 ,他的期望和方差的公式我们直接给出来:
期望:
方差:
我们可以结合上面的试验,用几种方法来验证一下上述结论:
代码片段:
import numpy as np from scipy.stats import binom binom_rv = binom(n=10, p=0.25) mean, var, skew, kurt = binom_rv.stats(moments='mvsk') binom_rvs = binom_rv.rvs(size=100000) E_sim = np.mean(binom_rvs) S_sim = np.std(binom_rvs) V_sim = S_sim * S_sim print('mean={},var={}'.format(mean,var)) print('E_sim={},V_sim={}'.format(E_sim,V_sim)) print('E=np={},V=np(1-p)={}'.format(10 * 0.25,10 * 0.25 * 0.75))
运行结果:
mean=2.5,var=1.875 E_sim=2.50569,V_sim=1.8735076238999997 E=np=2.5,V=np(1-p)=1.875
我们用三种方法计算了服从参数为 的二项分布随机变量的均值和方差,其中:
第04行~第05行: 是用函数包中的方法计算的分布的各个理论统计值;
第07行~第10行:从采样试验中得到的样本数据计算出来的均值和方差;
第14行:通过公式直接计算出来的理论值。
看的出,利用采样样本数据计算出来的值和理论值基本上是相等的。
看完了二项分布,我们再来看看泊松分布:
关注 @石溪 知乎账号,分享更多机器学习数学基础精彩内容。
2.1.泊松分布的应用场景
我们刚刚讲了, 次独立的伯努利试验成功的次数是一个服从二项分布的随机变量,其中参数为 和 ,期望为 。我们这里看一种非常特殊的情况:就是 非常大, 非常小,但是期望 结果适中。
现实生活中有没有这类情况?有,比如我们考虑任何一天内发生飞机事故的总数,记作随机变量 ,总共飞机飞行的次数 非常大,但是单架次飞机出现事故的概率 非常小。或者用随机变量 表示一本书中字印刷错误的次数, 表示一本书中的总字数,非常大,而 表示每个字印刷出错的概率,非常小。
这种情况下, 很大 很小,二项分布的分布列可以简化为我们这里谈到的泊松分布的分布列:
,其中, ,
期望和方差满足:
特别的,当我们的 ,且 时,对应的二项分布列:
就收敛于上面的泊松分布列了。
通俗点说把,就是只要当 ,且 非常大, 非常小,泊松分布就是二项分布的一个非常好的近似。计算简便就是他的一个很大的优势
2.2.泊松分布的PMF图
同样的,我们也用python代码来画一下他的PMF函数图,对应的观察一下指定参数下泊松分布的分布列。
正如我们所说,泊松分布的参数就是一个 ,我们分别绘制一个 和 的泊松分布PMF图,并获取他们的均值和方差。
代码片段:
from scipy.stats import poisson import matplotlib.pyplot as plt import seaborn seaborn.set() fig, ax = plt.subplots(2, 1) x = range(0, 20) params = [10, 2] for i in range(len(params)): poisson_rv = poisson(mu=params[i]) mean, var, skew, kurt = poisson_rv.stats(moments='mvsk') ax[i].plot(x, poisson_rv.pmf(x), 'bo', ms=8) ax[i].vlines(x, 0, poisson_rv.pmf(x), colors='b', lw=5) ax[i].set_title('$\lambda$={}'.format(params[i])) ax[i].set_xticks(x) print('lambda={},E[X]={},V[X]={}'.format(params[i], mean, var)) plt.show()
运行结果:
lambda=10,E[X]=10.0,V[X]=10.0 lambda=2,E[X]=2.0,V[X]=2.0
同样的,我们对 的泊松分布进行采样。
代码片段:
import numpy as np from scipy.stats import poisson import matplotlib.pyplot as plt import seaborn seaborn.set() lambda_ = 2 data = poisson.rvs(mu=lambda_, size=100000) plt.figure() plt.hist(data, normed=True) plt.gca().axes.set_xticks(range(0, 11)) print('mean=', np.mean(data)) print('var=', np.square(np.std(data))) plt.show()
运行结果:
mean= 2.00542 var= 2.0082906236
这也是我们通过10万次采样试验得出的统计结果,我们通过这个结果集计算了均值和方差,和模型的理论推导值是一致的。
当然还有《机器学习中的数学(全集)》系列专栏,欢迎大家阅读,配合食用,效果更佳~
有订阅的问题可咨询微信:zhangyumeng0422
英国:山姆先生,你好你好,有什么事情我们可以帮你的?
美国:我要说的事,你们千万别害怕
英国:我们是日不落的高卢鸡,我们不会怕,你请说。
美国:这次5G,我不能监控了!
英国:5G是哪一位?
美国:不是哪位,是那个我完全不能走后门的5g啊!
法国:(画画中,2g)
美国:没那么老,刚出的
法国:(画画中,3g)
美国:没那么慢,网速几百兆的
法国:(画画中,4g)
美国:延迟呢?
法国:(转过来,4g+)
美国:这……
英国:(添几笔,韩国5g)
美国:(愤怒!!)5g啊,就那个华为出的5g啊,新闻有没有看?最近很火,速度起飞的那个啊,明白吗
英国:明白了,你请说
美国:他们说我监控别人,搞五眼联盟,试问谁不知道,他们就独立搞研发去了,就在中国深圳,一下子就搞出来了,一个基站就茶杯那么大,茶杯啊。我刚要制裁他,他就退出美国,还起诉我,对我的制裁理都不理,我就像个傻…
法国:(憋不住笑了)
美国:你在笑什么?
法国:我想起高兴的事情
美国:什么高兴的事情?
法国:我刚卖了300架飞机
英国:(憋不住了)
美国:你又笑什么?
英国:我也卖了300架飞机
美国:你们卖的是同样的飞机?
法国:对对对
英国:啊,不是,刚才直播了一下,把一个月流量烧光了。
美国:我再重申一遍,我没有开玩笑!
法国:对对
美国:喂~~
英国:哎,我们言归正传啊,那个,这个华为5g,他厉害吗?
美国:他不是厉不厉害的问题,他真的是那种,那种很少见的那种,敢跟我对着搞,速度快,延迟低,而且我还走不了后门,遗憾的是不是我们搞出来的,我还没来得及体验……
法国:(偷笑)
美国:你个高卢鸡还在笑,我忍你很久了
法国:我卖了飞机
美国:你明明在笑都没停过!
法国:山姆先生,我们受过严格的训练,无论多好笑我们都不会笑,除非忍不住。
英国:不如这样,你先回去用4g发推,等我直播爽完了就来帮你谴责他们。
美国:行,你们赶紧直播,完了多带点水军,很危险的。
(转b站评论