SciPy的stats模块包含了多种概率分布的随机变量（随机变量是指概率论中的概念，不是Python中的变量），随机变量分为连续和离散两种。所有的连续随机变量都是rv_continuous的派生类的对象，而所有的离散随机变量都是rv_discrete的派生类的对象。

一、连续概率分布

可以使用下面的语句获得stats模块中所有的连续随机变量：

连续随机变量对象都有如下方法：

·rvs：对随机变量进行随机取值，可以通过size参数指定输出的数组的大小。

·pdf：随机变量的概率密度函数。

·cdf：随机变量的累积分布函数，它是概率密度函数的积分。

·sf：随机变量的生存函数，它的值是1-cdf(t)。

·ppf：累积分布函数的反函数。

·stat：计算随机变量的期望值和方差。

·fit：对一组随机取样进行拟合，找出最适合取样数据的概率密度函数的系数。

下面以正规分布为例，简单地介绍随机变量的用法。下面的语句获得默认正规分布的随机变量的期望值和方差，我们看到默认情况下它是一个均值为0、方差为1的随机变量：

stats.norm.stats( )
(array(0.0), array(1.0))

norm可以像函数一样调用，通过loc和scale参数可以指定随机变量的偏移和缩放参数。对于正态分布的随机变量来说，这两个参数相当于指定其期望值和标准差，标准差是方差的算术平方根，因此标准差为2.0时，方差为4.0：

X = stats.norm(loc=1.0, scale=2.0)
X.stats( )
(array(1.0), array(4.0))

下面调用随机变量X的rvs( )方法，得到包含一万次随机取样值的数组x，然后调用NumPy的mean( )和var( )计算此数组的均值和方差，其结果符合随机变量X的特性：

x = X.rvs(size=10000) # 对随机变量取10000个值
np.mean(x), np.var(x) # 期望值和方差
(1.0043406567303883, 3.8899572813426553)

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

stats.norm.fit(x) # 得到随机序列的期望值和标准差
(1.0043406567303883, 1.9722974626923433)

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

?其中histogram( )对数组x进行直方图统计，它将数组x的取值范围分为100个区间，并统计x中的每个值落入各个区间的次数。histogram( )返回两个数组pdf和t，其中pdf表示各个区间的取样值出现的频数，由于normed参数为True，因此pdf的值是正规化之后的结果，其结果应与随机变量的概率密度函数一致。

?t表示区间，由于其中包括区间的起点和终点，因此t的长度为101。计算每个区间的中间值，然后调用X.pdf(t)和X.cdf(t)计算随机变量的概率密度函数和累积分布函数的值，并与统计值比较。?计算样本的累积分布时，需要与区间的大小相乘，这样才能保证其结果与累积分布函数相同。

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.0217211429624, max cdf error: 0.0209887986472

图1（左）显示了概率密度函数和直方图统计的结果，可以看出二者是一致的。右图显示了随机变量X的累积分布函数和数组pdf的累加结果。

图1　正态分布的概率密度函数（左）和累积分布函数（右）

有些随机分布除了loc和scale参数之外，还需要额外的形状参数。例如伽玛分布可用于描述等待k个独立随机事件发生所需的时间，k就是伽玛分布的形状参数。下面计算形状参数k为1和2时的伽玛分布的期望值和方差：

print stats.gamma.stats(1.0)
print stats.gamma.stats(2.0)
(array(1.0), array(1.0))
(array(2.0), array(2.0))

伽玛分布的尺度参数θ和随机事件发生的频率相关，由scale参数指定：

stats.gamma.stats(2.0, scale=2)
(array(4.0), array(8.0))

根据伽玛分布的数学定义可知其期望值为kθ，方差为kθ2。上面的程序验证了这两个公式。

当随机分布有额外的形状参数时，它所对应的rvs( )、pdf( )等方法都会增加额外的参数来接收形状参数。例如下面的程序调用rvs( )对k=2、θ=2的伽玛分布取4个随机值：

x = stats.gamma.rvs(2, scale=2, size=4)
x
array([ 2.47613445,  1.93667652,  0.85723572,  9.49088092])

接下来调用pdf( )查看与上面4个随机值对应的概率密度：

stats.gamma.pdf(x, 2, scale=2)
array([ 0.17948513,  0.18384555,  0.13960273,  0.02062186])

