如何进行 F 检验来检查 Python 中两个向量的方差是否相等?
例如,如果我有
a = [1,2,1,2,1,2,1,2,1,2]
b = [1,3,-1,2,1,5,-1,6,-1,2]
有类似的东西吗
scipy.stats.ttest_ind(a, b)
我发现了
sp.stats.f(a, b)
但这似乎与 F 测试不同
等方差的检验统计量 F 检验很简单:
F = Var(X) / Var(Y)
其中
F
分布为 df1 = len(X) - 1, df2 = len(Y) - 1
您在问题中提到的scipy.stats.f
有一个CDF方法。这意味着您可以为给定的统计量生成 p 值,并测试该 p 值是否大于您选择的 alpha 水平。
因此:
alpha = 0.05 #Or whatever you want your alpha to be.
p_value = 1 - scipy.stats.f.cdf(F, df1, df2)
或者,可以使用 f 检验的 sf 方法计算 p 值:
p_value = scipy.stats.f.sf(F, df1, df2)
if p_value > alpha:
# Reject the null hypothesis that Var(X) == Var(Y)
请注意,F 检验对 X 和 Y 的非正态性极其敏感,因此您最好进行更稳健的检验,例如 Levene 检验 或 Bartlett 检验,除非您相当确定 X和 Y 呈正态分布。这些测试可以在
scipy
api 中找到:
对于来这里寻找 ANOVA F 检验或比较模型以进行特征选择的任何人
sklearn.feature_selection.f_classif
进行方差分析测试,并且sklearn.feature_selection.f_regression
对回归进行顺序测试要进行单向方差分析,您可以使用
import scipy.stats as stats
stats.f_oneway(a,b)
方差分析的一种方法是检查组之间的方差是否大于组内的方差,并使用 F 分布计算观察到该方差比的概率。可以在这里找到一个很好的教程:
这是一个使用 Python 和 SciPy 计算单边或双边 F 检验的简单函数。结果已根据 R 中
var.test()
函数的输出进行了检查。请记住其他答案中提到的有关 F 检验对非正态性的敏感性的警告。
import scipy.stats as st
def f_test(x, y, alt="two_sided"):
"""
Calculates the F-test.
:param x: The first group of data
:param y: The second group of data
:param alt: The alternative hypothesis, one of "two_sided" (default), "greater" or "less"
:return: a tuple with the F statistic value and the p-value.
"""
df1 = len(x) - 1
df2 = len(y) - 1
f = x.var() / y.var()
if alt == "greater":
p = 1.0 - st.f.cdf(f, df1, df2)
elif alt == "less":
p = st.f.cdf(f, df1, df2)
else:
# two-sided by default
# Crawley, the R book, p.355
p = 2.0*(1.0 - st.f.cdf(f, df1, df2))
return f, p
如果需要双尾测试,可以如下进行,我选择alpha=0.05:
a = [1,2,1,2,1,2,1,2,1,2]
b = [1,3,-1,2,1,5,-1,6,-1,2]
print('Variance a={0:.3f}, Variance b={1:.3f}'.format(np.var(a, ddof=1), np.var(b, ddof=1)))
fstatistics = np.var(a, ddof=1)/np.var(b, ddof=1) # because we estimate mean from data
fdistribution = stats.f(len(a)-1,len(b)-1) # build an F-distribution object
p_value = 2*min(fdistribution.cdf(f_critical), 1-fdistribution.cdf(f_critical))
f_critical1 = fdistribution.ppf(0.025)
f_critical2 = fdistribution.ppf(0.975)
print(fstatistics,f_critical1, f_critical2 )
if (p_value<0.05):
print('Reject H0', p_value)
else:
print('Cant Reject H0', p_value)
如果您想继续进行类似方差分析的测试,其中只有大值才会导致拒绝,您可以继续进行右尾测试,您需要注意方差的顺序(fstatistics = var1/var2 或 var2/var1):
a = [1,2,1,2,1,2,1,2,1,2]
b = [1,3,-1,2,1,5,-1,6,-1,2]
print('Variance a={0:.3f}, Variance b={1:.3f}'.format(np.var(a, ddof=1), np.var(b, ddof=1)))
fstatistics = max(np.var(a, ddof=1), np.var(b, ddof=1))/min(np.var(a, ddof=1), np.var(b, ddof=1)) # because we estimate mean from data
fdistribution = stats.f(len(a)-1,len(b)-1) # build an F-distribution object
p_value = 1-fdistribution.cdf(fstatistics)
f_critical = fd.ppf(0.95)
print(fstatistics, f_critical)
if (p_value<0.05):
print('Reject H0', p_value)
else:
print('Cant Reject H0', p_value)
左尾可以如下完成:
a = [1,2,1,2,1,2,1,2,1,2]
b = [1,3,-1,2,1,5,-1,6,-1,2]
print('Variance a={0:.3f}, Variance b={1:.3f}'.format(np.var(a, ddof=1), np.var(b, ddof=1)))
fstatistics = min(np.var(a, ddof=1), np.var(b, ddof=1))/max(np.var(a, ddof=1), np.var(b, ddof=1)) # because we estimate mean from data
fdistribution = stats.f(len(a)-1,len(b)-1) # build an F-distribution object
p_value = fdistribution.cdf(fstatistics)
f_critical = fd.ppf(0.05)
print(fstatistics, f_critical)
if (p_value<0.05):
print('Reject H0', p_value)
else:
print('Cant Reject H0', p_value)