常见概率分布的Python实现

更新于
概率分布，是概率论的基本概念之一，包括离散概率分布和连续概率分布，主要用以表述随机变量取值的概率规律。前面讲的抽样分布属于连续概率分布的范畴。

概率分布，是概率论的基本概念之一，包括离散概率分布和连续概率分布，主要用以表述随机变量取值的概率规律。前面讲的抽样分布属于连续概率分布的范畴。

离散概率分布

离散概率分布也称为概率质量函数（probability mass function），例子有伯努利分布、二项分布、泊松分布、几何分布、分类分布、多项式分布等。

1、伯努利分布（0-1分布）

伯努利分布（Bernoulli Distribution）又名0-1分布、两点分布。
介绍伯努利分布前首先需要引入伯努利试验（Bernoulli Trial）。伯努利试验是只有两种可能结果的单次随机试验，即对于一个随机变量X而言，伯努利试验都可以表达为”是或否”的问题。例如，抛一次硬币是正面向上吗？刚出生的小孩是个女孩吗？等等。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# coding:utf-8
'''
伯努利分布
'''
import random
import numpy as np
from matplotlib import pyplot as plt

def bernoulli(p, k):
    return p if k else 1 - p

n_experiment = 10
p = 0.6
x = np.arange(n_experiment)
y = []
for _ in range(n_experiment):
    pick = bernoulli(p, k=bool(random.getrandbits(1)))
    y.append(pick)

u, s = np.mean(y), np.std(y)
plt.scatter(x, y, label='mu=%.2f,sigma=%.2f' % (u, s))
plt.legend()
plt.show()

2、二项分布

如果某试验是一个伯努利试验，将该试验独立重复地进行n次，则称这一串重复的独立试验为n重伯努利试验。二项分布（Binomial Distribution）是n重伯努利试验成功k次的概率分布。
二项分布的典型例子是扔硬币，硬币正面朝上概率为p，重复扔n次硬币，k次为正面的概率即为一个二项分布概率。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# coding:utf-8
'''
二项分布
'''
import numpy as np
from matplotlib import pyplot as plt

import operator as op
from functools import reduce

def const(n, k):
    k = min(k, n-k)
    numer = reduce(op.mul, range(n, n-k, -1), 1)
    denom = reduce(op.mul, range(1, k+1), 1)
    return numer / denom

def binomial(n, p):
    q = 1 - p
    y = [const(n, k) * (p ** k) * (q ** (n-k)) for k in range(n)]
    return y, np.mean(y), np.std(y)

n_experiment = 50
p = 0.5
x = np.arange(n_experiment)
y, u, s = binomial(n_experiment, p)
plt.scatter(x, y, label='mu=%.2f,sigma=%.2f' % (u, s))
plt.legend()
plt.show()

3、泊松分布

泊松分布 （Poisson Distribution）是由二项分布推导来的，是二项分布的极限形态 。泊松分布中事件发生的概率非常小，不发生的概率非常大。
在实际事例中，当一个随机事件，例如某电话交换台收到的呼叫、来到某公共汽车站的乘客、某放射性物质发射出的粒子、显微镜下某区域中的白血球等等，以固定的平均瞬时速率λ（或称密度）随机且独立地出现时，那么这个事件在单位时间（面积或体积）内出现的次数或个数就近似地服从泊松分布P（λ）。因此，泊松分布在管理科学、运筹学以及自然科学的某些问题中都占有重要的地位。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# coding:utf-8
'''
泊松分布
'''
import random
import numpy as np
from matplotlib import pyplot as plt

import operator as op
from functools import reduce

def poisson(lamb, k): 
    if lamda <= 0:  # 参数合法判断
        return -1
    numer = lamb**k * np.exp(-lamb)
    denom = reduce(op.mul, range(1, k+1), 1)
    return  numer / denom

n_experiment = 30
p = 2
x = np.arange(n_experiment)
y = []
for k in range(n_experiment):
    pick = poisson(p, k)
    y.append(pick)

u, s = np.mean(y), np.std(y)
plt.scatter(x, y, label='mu=%.2f,sigma=%.2f' % (u, s))
plt.legend()
plt.show()

4、几何分布

几何分布（Geometric Distribution）是指在n重伯努利实验中，第k次实验才得到第一次成功的概率分布。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# coding:utf-8
'''
几何分布
'''
import random
import numpy as np
from matplotlib import pyplot as plt

def geometric(p, k):
    return (1 - p)**(k-1) * p