也可以先创建将形状参数和尺度参数固定的随机变量，然后再调用其pdf( )计算概率密度：

X = stats.gamma(2, scale=2)
X.pdf(x)
array([ 0.17948513,  0.18384555,  0.13960273,  0.02062186])

二、离散概率分布

当分布函数的值域为离散时我们称之为离散概率分布。例如投掷有六个面的骰子时，只能获得1到6的整数，因此所得到的概率分布为离散的。对于离散随机分布，通常使用概率质量函数(PMF)描述其分布情况。

在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))
dice.rvs(size=20)
array([1, 6, 3, 1, 2, 2, 4, 1, 1, 1, 2, 5, 6, 2, 4, 2, 5, 2, 1, 4])

下面我们用程序验证概率论中的中心极限定理：大量相互独立的随机变量，其均值的分布以正态分布为极限。我们计算上面那个特殊骰子投掷50次的平均值，由于每次投掷骰子都可以看作一个独立的随机事件，因此投掷50次的平均值可以看作“大量相互独立的随机变量”，其平均值的分布应该十分接近正态分布。仍然通过rvs( )获得取样值，其结果是一个形状为(20000, 50)的数组，沿着第一轴计算每行的平均值，得到samples_mean：

np.random.seed(42)
samples = dice.rvs(size=(20000, 50))
samples_mean = np.mean(samples, axis=1)

三、核密度估计

在上面的例子中，由于每个取样都是离散的，因此其平均值也是离散的，对这样的数据进行直方图统计很容易出现许多离散点恰巧聚集到同一区间的现象。为了更平滑地显示样本的概率密度，可以使用kde.gaussian_kde( )进行核密度估计。在图2中，直方图统计的结果有很大的起伏，而核密度估计与拟合的正态分布十分接近，因此验证了中心极限定理。

_, bins, step = pl.hist(
samples_mean, bins=100, normed=True, histtype="step", label=u"直方图统计")
kde = stats.kde.gaussian_kde(samples_mean)
x = np.linspace(bins[0], bins[-1], 100)
pl.plot(x, kde(x), label=u"核密度估计")
mean, std = stats.norm.fit(samples_mean)
pl.plot(x, stats.norm(mean, std).pdf(x), alpha=0.8, label=u"正态分布拟合")
pl.legend( )

图2　核密度估计能更准确地表示随机变量的概率密度函数

核密度估计算法在每个数据点处放置一条核函数曲线，最终的核密度估计就是所有这些核函数曲线的叠加。gaussian_kde( )的核函数为高斯曲线，其bw_method参数决定核函数的宽度，即高斯曲线的方差。bw_method参数可以是如下几种情况：

·当为'scott'、'silverman'时将采用相应的公式根据数据个数和维数决定核函数的宽度系数。

·当为函数时将调用此函数计算曲线宽度系数，函数的参数为gaussian_kde对象。

·当为数值时，将直接使用此数值作为宽度系数。

核函数的方差由数据的方差和宽度系数决定。

下面的程序比较宽度系数对核密度估计的影响。当宽度系数较小时，可以看到在三个数据点处的高斯曲线的峰值，而当宽度逐渐变大时，这些峰值就合并成一个统一的峰值了。

for bw in [0.2, 0.3, 0.6, 1.0]:
kde = stats.gaussian_kde([-1, 0, 1], bw_method=bw)
x = np.linspace(-5, 5, 1000)
y = kde(x)
pl.plot(x, y, lw=2, label="bw={}".format(bw), alpha=0.6)
pl.legend(loc="best")

图3显示bw_method参数越大，核密度估计曲线越平滑。

图3　bw_method参数越大，核密度估计曲线越平滑

四、二项分布、泊松分布、伽玛分布

本节用几个实例程序对概率论中的二项分布、泊松分布以及伽玛分布进行一些实验和讨论。

二项分布是最重要的离散概率分布之一。假设有一种只有两个结果的试验，其成功概率为p，那么二项分布描述了进行n次这样的独立试验，成功k次的概率。二项分布的概率质量函数公式如下：

例如，可以通过二项分布的概率质量公式计算投掷5次骰子出现3次6点的概率。投掷一次骰子，点数为6的概率（即试验成功的概率）为p=1/6，试验次数为n=5。使用二项分布的概率质量函数pmf( )可以很容易计算出现k次6点的概率。和概率密度函数pdf( )类似，pmf( )的第一个参数为随机变量的取值，后面的参数为描述随机分布所需的参数。对于二项分布来说，参数分别为n和p，而取值范围则为0到n之间的整数。下面的程序计算k为0到6时对应的概率：

