4.5 统计stats模块

NumPy库已经提供了一些基本的统计函数，如求期望、方差、中位数、最大值和最小值等。示例代码：

import numpy as np

#构建一个1000个随机变量的数组
x = np.random.randn(1000)

#对数组元素的值进行统计
mean = x.mean()
std = x.std()
var = x.var()

print(mean,std,var)

运行结果：

(0.02877273942510088, 0.97623362287515114, 0.95303208643194282)

mean是期望值，std是标准差，var是方差，使用numpy.array对象已有的方法获得统计指标快速有效，而SciPy库则提供了更高级的统计工具，它的Stats模块包含了多种概率分布的随机变量（随机变量是指概率论中的概念，不是Python中的变量），其中随机变量又分为连续和离散两种。所有的连续随机变量都是rv_continuous的派生类的对象，而所有的离散随机变量都是rv_discrete的派生类的对象。

4.5.1 连续概率分布

SciPy的stats模块提供了大约80种连续随机变量和10多种离散分布变量，这些分布都依赖于numpy.random函数。可以通过如下语句获得stats模块中所有的连续随机变量，示例代码：

from scipy import stats
[k for k, v in stats.__dict__.items() if isinstance(v, stats.rv_continuous)]

运行结果：

['ksone', 'kstwobign', 'norm', 'alpha', 'anglit', 'arcsine', 'beta', 'betaprime', 'bradford', 'burr', 'burr12', 'fisk', 'cauchy', 'chi', 'chi2', 'cosine', 'dgamma', 'dweibull', 'expon', 'exponnorm', 'exponweib', 'exponpow', 'fatiguelife', 'foldcauchy', 'f', 'foldnorm', 'frechet_r', 'weibull_min', 'frechet_l', 'weibull_max', 'genlogistic', 'genpareto', 'genexpon', 'genextreme', 'gamma', 'erlang', 'gengamma', 'genhalflogistic', 'gompertz', 'gumbel_r', 'gumbel_l', 'halfcauchy', 'halflogistic', 'halfnorm', 'hypsecant', 'gausshyper', 'invgamma', 'invgauss', 'invweibull', 'johnsonsb', 'johnsonsu', 'laplace', 'levy', 'levy_l', 'levy_stable', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'gilbrat', 'maxwell', 'mielke', 'kappa4', 'kappa3', 'nakagami', 'ncx2', 'ncf', 't', 'nct', 'pareto', 'lomax', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rdist', 'rayleigh', 'reciprocal', 'rice', 'recipinvgauss', 'semicircular', 'skewnorm', 'trapz', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'vonmises_line', 'wald', 'wrapcauchy', 'gennorm', 'halfgennorm']

连续随机变量对象主要使用如下方法，下表所示：

方法名	全称	功能
rvs	Random Variates of given type	对随机变量进行随机取值，通过size参数指定输出数组的大小
pdf	Probability Density Function	随机变量的概率密度函数
cdf	Cumulative Distribution Function	随机变量的累积分布函数，它是概率密度函数的积分
sf	Survival function	随机变量的生存函数，它的值是1-cdf(t)
ppf	Percent point function	累积分布函数的反函数
stat	statistics	计算随机变量的期望值和方差
fit	fit	对一组随机取样进行拟合，找出最适合取样数据的概率密度函数的系数

下面以标准正态分布（函数表示f(x)=(1/√2π)exp(-x^2/2)）为例，简单介绍随机变量的用法。示例代码：

from scipy import stats
# 设置正态分布参数，其中loc是期望值参数，scale是标准差参数
X = stats.norm(loc=1.0, scale=2.0)

# 计算随机变量的期望值和方差
print(X.stats())

运行结果：

(array(1.0), array(4.0))

以上代码说明，norm可以像函数一样调用，通过loc和scale参数可以指定随机变量的偏移和缩放参数。对于正态分布的随机变量来说，这两个参数相当于指定其期望值和标准差，标准差是方差的算术平方根。X的stats()方法，可以计算随机变量X分布的特征值，如期望值和方差。
此外，通过调用随机变量X的rvs()方法，可以得到包含一万次随机取样值的数组x，然后调用NumPy的mean()和var()计算此数组的均值和方差，其结果符合随机变量X的特性，示例代码：

#对随机变量取10000个值
x = X.rvs(size=10000)
print(np.mean(x), np.var(x))

运行结果

(1.0287787687588861, 3.9944276709242805)

使用fit()方法对随机取样序列x进行拟合，它返回的是与随机取样值最吻合的随机变量参数，示例代码：

#输出随机序列的期望值和标准差
print(stats.norm.fit(x))

运行结果：

(1.0287787687588861, 1.998606432223283)

在下面的例子中，计算取样值x的直方图统计以及累计分布，并与随机变量的概率密度函数和累积分布函数进行比较。示例代码：

pdf, t = np.histogram(x, bins=100, normed=True)
t = (t[:-1]+t[1:])*0.5
cdf = np.cumsum(pdf) * (t[1] - t[0])
p_error = pdf - X.pdf(t)
c_error = cdf - X.cdf(t)
print("max pdf error: {}, max cdf error: {}".format(np.abs(p_error).max(), np.abs(c_error).max()))

运行结果：

max pdf error: 0.0208405611169, max cdf error: 0.0126874590568

通过绘图的方式查看概率密度函数求得的理论值（theory value）和直方图统计值（statistic value）的结果，可以看出二者是一致的，示例代码：

import pylab as pl
pl.plot(t, pdf, color="green", label = "statistic value")
pl.plot(t, X.pdf(t), color="yellow", label ="theory value")
pl.legend(loc = "best")
pl.show()

绘图如下：

正态分布的概率密度函数

也可以用同样的方式显示随机变量X的累积分布和数组pdf的累加结果，示例代码：

import pylab as pl
pl.plot(t, cdf, color="green", label = "statistic value")
pl.plot(t, X.cdf(t), color="yellow", label ="theory value")
pl.legend(loc = "best")
pl.show()

绘图如下：

累积分布函数

4.5.2 离散概率分布

投掷有六个面的骰子时，只能获得1到6的整数，因此所得到的概率分布是离散的。我们以值域离散的分布称为离散概率分布，包括泊松分布、二项分布、几何分布等。通常使用概率质量函数（PMF）描述其分布情况，如几何分布函数PMF=(1-p)^(k-1)p。
在stats模块中所有描述离散分布的随机变量都从rv_discrete类继承，也可以直接用rv_discrete类自定义离散概率分布。假设有一个不均匀的骰子，其各点出现的概率不相等，我们用如下代码定义其分布，示例代码：

# 数组x保存骰子的所有可能值，数组p保存每个值出现的概率
x = range(1, 7)
p = (0.4, 0.2, 0.1, 0.1, 0.1, 0.1)

# 创建表示这个骰子的随机变量dice，调用其rvs()方法投掷此骰子20次，获得符合概率p的随机数
dice = stats.rv_discrete(values=(x, p))
print(dice.rvs(size=20))

运行结果：

array([3, 6, 4, 5, 5, 2, 1, 3, 3, 1, 1, 3, 1, 5, 1, 3, 4, 1, 2, 2])

除了自定义离散概率分布，我们也可以利用stats模块里的函数定义各种分布。下面以生成几何分布为例，其函数是geom()，示例代码：

import numpy as np
from scipy.stats import geom

# 设置几何分布的参数
p = 0.5
dist = geom(p)

# 设置样本区间  
x = np.linspace(0, 5, 1000)  

# 得到几何分布的 PMF 和CDF  
pmf = dist.pmf(x) 
cdf = dist.cdf(x)  

# 生成500个随机数  
sample = dist.rvs(500)

4.5.3 描述与检验函数

SciPy中有超过60种统计函数。stats模块包括了诸如kstest 和normaltest等样本测试函数，用来检测样本是否服从某种分布。在使用这些工具前，要对数据有较好的理解，否则可能会误读它们的结果。样本分布检验为例，示例代码：

import numpy as np 
from scipy import stats 

# 生成包括100个服从正态分布的随机数样本
sample = np.random.randn(100) 

# 用normaltest检验原假设
out = stats.normaltest(sample) 
print('normaltest output') 
print('Z-score = ' + str(out[0])) 
print('P-value = ' + str(out[1])) 

# kstest 是检验拟合度的Kolmogorov-Smirnov检验，这里针对正态分布进行检验
# D是KS统计量的值，越接近0越好
out = stats.kstest(sample, 'norm') 
print('\nkstest output for the Normal distribution') 
print('D = ' + str(out[0])) 
print('P-value = ' + str(out[1])) 

# 类似地可以针对其他分布进行检验，例如Wald分布
out = stats.kstest(sample, 'wald') 
print('\nkstest output for the Wald distribution') 
print('D = ' + str(out[0])) 
print('P-value = ' + str(out[1]))

SciPy的stats模块中还提供了一些描述函数，如几何平均（gmean）、偏度（skew
）、样本频数（itemfreq）等。示例代码：

import numpy as np 
from scipy import stats 

# 生成包括100个服从正态分布的随机数样本
sample = np.random.randn(100) 

# 调和平均数，样本值须大于0 
out = stats.hmean(sample[sample > 0]) 
print('Harmonic mean = ' + str(out)) 

# 计算-1到1之间样本的均值
out = stats.tmean(sample, limits=(-1, 1)) 
print('\nTrimmed mean = ' + str(out)) 

# 计算样本偏度
out = stats.skew(sample) 
print('\nSkewness = ' + str(out)) 

# 函数describe可以一次给出样本的多种描述统计结果
out = stats.describe(sample) 
print('\nSize = ' + str(out[0])) 
print('Min = ' + str(out[1][0])) 
print('Max = ' + str(out[1][1])) 
print('Mean = ' + str(out[2])) 
print('Variance = ' + str(out[3])) 
print('Skewness = ' + str(out[4])) 
print('Kurtosis = ' + str(out[5]))