n_experiment = 100
p = 0.2
x = np.arange(n_experiment)
y = []
for k in range(n_experiment):
    pick = geometric(p, k)
    y.append(pick)

u, s = np.mean(y), np.std(y)
plt.scatter(x, y, label='mu=%.2f,sigma=%.2f' % (u, s))
plt.legend()
plt.show()

5、分类分布

分类分布（Categorical Distribution）是伯努利分布的泛化，单次随机试验的可能结果从两种拓展到m种，m个结果发生的概率互斥且和为1。
典型的例子是扔骰子，不同于扔硬币，骰子有6个面对应6个不同的点数，每个点数朝上的概率都是1/6（如果是不规则的骰子，各个面的概率可能不同）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# coding:utf-8
'''
分类分布
'''
import random
import numpy as np
from matplotlib import pyplot as plt

def categorical(p, k):
    return p[k]

n_experiment = 100
p = [0.2, 0.1, 0.7]
x = np.arange(n_experiment)
y = []
for _ in range(n_experiment):
    pick = categorical(p, k=random.randint(0, len(p) - 1))
    y.append(pick)

u, s = np.mean(y), np.std(y)
plt.scatter(x, y, label='mu=%.2f, sigma=%.2f' % (u, s))
plt.legend()
plt.show()

6、多项式分布

多项式分布（Multinomial Distribution）是二项式分布在类别上的泛化。做n次试验，每次试验的结果可以有m个且m个结果发生的概率互斥且和为1，则发生其中一个结果k次的概率就是多项式分布。
还是拿扔骰子举例，重复扔n次，k次都是点数6朝上的概率就是一个多项式分布。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
# coding:utf-8
'''
多项式分布
'''
import numpy as np
from matplotlib import pyplot as plt

import operator as op
from functools import reduce

def factorial(n):
    return reduce(op.mul, range(1, n + 1), 1)

def const(n, a, b, c):
    assert  a + b + c == n

    numer = factorial(n)
    denom = factorial(a) * factorial(b) * factorial(c)
    return numer / denom

def multinomial(n):
    ls = []
    for i in range(1, n + 1):
        for j in range(i, n + 1):
            for k in range(j, n + 1):
                if i + j + k == n:
                    ls.append([i, j, k])

    y = [const(n, l[0], l[1], l[2]) for l in ls]
    x = np.arange(len(y))
    return x, y, np.mean(y), np.std(y)

for n_experiment in [20, 21, 22]:
    x, y, u, s = multinomial(n_experiment)
    plt.scatter(x, y, label=r'$trial=%d$' % (n_experiment))

plt.legend()
plt.show()

连续概率分布

连续概率分布也称为概率密度函数（probability density function），例如正态分布、伽马分布、指数分布、卡方分布、学生t-分布、F分布、贝塔分布、均匀分布等。

1、正态分布

正态分布（Normal Distribution）又名高斯分布（Gaussian Distribution），是一个在数学、物理及工程等领域都非常重要的概率分布，正态分布的两个参数μ和σ决定了正态分布的位置和形态。期望μ=0，方差σ=1的正态分布是标准正态分布。
在社会经济问题中，有许多随机变量的概率分布都服从正态分布。例如，某地区同年龄组儿童的发育特征，如身高、体重、肺活量；某公司年销售量等等。
从同一分布大量取样然后相加，样本的和遵循（近似的）正态分布。取样数越大，样本之和就约接近正态分布。无论原分布是何种分布，这一点均成立，真是令人惊奇，这种现象被称为中心极限定理。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# coding:utf-8
'''
正态分布
'''
import numpy as np
from matplotlib import pyplot as plt

def normal(x, n):
    u = x.mean()
    s = x.std()

    x = (x - u) / s

    x = np.linspace(x.min(), x.max(), n)

    a = ((x - 0) ** 2) / (2 * (1 ** 2))
    y = 1 / (s * np.sqrt(2 * np.pi)) * np.exp(-a)

    return x, y, x.mean(), x.std()

x = np.arange(-100, 100)
x, y, u, s = normal(x, 10000)

plt.plot(x, y, label='mu=%.2f,sigma=%.2f' % (u, s))
plt.legend()
plt.show()

2、伽马分布

伽马分布（Gamma Distribution）也是一种统计学中的常见连续型分布，解决的问题是”要等到n个随机事件都发生，需要经历多久时间”。
伽马分布有两个参数，参数α称为形状参数（shape parameter），参数β称为尺度参数（scale parameter）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# coding:utf-8
'''
伽马分布
'''
import numpy as np
from matplotlib import pyplot as plt