stats.binom.pmf(range(6), 5, 1/6.0)
array([  4.01877572e-01,   4.01877572e-01,   1.60751029e-01,
3.21502058e-02,   3.21502058e-03,   1.28600823e-04])

由结果可知：出现0或1次6点的概率为40.2%，而出现3次6点的概率为3.215%。

在二项分布中，如果试验次数n很大，而每次试验成功的概率p很小，乘积np比较适中，那么试验成功次数的概率可以用泊松分布近似描述。

在泊松分布中使用λ描述单位时间（或单位面积）中随机事件的平均发生率。如果将二项分布中的试验次数n看作单位时间中所做的试验次数，那么它和事件出现的概率p的乘积就是事件的平均发生率λ，即λ=n·p。泊松分布的概率质量函数公式如下：

下面的程序分别计算二项分布和泊松分布的概率质量函数，结果如图4所示。可以看出当n足够大时，二者是十分接近的。程序中的事件平均发生率λ恒等于10。根据二项分布的试验次数n，计算每次事件出现的概率p=λ/n。

lambda_ = 10.0
x = np.arange(20)

n1, n2 = 100, 1000

y_binom_n1 = stats.binom.pmf(x, n1, lambda_ / n1)
y_binom_n2 = stats.binom.pmf(x, n2, lambda_ / n2)
y_poisson = stats.poisson.pmf(x, lambda_)
print np.max(np.abs(y_binom_n1 - possion))
print np.max(np.abs(y_binom_n2 - possion))
0.00675531110335
0.000630175404978

图4　当n足够大时二项分布和泊松分布近似相等

泊松分布适合描述单位时间内随机事件发生的次数的分布情况。例如某个设施在一定时间内的使用次数、机器出现故障的次数、自然灾害发生的次数等。

为了加深读者对泊松分布概念的理解，下面我们使用随机数模拟泊松分布，并与概率质量函数进行比较，结果如图5所示。图中，每秒内事件的平均发生次数为10，即λ=10。其中左图的观察时间为1000秒，而右图的观察时间为50000秒。可以看出观察时间越长，每秒内事件发生的次数越符合泊松分布。

np.random.seed(42)

def sim_poisson(lambda_, time):
t = np.random.uniform(0, time, size=lambda_ * time)?
count, time_edges = np.histogram(t, bins=time, range=(0, time))  ?
dist, count_edges = np.histogram(count, bins=20, range=(0, 20), density=True) ?
x = count_edges[:-1]
poisson = stats.poisson.pmf(x, lambda_)
return x, poisson, dist

lambda_ = 10      
times = 1000, 50000
x1, poisson1, dist1 = sim_poisson(lambda_, times[0])
x2, poisson2, dist2 = sim_poisson(lambda_, times[1])
max_error1 = np.max(np.abs(dist1 - poisson1))
max_error2 = np.max(np.abs(dist2 - poisson2))         
print "time={}, max_error={}".format(times[0], max_error1)
print "time={}, max_error={}".format(times[1], max_error2)
time=1000, max_error=0.019642302016
time=50000, max_error=0.00179801289496

图5　模拟泊松分布

?可以用NumPy的随机数生成函数uniform( )产生平均分布于0到time之间的lambda_*time个事件所发生的时刻。?用histogram( )可以统计数组t中每秒之内的事件发生的次数count，根据泊松分布的定义，count数组中的数值的分布情况应该符合泊松分布。?接下来统计事件次数在0到20区间内的概率分布。当histogram( )的density参数为True时，结果和概率质量函数相等。

还可以换个角度看随机事件的分布问题。我们可以观察相邻两个事件之间的时间间隔的分布情况，或者隔k个事件的时间间隔的分布情况。根据概率论，事件之间的时间间隔应符合伽玛分布，由于时间间隔可以是任意数值，因此伽玛分布是连续概率分布。伽玛分布的概率密度函数公式如下，它描述第k个事件发生所需的等待时间的概率分布。Γ(k)是伽玛函数，当k为整数时，它的值和k的阶乘k!相等。

