常见概率分布的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%