def gamma_function(n):
    cal = 1
    for i in range(2, n):
        cal *= i
    return cal

def gamma(x, a, b):
    c = (b ** a) / gamma_function(a)
    y = c * (x ** (a - 1)) * np.exp(-b * x)
    return x, y, np.mean(y), np.std(y)

for ls in [(1, 1), (2, 1), (3, 1), (2, 2)]:
    a, b = ls[0], ls[1]
    x = np.arange(0, 20, 0.01, dtype=np.float)
    x, y, u, s = gamma(x, a=a, b=b)
    plt.plot(x, y, label='mu=%.2f, sigma=%.2f, alpha=%d, beta=%d' % (u, s, a, b))
plt.legend()
plt.show()

3、指数分布

指数分布（Exponential Distribution）是伽马分布的特殊形态，形状参数α=1时的伽马分布就是指数分布。它解决的问题是”要等到一个随机事件发生，需要经历多久时间”，比如旅客进入机场的时间间隔、打进客服中心电话的时间间隔等等。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# coding:utf-8
'''
指数分布
'''
import numpy as np
from matplotlib import pyplot as plt

def exponential(x, lamb):
    y = lamb * np.exp(-lamb * x)
    return x, y, np.mean(y), np.std(y)

lamb = 0.5
x = np.arange(0, 20, 0.01)
x, y, u, s = exponential(x, lamb=lamb)
plt.plot(x, y, label='mu=%.2f, sigma=%.2f, lambda=%d' % (u, s, lamb))
plt.legend()
plt.show()

4、卡方分布

卡方分布（chi-squared Distribution）是卡方检验的基础，它是伽马分布的特殊形态（形状参数α=n/2，尺度参数β=1/2时的伽马分布就是自由度为n的卡方分布），是统计推断中应用最为广泛的概率分布之一，例如假设检验和置信区间的计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np
from matplotlib import pyplot as plt

def gamma_function(n):
    cal = 1
    for i in range(2, n):
        cal *= i
    return cal