下面的程序模拟了事件的时间间隔的伽玛分布，结果如图6所示。图中的观察时间为1000秒，平均每秒产生10个事件。左图中k=1，它表示相邻两个事件间的间隔的分布，而k=2则表示相隔一个事件的两个事件间的间隔的分布，可以看出它们都符合伽玛分布。

def sim_gamma(lambda_, time, k):
t = np.random.uniform(0, time, size=lambda_ * time) ?
t.sort( )  ?
interval = t[k:] - t[:-k] ?
dist, interval_edges = np.histogram(interval, bins=100, density=True) ?
x = (interval_edges[1:] + interval_edges[:-1])/2  ?
gamma = stats.gamma.pdf(x, k, scale=1.0/lambda_) ?
return x, gamma, dist

lambda_ = 10.0
time = 1000
ks = 1, 2
x1, gamma1, dist1 = sim_gamma(lambda_, time, ks[0])
x2, gamma2, dist2 = sim_gamma(lambda_, time, ks[1])

图6　模拟伽玛分布

?首先在1000秒之内产生10000个随机事件发生的时刻。因此事件的平均发生次数为每秒10次。?为了计算事件前后的时间间隔，需要先对随机时刻进行排序，?然后再计算k个事件之间的时间间隔。?对该时间间隔调用histogram( )进行概率统计，设置density为True可以直接计算概率密度。histogram( )返回的第二个值为统计区间的边界，?接下来用gamma.pdf( )计算伽玛分布的概率密度时，使用各个区间的中值进行计算。pdf( )的第二个参数为k值，scale参数为1/λ。

接下来我们看一道关于伽玛分布的概率题：有A和B两路公交车，平均发车间隔时间分别是5分钟和10分钟，某乘客在站点S可以任意选择两者之一乘坐，假设A和B到达S的时刻无法确定，计算该乘客的平均等待公交车的时间。

可以将“假设A和B到达S的时刻无法确定”理解为公交车到达S站点的时刻是完全随机的，因此单位时间之内到达S站点的公交车次数符合泊松分布，而前后两辆公交车的时间差符合k=1的伽玛分布。下面我们先用随机数模拟的方法求出近似解，然后推导出解的公式。

T = 100000
A_count = T / 5
B_count = T / 10

A_time = np.random.uniform(0, T, A_count) ?
B_time = np.random.uniform(0, T, B_count)

bus_time = np.concatenate((A_time, B_time)) ?
bus_time.sort( )

N = 200000
passenger_time = np.random.uniform(bus_time[0], bus_time[-1], N) ?

idx = np.searchsorted(bus_time, passenger_time) ?
np.mean(bus_time[idx] - passenger_time) * 60    ?
199.12512768644049

模拟的总时间为T分钟，在这段之间之内，应该有A_count次A路公交车和B_count次B路公交车到达S站点。?可以用均匀分布uniform( )产生两路公交车到达S站点的时刻，?将这两个保存时刻的数组连接起来，并进行排序。

?在第一趟和最后一趟公交车的到达时间之间，产生乘客随机到达S站点的时刻数组。?在已经排序的公交车到达时刻数组bus_time中使用二分法搜索每个乘客到达时刻所在的下标数组idx。?bus_time[idx]就是乘客到达车站之后第一个到达车站的公交车的时刻，因此只需要计算其差值，并求平均值即可。通过随机数模拟得出的平均等待时间约为200秒。

将A和B两路汽车一起考虑，前后两个车次的平均间隔也为200秒。这似乎有些不可思议，直觉上我们可能期待一个小于平均间隔的等待时间。

np.mean(np.diff(bus_time)) * 60
199.98208112933918

这是因为存在观察者偏差，即会有更多的乘客出现在时间间隔较长的时间段。我们可以想象如果公交车因为事故晚点很长时间，那么通常车站上会挤满等待的人。在图7（上）中，蓝色竖线代表公交车的到站时刻，红色竖线代表乘客的到站时刻。可以看出，两条蓝色竖线之间的距离越大，其间的红色竖线就会越多。图7（下）的横轴是前后两辆公交车的时间差，纵轴是这段时间差之内的等待人数，可以看出二者成正比关系。

图7　观察者偏差

通过以上分析，不难写出计算平均等待时间的计算公式：

