Numpy计算大数的二项式概率分布

在本文中，我们将介绍如何使用NumPy计算大数的二项式概率分布。二项式分布是用来计算二元事件发生次数的概率分布。在实际问题中，我们经常需要计算二元事件发生的概率。如服从二项式分布的事件可以是扔硬币正面朝上的次数，或者在一次抛掷两个骰子时，其中一个骰子的点数为6的次数等等。下面，我们将来看看二项式分布的计算方法。

阅读更多：Numpy 教程

二项式概率分布的计算方法

二项式分布的定义式为：
$P(X=k)=C_n^kp^k(1-p)^{n-k}$

其中， $P(X=k)$ 是二项式分布的概率。 $C_n^k$ 表示从 $n$ 个元素中选择 $k$ 个元素的组合数，即：
$C_n^k=\frac{n!}{k!(n-k)!}$

$p$ 是事件发生的概率， $1-p$ 则是事件不发生的概率。 $n$ 是试验次数， $k$ 是事件发生的次数。

我们来看一个具体的例子，假设我们翻20次硬币，求硬币正面朝上恰好13次的概率。那么我们可以使用下面的代码：

import numpy as np

n = 20
p = 0.5
k = 13

P = np.math.comb(n, k) * (p ** k) * ((1 - p) ** (n - k))
print(P)

运行结果为：

0.09802913665771484

即硬币正面朝上恰好13次的概率为0.098。

不过，当试验次数很大，或者事件发生的次数很多时，直接计算组合数和幂次计算会非常困难甚至不可能完成。这时，我们可以使用NumPy中提供的二项式分布函数来计算。

NumPy二项式分布的使用方法

NumPy中的二项式分布函数是numpy.random.binomial(n, p, size=None)，它可以生成一个满足二项式分布的随机样本。

使用方法很简单，我们来看一个例子：我们希望模拟10万次抛掷硬币，每次抛掷的概率为0.5，求硬币正面朝上恰好50000次的概率。

import numpy as np

n = 100000
p = 0.5
k = 50000

x = np.random.binomial(n, p, size=10000)
P = sum(x == k) / len(x)
print(P)

首先，我们通过np.random.binomial(n, p, size=10000)生成了一万个随机样本。然后，我们统计其中恰好有50000个正面朝上的样本个数，并除以总样本数，得到概率0.0269。这说明在抛掷硬币十万次中，硬币正面朝上恰好50000次的概率非常小。这样的模拟可以帮助我们更好地理解和预测真实的问题。

不过，当试验次数更大时，即使模拟也会非常困难。这时，我们可以使用Stirling公式和中心极限定理来近似计算。

Stirling公式和中心极限定理在计算二项式概率分布中的应用

Stirling公式是关于 $n!$ 的一个近似公式，它的精确程度随着 $n$ 的增大而提高。其中，一个常用形式的Stirling公式为：
$n! \approx \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$

中心极限定理是一种概率论中的重要定理，它表明当试验次数 $n$ 趋近于无穷大时，二项式分布可以近似地看作正态分布。即二项式分布 $B(n,p)$ 在 $n$ 很大时可以近似为：
$B(n,p) \approx N(np, np(1-p))$

其中， $N(\mu, \sigma^2)$ 表示均值为 $\mu$ ，标准差为 $\sigma$ 的正态分布。

我们来看一个例子，假设我们翻100次硬币，求硬币正面朝上恰好50次的概率。我们可以用Stirling公式和中心极限定理的近似方法计算。

import numpy as np

n = 100
p = 0.5
k = 50

P = (2 * np.pi * n * p * (1 - p)) ** (-0.5) * np.exp((-1 * (k - n * p) ** 2) / (2 * n * p * (1 - p)))
print(P)

运行结果为：