进行双样本t检验

Question

我有两个样本的平均值、标准差和样本数量n，这两个样本是从同一个总体中抽取的，但是由不同的实验室测量的。

样本1和样本2的样本数量n是不一样的。我想进行一个加权的双尾t检验，也就是要考虑样本数量n。

我试着使用了scipy.stat模块，通过np.random.normal生成我的数据，因为这个模块只接受数据，而不接受像平均值和标准差这样的统计值（有没有办法直接使用这些值呢）。但是这个方法不行，因为数据数组的大小必须相等。

如果能帮我找到计算p值的方法，我将非常感激。

Answer 1

在最近版本的Scipy 0.12.0中，这个功能已经内置了（实际上可以处理不同大小的样本）。在scipy.stats模块中，ttest_ind这个函数在equal_var参数设置为False时，会执行Welch的t检验。

举个例子：

>>> import scipy.stats as stats
>>> sample1 = np.random.randn(10, 1)
>>> sample2 = 1 + np.random.randn(15, 1)
>>> t_stat, p_val = stats.ttest_ind(sample1, sample2, equal_var=False)
>>> t_stat
array([-3.94339083])
>>> p_val
array([ 0.00070813])

Answer 2

如果你有原始数据，分别存放在数组 a 和 b 中，你可以使用 scipy.stats.ttest_ind 这个工具，并把参数 equal_var=False 加上。

t, p = ttest_ind(a, b, equal_var=False)

如果你只有这两个数据集的总结统计信息，你可以用 scipy.stats.ttest_ind_from_stats 来计算 t 值（这个功能在 scipy 0.16 版本中新增）或者使用公式来计算（可以参考这个链接：http://en.wikipedia.org/wiki/Welch%27s_t_test）。

下面的脚本展示了这些可能性。

from __future__ import print_function

import numpy as np
from scipy.stats import ttest_ind, ttest_ind_from_stats
from scipy.special import stdtr

np.random.seed(1)

# Create sample data.
a = np.random.randn(40)
b = 4*np.random.randn(50)

# Use scipy.stats.ttest_ind.
t, p = ttest_ind(a, b, equal_var=False)
print("ttest_ind:            t = %g  p = %g" % (t, p))

# Compute the descriptive statistics of a and b.
abar = a.mean()
avar = a.var(ddof=1)
na = a.size
adof = na - 1

bbar = b.mean()
bvar = b.var(ddof=1)
nb = b.size
bdof = nb - 1

# Use scipy.stats.ttest_ind_from_stats.
t2, p2 = ttest_ind_from_stats(abar, np.sqrt(avar), na,
                              bbar, np.sqrt(bvar), nb,
                              equal_var=False)
print("ttest_ind_from_stats: t = %g  p = %g" % (t2, p2))

# Use the formulas directly.
tf = (abar - bbar) / np.sqrt(avar/na + bvar/nb)
dof = (avar/na + bvar/nb)**2 / (avar**2/(na**2*adof) + bvar**2/(nb**2*bdof))
pf = 2*stdtr(dof, -np.abs(tf))

print("formula:              t = %g  p = %g" % (tf, pf))

输出结果：

ttest_ind:            t = -1.5827  p = 0.118873
ttest_ind_from_stats: t = -1.5827  p = 0.118873
formula:              t = -1.5827  p = 0.118873