在公式中，x是两辆公交车之间的间隔时间，f(x)dx是时间间隔为x出现的概率。由于观察者效应，乘客出现在较长时间间隔的概率也较大，因此xf(x)dx可以看作与乘客出现在时间间隔为x时段的概率成比例的量，分母的积分将其归一化。而分子中的x/2是在该时间间隔段到达车站所需的平均等待时间。下面我们计算该公式，由图3-17可知，公交车的间隔几乎不会超过30分钟，因此虽然公式中的积分上限为+∞，但在实际计算时只需要指定一个较大的数即可。在本章后续的小节中会详细介绍数值积分quad( )的用法。

from scipy import integrate
t = 10.0 / 3  # 两辆公交车之间的平均时间间隔
bus_interval = stats.gamma(1, scale=t)
n, _ = integrate.quad(lambda x: 0.5 * x * x * bus_interval.pdf(x), 0, 1000)
d, _ = integrate.quad(lambda x: x * bus_interval.pdf(x), 0, 1000)
n / d * 60
200.0

五、学生t-分布与t检验

从均值为μ的正态分布中，抽取有n个值的样本，计算样本均值和样本方差s：

则符合df=n-1的学生t-分布。t值是抽选的样本的平均值与整体样本的期望值之差经过正规化之后的数值，可以用来描述抽取的样本与整体样本之间的差异。

下面的程序模拟学生t-分布（参见图8），?创建一个形状为(100000, 10)的正态分布的随机数数组，?使用上面的公式计算t值，?统计t值的分布情况并与stats.t的概率密度函数进行比较。如果我们使用10个样本计算t值，则它应该符合df=9的学生t-分布。

mu = 0.0
n = 10
samples = stats.norm(mu).rvs(size=(100000, n)) ?
t_samples = (np.mean(samples, axis=1) - mu) / np.std(samples, ddof=1, axis=1) * n**0.5 ?
sample_dist, x = np.histogram(t_samples, bins=100, density=True) ?
x = 0.5 * (x[:-1] + x[1:])
t_dist = stats.t(n-1).pdf(x)
print "max error:", np.max(np.abs(sample_dist - t_dist))
max error: 0.00658734287935

图8　模拟学生t-分布

图9绘制df为5和39的概率密度函数和生存函数，当df增大时，学生t-分布趋向于正态分布。

图9　当df增大时，学生t-分布趋向于正态分布

学生t-分布可以用于检测样本的平均值，下面我们从一个期望值为1的正态分布的随机变量中取30个数值：

n = 30
np.random.seed(42)
s = stats.norm.rvs(loc=1, scale=0.8, size=n)

我们建立整体样本的期望值为0.5的零假设，并用stats.ttest_1samp( )检验零假设是否能够被推翻。stats.ttest_1samp( )返回的第一个值为使用前述公式计算的t值，第二个值被称为p值。当p值小于0.05时，通常我们认为零假设不成立。因此下面的测试表明我们可以拒绝整体样本的期望值为0.5的假设。

零假设

在统计学中，零假设或虚无假设(null hypothesis)是做统计检验时的一类假设。零假设的内容一般是希望被证明为错误的假设或是需要着重考虑的假设。

t = (np.mean(s) - 0.5) / (np.std(s, ddof=1) / np.sqrt(n))
print t, stats.ttest_1samp(s, 0.5)
2.65858434088 (2.6585843408822241, 0.012637702257091229)

下面我们检验期望值是否为1.0，由于p值大于0.05，我们不能推翻期望值为1.0的零假设，但这并不意味着可以接受该假设，因为期望值为0.9的假设对应的p值也大于0.05，甚至比1.0的t值还大。

print (np.mean(s) - 1) / (np.std(s, ddof=1) / np.sqrt(n))
print stats.ttest_1samp(s, 1), stats.ttest_1samp(s, 0.9)
-1.14501736704
(-1.1450173670383303, 0.26156414618801477) (-0.38429702545421962, 0.70356191034252025)

通过ttest_1samp( )计算的p值就是图10中红色部分的面积。可以这样理解p值的含义：如果随机变量的期望值真的和假设的相同，那么从这个随机变量中随机抽取n个数值，其t值比测试样本的t值还极端（绝对值大）的可能性为p。因此当p很小时，我们可以推断假设不太可能成立。反过来当p值较大时，则不能推翻零假设，注意不能推翻假设并不代表能接受该假设。

?图10　红色部分为ttest_1samp( )计算的p值

