Python科学计算(第2版)
上QQ阅读APP看书,第一时间看更新

3.4.5 学生t-分布与t检验

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

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

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

图3-18 模拟学生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

图3-19绘制df为5和39的概率密度函数和生存函数,当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值就是图3-20中红色部分的面积。可以这样理解p值的含义:如果随机变量的期望值真的和假设的相同,那么从这个随机变量中随机抽取n个数值,其t值比测试样本的t值还极端(绝对值大)的可能性为p。因此当p很小时,我们可以推断假设不太可能成立。反过来当p值较大时,则不能推翻零假设,注意不能推翻假设并不代表能接受该假设。

图3-20 红色部分为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)