def chi_squared(x, k):

    c = 1 / (2 ** (k/2)) * gamma_function(k//2)
    y = c * (x ** (k/2 - 1)) * np.exp(-x /2)

    return x, y, np.mean(y), np.std(y)

for k in [2, 3, 4, 6]:
    x = np.arange(0, 10, 0.01, dtype=np.float)
    x, y, _, _ = chi_squared(x, k)
    plt.plot(x, y, label='k=%d' % (k))

plt.legend()
plt.show()

5、学生t-分布

t分布（t-distribution）是t检验的基础，与正态分布相比多了自由度参数（随着自由度参数的增加而更加接近正态分布），用于根据小样本（样本容量小于30或小于50时）来估计呈正态分布且方差未知的总体的均值。
可以说t分布是正态分布的小样本形态，如果某变量服从正态分布，当样本容量小于30或小于50时，该变量呈t分布。
t分布也是对称的倒钟型分布，就如同正态分布一样，但它的长尾占比更多，这意味着t分布更容易产生远离均值的样本。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
from matplotlib import pyplot as plt

def gamma_function(n):
    cal = 1
    for i in range(2, n):
        cal *= i
    return cal

def student_t(x, m):

    c = gamma_function((m + 1) // 2) / np.sqrt(m * np.pi) * gamma_function(m // 2)
    y = c * (1 + x**2 / m) ** (-((m + 1) / 2))

    return x, y, np.mean(y), np.std(y)

for freedom in [1, 2, 5]:

    x = np.arange(-10, 10, 0.01)
    x, y, _, _ = student_t(x, m=freedom)
    plt.plot(x, y, label='v=%d' % (freedom))

plt.legend()
plt.show()

7、F分布

F分布（F Distribution）是在t分布的基础上引申而来，为了纪念统计学家费希尔（Fisher）而用他的首字母命名，是F检验的基础。如果说t分布是正态分布的儿子，那么F分布就是正态分布的孙子。
需要注意的是F分布有两个自由度且位置不可互换。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# coding:utf-8
'''
F分布
'''
import numpy as np
from matplotlib import pyplot as plt

def gamma_function(n):
    cal = 1
    for i in range(2, n):
        cal *= i
    return cal

def f(x, m, n):

    c = gamma_function((n + m) // 2) / gamma_function(n // 2) * gamma_function(m // 2)
    y = c * (m**(m // 2) * n**(n // 2) * x**(m//2 -1) * (n + m*x)**(-((n + m) / 2)))

    return x, y, np.mean(y), np.std(y)

for freedom in [1, 3, 10, 50]:

    x = np.arange(0, 10, 0.01)
    # 给定 m = 10, n 取不同值时看分布的形状，可以看到曲线是偏态的，n 越小偏态越严重
    x, y, _, _ = f(x, m=10, n=freedom)
    plt.plot(x, y, label='v=%d' % (freedom))

plt.legend()
plt.show()

8、贝塔分布

贝塔分布（Beta Distribution）是伽马分布的一种特殊情况，它有两个参数分别代表某件事成功与失败的次数，Gamma(a,1) / Gamma(a,1) + Gamma(b,1) 就等价于 Beta(a, b) 。
贝塔分布可以看作一个概率的概率分布，当你不知道一个东西的具体概率是多少时，它可以给出了所有概率出现的可能性大小。
举一个简单的例子，熟悉棒球运动的都知道有一个指标就是棒球击球率（batting average），就是用一个运动员击中的球数除以击球的总数，我们一般认为0.266是正常水平的击球率，而如果击球率高达0.3就被认为是非常优秀的。现在有一个棒球运动员，我们希望能够预测他在这一赛季中的棒球击球率是多少。
传统的频率学派会直接计算棒球击球率，用击中的数除以击球数，但是如果这个棒球运动员只打了一次，而且还命中了，那么他就击球率就是100%了，这显然是不合理的，因为根据棒球的历史信息，我们知道这个击球率应该是0.215到0.36之间才对。对于这个问题，我们可以用一个二项分布表示（一系列成功或失败），一个最好的方法来表示这些经验（在统计中称为先验信息）就是用贝塔分布，这表示在我们没有看到这个运动员打球之前，我们就有了一个大概的范围。
贝塔分布的定义域是（0,1）这就跟概率的范围是一样的。接下来我们将这些先验信息转换为贝塔分布的参数，我们知道一个击球率应该是平均0.27左右，而他的范围是0.21到0.35，那么根据这个信息，我们可以取α=81，β=219作为贝塔分布的参数。
那么有了先验信息后，现在我们考虑一个运动员只打一次球，那么他现在的数据就是”1中;1击”。这时候我们就可以更新我们的分布了，让这个曲线做一些移动去适应我们的新信息。在这一例子里， α增加了1（击中了一次）。 β没有增加（没有漏球）。这就是我们的新的贝塔分布：Beta（81+1,219）。
只打了1次球并不能说明什么问题，但是如果我们得到了更多的数据……假设一共打了300次，其中击中了100次，200次没击中，那么这一新分布就是：Beta（81+100,219+200）
一个有趣的事情是，根据这个新的贝塔分布，我们可以得出他的数学期望为：α / （α+β）= 0.302，这一结果要比直接的估计要小100 / （100+200）= 0.333。
因此，对于一个我们不知道概率是什么，而又有一些合理的猜测时，贝塔分布能很好的作为一个表示概率的概率分布。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
# coding:utf-8
'''
贝塔分布
'''
import numpy as np
from matplotlib import pyplot as plt

def gamma_function(n):
    cal = 1
    for i in range(2, n):
        cal *= i
    return cal

def beta(x, a, b):

    gamma = gamma_function(a + b) / (gamma_function(a) * gamma_function(b))
    y = gamma * (x ** (a - 1)) * ((1 - x) ** (b - 1))
    return x, y, np.mean(y), np.std(y)

for ls in [(1, 3), (5, 1), (2, 2), (2, 5)]:
    a, b = ls[0], ls[1]

    x = np.arange(0, 1, 0.001, dtype=np.float)
    x, y, u, s = beta(x, a=a, b=b)
    plt.plot(x, y, label=r'mu=%.2f,sigma=%.2f,alpha=%d,beta=%d' % (u, s, a, b))
plt.legend()
plt.show()

9、均匀分布

均匀分布（Uniform Distribution）是贝塔分布的特殊形态，可以用Beta（1,1）来表示均匀分布。均匀分布由两个参数a和b定义，它们是数轴上的最小值和最大值，通常缩写为U（a,b）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# coding:utf-8
'''
均匀分布
'''
import numpy as np
from matplotlib import pyplot as plt

def uniform(x, a, b):

    y = [1 / (b - a) if a <= val and val <= b
                    else 0 for val in x]

    return x, y, np.mean(y), np.std(y)

x = np.arange(-100, 100)
ls = (-50, 50)
a, b = ls[0], ls[1]
x, y, u, s = uniform(x, a, b)
plt.plot(x, y, label='mu=%.2f,sigma=%.2f' % (u, s))
plt.legend()
plt.show()
0%