拿上面期望值为0.5的测试为例：如果整体样本的期望值真的是0.5，那么抽取到t值大于2.65858或小于-2.65858的样本的概率为0.0126377。由于这个概率太小，因此整体样本的期望值应该不是0.5。也可以这样理解：如果整体样本的期望值为0.5，那么随机抽取比s更极端的样本的概率为0.0126377。

x = np.linspace(-5, 5, 500)
y = stats.t(n-1).pdf(x)
plt.plot(x, y, lw=2)
t, p = stats.ttest_1samp(s, 0.5)
mask = x > np.abs(t)
plt.fill_between(x[mask], y[mask], color="red", alpha=0.5)
mask = x < -np.abs(t)
plt.fill_between(x[mask], y[mask], color="red", alpha=0.5)
plt.axhline(color="k", lw=0.5)
plt.xlim(-5, 5)

下面用scipy.integrate.trapz( )积分验证p值，由于左右两块红色面积是相等的，下面的积分只需要计算其中一块的面积：

from scipy import integrate
x = np.linspace(-10, 10, 100000)
y = stats.t(n-1).pdf(x)
mask = x >= np.abs(t)
integrate.trapz(y[mask], x[mask])*2
0.012633433707685974

下面我们用随机数验证前面计算的p值，我们创建m组随机数，每组都有n个数值，然后计算假设总体样本期望值为0.5时每组随机数对应的t值tr，它是一个长度为m的数组。将tr与样本的t值ts进行绝对值大小比较，当|tr|>|ts|时，就说明该组随机数比样本更极端，统计极端组出现的概率，可以看到它和p值是相同的。

m = 200000
mean = 0.5
r = stats.norm.rvs(loc=mean, scale=0.8, size=(m, n))
ts = (np.mean(s) - mean) / (np.std(s, ddof=1) / np.sqrt(n))
tr = (np.mean(r, axis=1) - mean) / (np.std(r, ddof=1, axis=1) / np.sqrt(n))
np.mean(np.abs(tr) > np.abs(ts))
0.012695

如果s1和s2是两个独立的来自正态分布总体的样本，可以通过ttest_ind( )检验这两个总体的均值是否存在差异。通过equal_var参数指定两个总体的方差是否相同。在下面的例子中，?由于s1和s2样本来自不同方差的总体，因此equal_var参数为False。由于p<0.05，因此认为两个总体的均值存在差异。?s2和s3来自相同方差的总体，因此equal_var参数为True，所得到的p值很大，因此无法推翻零假设，也就无法否定两个总体的均值相同的零假设。

np.random.seed(42)

s1 = stats.norm.rvs(loc=1, scale=1.0, size=20)
s2 = stats.norm.rvs(loc=1.5, scale=0.5, size=20)
s3 = stats.norm.rvs(loc=1.5, scale=0.5, size=25)

print stats.ttest_ind(s1, s2, equal_var=False) ?
print stats.ttest_ind(s2, s3, equal_var=True)  ?
(-2.2391470627176755, 0.033250866086743665)
(-0.59466985218561719, 0.55518058758105393)

六、卡方分布和卡方检验

卡方分布（χ2）是概率论与统计学中常用的一种概率分布。k个独立的标准正态分布变量的平方和服从自由度为k的卡方分布。下面通过随机数验证该分布，结果如图11所示：

a = np.random.normal(size=(300000, 4))
cs = np.sum(a**2, axis=1)

sample_dist, bins = np.histogram(cs, bins=100, range=(0, 20), density=True)
x = 0.5 * (bins[:-1] + bins[1:])
chi2_dist = stats.chi2.pdf(x, 4)
print "max error:", np.max(np.abs(sample_dist - chi2_dist))
max error: 0.00340194486328

图11　使用随机数验证卡方分布

卡方分布可以用来描述这样的概率现象：袋子里有5种颜色的球，抽到每种球的概率相同，从中选N次，并统计每种颜色的次数Oi。则下面的χ2符合自由度为4的卡方分布，其中E=N/5为每种球被抽选的期望次数：

?

下面用程序模拟这个过程，结果如图12所示。?首先调用randint( )创建从0到5的随机数，其结果ball_ids的第0轴表示实验次数，第1轴为每次实验抽取的100个球的编号。?使用bincount( )统计每次实验中每个编号出现的次数，由于它不支持多维数组，因此这里使用apply_along_axis( )对第1轴上的数据循环调用bincount( )。为了保证每行的统计结果的长度相同，设置minlength参数为5，apply_along_axis( )会将所有关键字参数传递给进行实际运算的函数。?使用上面的公式计算χ2统计量cs2，?并用gaussian_kde( )计算cs2的分布情况。

