Python模拟二项分布

Python的numpy库中有二项分布采样的函数：

numpy.random.binomial(n,p,size=None)

3个参数：n表示伯努利试验次数，p表示伯努利试验得到正例的概率，size表示采样次数；返回结果为出现正例的次数k。

可以使用该函数来模拟二项分布，借用Inside_Zhang[1]的例子：

野外正在进行9（n=9）口石油勘探井的发掘工作，每一口井能够开发出油的概率是0.1（p=0.1）。请问，最终所有的勘探井都勘探失败的概率？

Python代码：

import numpy as np
from matplotlib import pyplot as plt

# 统计一个列表中每个元素出现的次数
# 参考https://blog.csdn.net/qq_42467563/article/details/86182266.[2]
def count(lis):
    lis=np.array(lis)
    key=np.unique(lis)
    x = []
    y = []
    for k in key:
        mask =(lis == k)
        list_new=lis[mask]
        v=list_new.size
        x.append(k)
        y.append(v)
    return x,y

n,p=9,0.1
size=10000
result=np.random.binomial(n,p,size) #得到采样结果

a,b=count(result)
print(a)
print(b)

# 绘制结果分布图
# 参考https://www.jianshu.com/p/8c0fe1240e78.[3]
plt.figure(1)
plt.subplot(2,1,1)
plt.bar(a,height=b,width = 0.5, color = 'lightskyblue',edgecolor = 'white')
for x,y in zip(a,b):
    plt.text(x, y+0.1, '%d' % y, ha='center', va= 'bottom')
plt.ylim(top=size)
plt.xlim(right=n)
x_ticks = np.arange(0, n+1, 1)
plt.xticks(x_ticks)
plt.ylabel('counter')

# 绘制pmf柱状图
# 参考https://blog.csdn.net/howhigh/article/details/78007317.[4]
plt.subplot(2,1,2)
plt.hist(result, bins=x_ticks, align='left', density=True, rwidth=0.5)  # 绘制直方图
plt.xticks(x_ticks)
plt.xlabel('k')
plt.ylabel('probability')

plt.show()

结果：

[0, 1, 2, 3, 4, 5, 6]
[3885, 3894, 1692, 457, 66, 5, 1]

最终所有的勘探井都勘探失败的概率可以通过"sum(result==0)/size"来近似，结果为：

0.3885

 

参考资料：

[1] Inside_Zhang, “二项分布的实现（np.random.binomial）”, https://blog.csdn.net/lanchunhui/article/details/50172659.

[2] 三尺秋水一点飞鸿, “Python统计一个列表中每个元素出现的次数。四种方法，总有一款适合你”, https://blog.csdn.net/qq_42467563/article/details/86182266.

[3] Kedi, “python_使用matplotlib画柱状图（bar），散点图（scatter）”, https://www.jianshu.com/p/8c0fe1240e78.

[4] howhigh, “9.7 python模拟常用分布”, https://blog.csdn.net/howhigh/article/details/78007317.