python实现常见一元随机变量的概率分布

一. 随机变量

随机变量是一个从样本空间$\Omega$到实数空间$R$的函数，比如随机变量$X$可以表示投骰子的点数。随机变量一般可以分为两类：

• 离散型随机变量：随机变量的取值为有限个。
• 连续型随机变量：随机变量的取值是连续的，有无限多个。

scipy.stat模块中包含了多种概率分布的随机变量,包含离散型随机变量和连续型随机变量。离散型随机变量的常见接口如下：

连续型随机变量的常见接口如下：

二. 常见离散分布

1. 二项分布

如果随机变量$X$的分布律为 P ( X = k ) = Cnkpkq n − k， k = 0 , 1 , . . . n ，P(X=k) = C^k_np^kq^{n-k}，k = 0,1,…n，P(X=k)=Cnk​pkqn−k，k=0,1,…n，其中 $p+q=1$ ,则称$X$服从参数为$n,p$的二项分布，记为$X\sim B(n,p)$。

• 期望：$E(X)=np$
• 方差:$D(X)=np(1-p)$

1. 画出不同参数下的二项分布，$n,p$分别为$(10，0.3),（10，0.5）,（10，0.7）$

   import numpy as npfrom scipy.stats import binomimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(3, 1, figsize = (10, 10))# 调整子图间距fig.subplots_adjust(hspace = 0.5)params = [(10, 0.3), (10, 0.5), (10, 0.7)]for i in range(len(params)):n = params[i][0]p = params[i][1]x = np.arange(0, n + 1)y = binom(n, p).pmf(x)# 计算随机变量的期望,方差mean, var = binom.stats(n, p, moments='mv')ax[i].scatter(x, y, color = 'blue', marker = 'o')ax[i].set_title('n = {}, p = {}'.format(n, p))ax[i].set_xticks(x)ax[i].text(1, 0.2, '期望: {:.2f}\n方差: {:.2f}'.format(mean, var))ax[i].grid()plt.show()

   运行结果：
   

2. 生成服从不同参数二项分布的随机数组(采样100000次)，然后查看数组的频率分布

   import numpy as npfrom scipy.stats import binomimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(3, 1, figsize = (10, 10))# 调整子图间距fig.subplots_adjust(hspace = 0.5)params = [(10, 0.3), (10, 0.5), (10, 0.7)]for i in range(len(params)):n = params[i][0]p = params[i][1]x = np.arange(0, 11)# 抽样10万次sample = binom.rvs(n = n, p = p, size = 100000)print(sample)ax[i].hist(sample, color = 'blue', density=True, bins = 50)ax[i].set_title('n = {}, p = {}'.format(n, p))ax[i].set_xticks(x)ax[i].grid()plt.show()

   运行结果：
   

2. 几何分布

若随机变量$X$的分布律为 P ( X = k ) = ( 1 − p ) k − 1p ， k = 1 , 2 , . . . ，P(X = k) = (1 – p)^{k – 1}p，k = 1, 2, …，P(X=k)=(1−p)k−1p，k=1,2,…，其中 $0<p<1$,则称$X$服从参数为$p$的几何分布，记为$X\sim Ge(p)$。

• 期望：E(X)= 1 pE(X) = \frac{1}{p} E(X)=p1​
• 方差：D(X)= 1−pp2 D(X) = \frac{1 – p}{p^2} D(X)=p21−p​

1. 画出不同参数下的几何分布，$p$分别为$(0.3，0.5，0.7)$

   import numpy as npfrom scipy.stats import geomimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(3, 1, figsize = (10, 10))# 调整子图间距fig.subplots_adjust(hspace = 0.5)params = [0.3,0.5,0.7]for i in range(len(params)):p = params[i]x = np.arange(1, 15)y = geom(p = p).pmf(x)print(y)# 计算随机变量的期望,方差mean, var = geom.stats(p = p, moments='mv')ax[i].scatter(x, y, color = 'blue', marker = 'o')ax[i].set_title('p = {}'.format(p))ax[i].set_xticks(x)ax[i].text(5, 0.2, '期望: {:.2f}\n方差: {:.2f}'.format(mean, var))ax[i].grid()plt.show()

   运行结果：
   