repeat_count = 60000
n, k = 100, 5

np.random.seed(42)
ball_ids = np.random.randint(0, k, size=(repeat_count, n)) ?
counts = np.apply_along_axis(np.bincount, 1, ball_ids, minlength=k) ?
cs2 = np.sum((counts - n/k)**2.0/(n/k), axis=1) ?
k = stats.kde.gaussian_kde(cs2) ?
x = np.linspace(0, 10, 200)
pl.plot(x, stats.chi2.pdf(x, 4), lw=2, label=u"$\chi ^{2}$分布")
pl.plot(x, k(x), lw=2, color="red", alpha=0.6, label=u"样本分布")
pl.legend(loc="best")
pl.xlim(0, 10)

图12　模拟卡方分布

卡方检验可以用来评估观测值与理论值的差异是否只是因为随机误差造成的。在下面的例子中，袋子中各种颜色的球按照probabilities参数指定的概率分布，choose_balls(probabilities, size)从袋中选择size次并返回每种球被选中的次数。袋子1中的球的概率分布为：0.18、0.24、0.25、0.16、0.17。袋子2中各种颜色的球的个数一样多。通过调用choose_balls( )得到两组数字：80 93 97 64 66和89 76 79 71 85。现在需要判断袋子中的球是否是平均分布的。?

def choose_balls(probabilities, size):
r = stats.rv_discrete(values=(range(len(probabilities)), probabilities))
s = r.rvs(size=size)
counts = np.bincount(s)    
return counts

np.random.seed(42)
counts1 = choose_balls([0.18, 0.24, 0.25, 0.16, 0.17], 400)
counts2 = choose_balls([0.2]*5, 400)
counts1               counts2       
--------------------  --------------------
[80, 93, 97, 64, 66]  [89, 76, 79, 71, 85]

使用chisquare( )进行卡方检验，它的参数为每种球被选中次数的列表，如果没有设置检验的目标概率，就测试它们是否符合平均分布。卡方检验的零假设为样本符合目标概率，由下面的检验结果可知，第一个袋子对应的p值只有0.02，也就是说如果第一个袋子中的球真的符合平均分布，那么得到的观测结果80 93 97 64 66的概率只有2%，因此可以推翻零假设，即袋子中的球不太可能是平均分布的。第二个袋子对应的p值为0.64，无法推翻零假设，即我们的结论是不能否定第二个袋子中的球是平均分布的。注意，和前面介绍的t检验一样，零假设只能用来否定，因此不能根据观测结果89 76 79 71 85得出袋子中的球是符合平均分布的结论。?

chi1, p1 = stats.chisquare(counts1)
chi2, p2 = stats.chisquare(counts2)

print "chi1 =", chi1, "p1 =", p1
print "chi2 =", chi2, "p2 =", p2
chi1 = 11.375 p1 = 0.0226576012398
chi2 = 2.55 p2 = 0.635705452704

卡方检验是通过卡方分布进行计算的。图13显示自由度为4的卡方分布的概率密度函数，以及chi1和chi2对应的位置和。p1是右侧部分的面积，而p2是右侧的面积。

图13　卡方检验计算的概率为阴影部分的面积

卡方检验还可以用于二维数据。前面介绍的彩色球的例子中，只是按照球的颜色分组，而二维数据则按照样本的两个属性分组，统计学上称之为列联表。例如下面是性别与惯用手的列联表，我们希望知道这个统计结果能否说明性别与惯用手之间存在某种联系。

下面使用chi2_contingency( )对列联表进行卡方检验，零假设为性别与惯用手之间不存在联系，即男性与女性惯用左右手的概率相同。由于p值为0.3，因此不能推翻零假设，即该实验数据中没有明显证据表明男性和女性在使用左右手习惯上存在区别。

table = [[43, 9], [44, 4]]
chi2, p, dof, expected = stats.chi2_contingency(table)
chi2                  p         
------------------  -------------------
1.0724852071005921  0.30038477039056899

对于上面的2×2的数值较小的列联表，可以使用fisher_exact( )计算出精确的p值：?

stats.fisher_exact(table)
(0.43434343434343436, 0.23915695682225618)