概率统计方法简介
Python
中常用的统计工具有 Numpy, Pandas, PyMC, StatsModels
等。
Scipy
中的子库 scipy.stats
中包含很多统计上的方法。
导入 numpy
和 matplotlib
：
1%pylab inline
1Populating the interactive namespace from numpy and matplotlib
1heights = array([1.46, 1.79, 2.01, 1.75, 1.56, 1.69, 1.88, 1.76, 1.88, 1.78])
Numpy
自带简单的统计方法：
1print('均值mean: ', heights.mean())
2print('中位数median: ', np.median(heights))
3print('忽略nan值之后的中位数nanmedian: ', np.nanmedian(heights))
4print('最小值min: ', heights.min())
5print('最大值max: ', heights.max())
6print('标准差std: ', heights.std())
1均值mean: 1.7559999999999998
2中位数median: 1.77
3忽略nan值之后的中位数nanmedian: 1.77
4最小值min: 1.46
5最大值max: 2.01
6标准差std: 0.15081114017207078
导入 Scipy
的统计模块：
1import scipy.stats.stats as st
其他统计量：
1print('mode, ', st.mode(heights)) # 众数及其出现次数
2print('skewness, ', st.skew(heights)) # 偏度
3print('kurtosis, ', st.kurtosis(heights)) # 峰度
4print('and so many more...')
1mode, ModeResult(mode=array([1.88]), count=array([2]))
2skewness, -0.3935244564726347
3kurtosis, -0.33067209772439865
4and so many more...
概率分布
常见的连续概率分布有：
离散概率分布：
这些都可以在 scipy.stats
中找到。
连续分布
正态分布
以正态分布为例，先导入正态分布：
1from scipy.stats import norm
它包含四类常用的函数：
从正态分布产生500个随机点：
1x_norm = norm.rvs(size=500)
2x_norm[:5]
1array([-1.33352254, 0.0858442 , 0.29593714, 0.39278228, -0.60690144])
直方图：
1p = hist(x_norm)
2print('counts: ', p[0])
3print('bin centers:', p[1])
1counts: [ 7. 15. 42. 69. 96. 106. 92. 40. 23. 10.]
2bin centers: [-2.71399184 -2.18462072 -1.65524961 -1.12587849 -0.59650737 -0.06713626
3 0.46223486 0.99160597 1.52097709 2.05034821 2.57971932]
归一化直方图（用出现频率代替次数），将划分区间变为 20
（默认 10
）：
1h = hist(x_norm, density=True, bins=20)
在这组数据下，正态分布参数的最大似然估计值为：
1x_mean, x_std = norm.fit(x_norm)
2
3print('mean, ', x_mean)
4print('x_std, ', x_std)
1mean, 0.01748801315861445
2x_std, 0.9788390100195437
将真实的概率密度函数与直方图进行比较：
1h = hist(x_norm, density=True, bins=20)
2
3x = linspace(-3, 3, 50)
4p = plot(x, norm.pdf(x), 'r-')
导入积分函数：
1from scipy.integrate import trapz
通过积分，计算落在某个区间的概率大小：
1x1 = linspace(-2, 2, 108)
2p = trapz(norm.pdf(x1), x1)
3print(f'{p:.2%}的值落在-2到2之间')
4
5fill_between(x1, norm.pdf(x1), color='red')
6plot(x, norm.pdf(x), 'k-')
195.45%的值落在-2到2之间
2[]
默认情况，正态分布的参数为均值0，标准差1，即标准正态分布。
可以通过 loc
和 scale
来调整这些参数，一种方法是调用相关函数时进行输入：
1p = plot(x, norm.pdf(x, loc=0, scale=1), label="N(0,1)")
2p = plot(x, norm.pdf(x, loc=0.5, scale=2), label="N(0.5,2)")
3p = plot(x, norm.pdf(x, loc=-0.5, scale=.5), label="N(-0.5,0.5)")
4plt.legend()
5plt.show()
另一种则是将 loc, scale
作为参数直接输给 norm
生成相应的分布：
1p = plot(x, norm(loc=0, scale=1).pdf(x), label="N(0,1)")
2p = plot(x, norm(loc=0.5, scale=2).pdf(x), label="N(0.5,2)")
3p = plot(x, norm(loc=-0.5, scale=.5).pdf(x), label="N(-0.5,0.5)")
4plt.legend()
5plt.show()
其他连续分布
1from scipy.stats import lognorm, t, dweibull
支持与 norm
类似的操作，如概率密度函数等。
不同参数的对数正态分布：
lognorm.pdf要得到一般意义上符合对数正态分布的随机变量X（logX服从n(mu,sigma^2)），需要令lognorm中的参数s=sigma,loc=0,scale=exp(mu)。
1x = linspace(-3, 3, 300)
2
3plot(x, lognorm.pdf(x, .1), label='$\sigma$=0.1')
4plot(x, lognorm.pdf(x, .5), label='$\sigma$=0.5')
5plot(x, lognorm.pdf(x, .8), label='$\sigma$=0.8')
6plot(x, lognorm.pdf(x, 1), label='$\sigma$=1')
7
8legend()
1<matplotlib.legend.Legend at 0x174c65f8>
不同的韦氏分布：
1x = linspace(0.01, 3, 100)
2
3plot(x, dweibull.pdf(x, 1), label='s=1, constant failure rate')
4plot(x, dweibull.pdf(x, 2), label='s>1, increasing failure rate')
5plot(x, dweibull.pdf(x, .1), label='0)
6
7legend()
1<matplotlib.legend.Legend at 0x175d47b8>
不同自由度的学生 t
分布：
1x = linspace(-3, 3, 100)
2
3plot(x, t.pdf(x, 1), label='df=1')
4plot(x, t.pdf(x, 2), label='df=2')
5plot(x, t.pdf(x, 100), label='df=100')
6plot(x[::5], norm.pdf(x[::5]), 'kx', label='normal')
7
8legend()
1<matplotlib.legend.Legend at 0x1737ef60>
离散分布
导入离散分布：
1from scipy.stats import binom, poisson, randint
离散分布没有概率密度函数，但是有概率质量函数。
离散均匀分布的概率质量函数（PMF）：
1high = 10
2low = -10
3
4x = arange(low, high+1, 0.5)
5p = stem(x, randint(low, high).pmf(x), use_line_collection=True) # 杆状图
二项分布：
1num_trials = 60
2x = arange(num_trials)
3
4plot(x, binom(num_trials, 0.5).pmf(x), 'o-', label='p=0.5')
5plot(x, binom(num_trials, 0.2).pmf(x), 'o-', label='p=0.2')
6
7legend()
1<matplotlib.legend.Legend at 0x1877d668>
泊松分布：
1x = arange(0,21)
2
3plot(x, poisson(1).pmf(x), 'o-', label=r'$\lambda$=1')
4plot(x, poisson(4).pmf(x), 'o-', label=r'$\lambda$=4')
5plot(x, poisson(9).pmf(x), 'o-', label=r'$\lambda$=9')
6
7legend()
1<matplotlib.legend.Legend at 0x18813eb8>
自定义离散分布
导入要用的函数：
1from scipy.stats import rv_discrete
一个不均匀的骰子对应的离散值及其概率：
1xk = [1, 2, 3, 4, 5, 6]
2pk = [.3, .35, .25, .05, .025, .025]
定义离散分布：
1loaded = rv_discrete(values=(xk, pk))
此时， loaded
可以当作一个离散分布的模块来使用。
产生两个服从该分布的随机变量：
1loaded.rvs(size=2)
1array([1, 2])
产生100个随机变量，将直方图与概率质量函数进行比较：
1samples = loaded.rvs(size=100)
2bins = linspace(.5, 6.5, 7)
3
4hist(samples, bins=bins, density=True)
5stem(xk,
6 loaded.pmf(xk),
7 markerfmt='ro',
8 linefmt='r-',
9 use_line_collection=True)
1object of 3 artists>
假设检验
导入相关的函数：
t
检验的相关内容请参考：
独立样本T检验用于检验两个正态分布总体的均值是否相等，两个正态分布总体的方差一致与不一致时计算方式不一样，即检验假设Ho：μ1=μ2是否成立。
配对样本T检验用于检验两个相关的样本是否来自具有相同均值的正态分布。即检验假设Ho：μ=μ1-μ2=0，实质就是检验差值的均值和零均值之间差异的显著性。
1from scipy.stats import norm
2from scipy.stats import ttest_ind, ttest_rel, ttest_1samp
3from scipy.stats import t
独立样本 t 检验
两组参数不同的正态分布：
1n1 = norm(loc=0.3, scale=1.0)
2n2 = norm(loc=0, scale=1.0)
从分布中产生两组随机样本：
1n1_samples = n1.rvs(size=100)
2n2_samples = n2.rvs(size=100)
将两组样本混合在一起：
1samples = hstack((n1_samples, n2_samples))
最大似然参数估计：
1loc, scale = norm.fit(samples)
2n = norm(loc=loc, scale=scale)
比较：
1x = linspace(-3, 3, 100)
2
3hist([samples, n1_samples, n2_samples],
4 color=['royalblue', 'orange', 'seagreen'],
5 density=True)
6plot(x, n.pdf(x), '-', color="royalblue", label="mix")
7plot(x, n1.pdf(x), '-', color="orange", label="loc=0.3")
8plot(x, n2.pdf(x), '-', color="seagreen", label="loc=0")
9plt.legend()
10plt.show()
独立双样本 t
检验的目的在于判断两组样本之间均值是否有显著差异：
1t_val, p = ttest_ind(n1_samples, n2_samples)
2
3print(f't = {t_val}')
4print(f'p-value = {p}')
1t = 2.1753813477310393
2p-value = 0.030786184428796288
p
值小，两组样本之间对应的正态分布总体的均值相等的概率大于5%，有显著性差异。
配对样本 t 检验
配对样本指的是两组样本之间的元素一一对应，例如，假设我们有一组病人的数据：
1pop_size = 35
2
3pre_treat = norm(loc=0, scale=1)
4n0 = pre_treat.rvs(size=pop_size)
经过某种治疗后，对这组病人得到一组新的数据：
1effect = norm(loc=0.05, scale=0.2)
2eff = effect.rvs(size=pop_size)
3
4n1 = n0 + eff
新数据的最大似然估计：
1loc, scale = norm.fit(n1)
2post_treat = norm(loc=loc, scale=scale)
画图：
1fig = figure(figsize=(10, 4))
2
3ax1 = fig.add_subplot(1, 2, 1)
4h = ax1.hist([n0, n1], color=['dodgerblue', 'orange'], density=True)
5x = linspace(-3, 3, 100)
6p = ax1.plot(x, pre_treat.pdf(x), color='dodgerblue', linefont-size: inherit;line-height: inherit;color: rgb(162, 252, 162);overflow-wrap: inherit !important;word-break: inherit !important;">'-', label='pre_treat')
7p = ax1.plot(x, post_treat.pdf(x), color='orange', linefont-size: inherit;line-height: inherit;color: rgb(162, 252, 162);overflow-wrap: inherit !important;word-break: inherit !important;">'-', label='post_treat')
8plt.legend()
9
10ax2 = fig.add_subplot(1, 2, 2)
11h = ax2.hist(eff, color='dodgerblue', density=True)
12x = linspace(-0.5, 0.5, 100)
13p = ax2.plot(x, effect.pdf(x), color='dodgerblue', linefont-size: inherit;line-height: inherit;color: rgb(162, 252, 162);overflow-wrap: inherit !important;word-break: inherit !important;">'-', label='effect')
14plt.legend()
1<matplotlib.legend.Legend at 0x18997fd0>
独立 t
检验：
1t_val, p = ttest_ind(n0, n1)
2
3print(f't = {t_val}')
4print(f'p-value = {p}')
1t = -0.23922148552080355
2p-value = 0.811653400430429
高 p
值说明两组样本之间对应的正态分布总体的均值相等的概率大于5%，没有显著性差异。
配对 t
检验：
1t_val, p = ttest_rel(n0, n1)
2
3print(f't = {t_val}')
4print(f'p-value = {p}')
1t = -1.3183846640746473
2p-value = 0.19618772919021144
配对 t
检验的结果说明，两组样本之间均值无差异的概率大于5%，存在显著性差异，说明治疗效果没有达到预期。
p 值计算原理
p
值对应的部分是下图中的红色区域，边界范围由 t
值决定。
1my_t = t(pop_size) # 传入参数为自由度，这里自由度为50
2
3x = linspace(-3, 3, 300)
4p = plot(x, my_t.pdf(x), 'b-')
5lower_x = x[x<= -abs(t_val)]
6upper_x = x[x>= abs(t_val)]
7
8p = fill_between(lower_x, my_t.pdf(lower_x), color='red')
9p = fill_between(upper_x, my_t.pdf(upper_x), color='red')
京公网安备 11010802041100号