2. 生成服从不同参数几何分布的随机数组(采样100000次)，然后查看数组的频率分布

   import numpy as npfrom scipy.stats import geomimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(3, 1, figsize = (10, 10))# 调整子图间距fig.subplots_adjust(hspace = 0.5)params = [0.3, 0.5, 0.7]for i in range(len(params)):p = params[i]x = np.arange(0, 15)# 抽样sample = geom.rvs(p = p, size = 100000)print(sample)ax[i].hist(sample, color = 'blue', density=True, bins = 50)ax[i].set_title('p = {}'.format(p))ax[i].set_xlim(0,15)ax[i].set_xticks(x)ax[i].grid()plt.show()

   运行结果：
   

3. 泊松分布

若随机变量$X$的分布律为 P ( X = k ) =λ k k !e − λ， k = 0 , 1 , 2… ，P(X=k) = \frac{\lambda^k}{k!}e^{-\lambda}，k = 0, 1, 2 …，P(X=k)=k!λk​e−λ，k=0,1,2…，其中$\lambda >0，$则称$X$服从参数为$\lambda$的泊松分布，记为$X\sim P(\lambda )$。

• 期望：$E(X)=\lambda$
• 方差：$D(X)=\lambda$

1. 画出不同参数下的泊松分布，$\lambda$分别为$(2,6,8)$

   import numpy as npfrom scipy.stats import poissonimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(3, 1, figsize = (10, 10))# 调整子图间距fig.subplots_adjust(hspace = 0.5)params = [2,6,8]for i in range(len(params)):numda = params[i]x = np.arange(1, 15)y = poisson(numda).pmf(x)# 计算随机变量的期望,方差mean, var = poisson.stats(numda, moments='mv')ax[i].scatter(x, y, color = 'blue', marker = 'o')ax[i].set_title('lambda = {}'.format(numda))ax[i].set_xticks(x)ax[i].set_yticks([0, 0.1, 0.2, 0.3, 0.4])ax[i].text(5, 0.2, '期望: {:.2f}\n方差: {:.2f}'.format(mean, var))ax[i].grid()plt.show()

   运行结果：
   

2. 生成服从不同参数泊松分布的随机数组(采样100000次)，然后查看数组的频率分布

   import numpy as npfrom scipy.stats import poissonimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(3, 1, figsize = (10, 10))# 调整子图间距fig.subplots_adjust(hspace = 0.5)params = [2, 6, 8]for i in range(len(params)):numda = params[i]x = np.arange(0, 16)# 抽样sample = poisson.rvs(numda, size = 1000000)print(sample)ax[i].hist(sample, color = 'blue', density=True, bins = 50)ax[i].set_title('lamdba = {}'.format(numda))ax[i].set_xticks(x)ax[i].set_xlim(0, 16)ax[i].grid()plt.show()

   运行结果：
   

三. 常见连续分布

1. 正太分布

若随机变量$X$的概率密度函数为 f ( x ) = 12 πδe −( x − μ )22 δ2， ( − ∞ < x < + ∞ )f(x) = \frac{1}{\sqrt{2\pi}\delta}e^{- \frac{(x – \mu)^2}{2\delta^2}}，( -\infty< x < +\infty)f(x)=2π ​δ1​e−2δ2(x−μ)2​，(−∞<x<+∞)，则称$X$服从参数为 ( μ ， δ2)(\mu，\delta^2)(μ，δ2)的正太分布，记为 X ∼ N ( μ ， δ2)X \sim N(\mu，\delta^2)X∼N(μ，δ2)。当$\mu =0，\delta =1$时称$X$服从标准正太分布。

• 期望：$E(X)=\mu$
• 方差：D(X)= δ 2D(X) = \delta^2 D(X)=δ2

1. 画出不同参数下的正太分布，$\mu ，\delta$分别为$(0,1),(0,3)$

   import numpy as npfrom scipy.stats import normimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(figsize=(10, 8))params = [(0, 1, 'red'), (0, 3, 'blue')]x = np.linspace(-20, 20, 1000)for i in range(0, len(params)):loc = params[i][0]scale = params[i][1]color = params[i][2]mean, var = norm.stats(loc, scale, moments='mv')ax.plot(x, norm(loc = loc, scale = scale).pdf(x), color = color, label = 'loc={},scale={},均值={},方差={}'.format(loc, scale,mean,var))ax.set_xticks(np.arange(-20, 21))ax.grid()ax.legend()plt.show()

2. 生成服从不同参数正太分布的随机数组(采样100000次)，然后查看数组的频率分布

   import numpy as npfrom scipy.stats import normimport matplotlib.pyplot as pltplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(2, 1, figsize=(10, 8))params = [(0, 1, 'red'), (0, 3, 'blue')]x = np.linspace(-20, 20, 1000)# 采样for i in range(0, len(params)):loc = params[i][0]scale = params[i][1]color = params[i][2]# 画出分布图ax[i].plot(x, norm(loc = loc, scale = scale).pdf(x), color = color, label = 'loc={},scale={}'.format(loc, scale))# 画出随机抽样的频率分布直方图ax[i].hist(norm(loc = loc, scale = scale).rvs(size = 100000), density=True, bins = 100)ax[i].set_xticks(np.arange(-20, 21))ax[i].grid()ax[i].legend()plt.show()

2. 指数分布

若随机变量$X$的概率密度函数为 f ( x ) = { λ e − λ xx ≥ 00 x 0 )f(x) = \begin{cases} {\lambda}e^{-{\lambda}x} & x \ge 0\\0 & x 0)f(x)={λe−λx0​x≥0x<0​(λ>0)，则称$X$服从参数为$\lambda$的指数分布，记为$X\sim E(\lambda )$。

• 期望：E(X)= 1 λE(X) = \frac{1}{\lambda} E(X)=λ1​
• 方差：D(X)= 1 λ2 D(X) = \frac{1}{{\lambda}^2} D(X)=λ21​

scipy中指数分布expon的参数传入$\lambda$的倒数。

A common parameterization for expon is in terms of the rate parameter lambda, such that pdf = lambda * exp(-lambda * x). This parameterization corresponds to using scale = 1 / lambda.

1. 画出不同参数下的指数分布，$\lambda$分别为$(0.5,1,1.5)$

   import numpy as npimport matplotlib.pyplot as pltfrom scipy.stats import exponplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(figsize = (10, 8))params = [(0.5, 'red'), (1, 'blue'), (1.5, 'green')]x = np.linspace(0, 15, 1000)for i in range(0, len(params)):numda = params[i][0]color = params[i][1]mean, var = expon.stats(loc = 0, scale = 1 / numda, moments='mv')ax.plot(x, expon(scale = 1 / numda).pdf(x), color = color, label = 'lambda = {:.2f}, 均值:{:.2f}, 方差: {:.4f}'.format(numda, mean, var))ax.grid()ax.legend()plt.show()

2. 生成服从不同参数指数分布的随机数组(采样100000次)，然后查看数组的频率分布

   import numpy as npimport matplotlib.pyplot as pltfrom scipy.stats import exponplt.rcParams["font.family"] = "SimHei"# 设置字体plt.rcParams["axes.unicode_minus"] = False# 正常显示负号if __name__ == '__main__':fig, ax = plt.subplots(3, 1, figsize = (10, 8))params = [(0.5, 'red'), (1, 'blue'), (1.5, 'green')]x = np.linspace(0, 15, 1000)# 采样for i in range(0, len(params)):numda = params[i][0]color = params[i][1]ax[i].plot(x, expon(scale = 1/numda).pdf(x), color = color, label = 'lambda={}'.format(numda))ax[i].hist(expon(scale = 1/numda).rvs(size = 10000), density=True, bins = 100)ax[i].set_xticks(np.arange(0, 15))ax[i].set_xlim(0, 15)ax[i].grid()ax[i].legend()plt.show()