44

How do I do an F-test to check if the variance is equivalent in two vectors in Python?

For example if I have

a = [1,2,1,2,1,2,1,2,1,2]
b = [1,3,-1,2,1,5,-1,6,-1,2]

is there something similar to

scipy.stats.ttest_ind(a, b)

I found

sp.stats.f(a, b)

But it appears to be something different to an F-test

1
  • 2
    If you just want to test for equal variance, scipy has Bartlett and Levene tests.
    – BrenBarn
    Feb 1, 2014 at 5:08

5 Answers 5

57

The test statistic F test for equal variances is simply:

F = Var(X) / Var(Y)

Where F is distributed as df1 = len(X) - 1, df2 = len(Y) - 1

scipy.stats.f which you mentioned in your question has a CDF method. This means you can generate a p-value for the given statistic and test whether that p-value is greater than your chosen alpha level.

Thus:

alpha = 0.05 #Or whatever you want your alpha to be.
p_value = scipy.stats.f.cdf(F, df1, df2)
if p_value > alpha:
    # Reject the null hypothesis that Var(X) == Var(Y)

Note that the F-test is extremely sensitive to non-normality of X and Y, so you're probably better off doing a more robust test such as Levene's test or Bartlett's test unless you're reasonably sure that X and Y are distributed normally. These tests can be found in the scipy api:

3
  • 40
    This answer is mathematically correct but conceptually wrong. The p_value is rather 1 - CDF. This is the area under the probability distribution above your F_critical. I would write it this way: p_value = scipy.stats.f.sf(F, df1, df2) I think the survival function calculates this slightly different so you want to double check in case you are handling very small p_values (due to multiple comparisons) Then you reject the null hypothesis if p_value < alpha
    – Diego
    Nov 6, 2014 at 21:38
  • 1
    Are the variances not significantly different if the p-value is less than 0.05 OR greater than 0.95? If F = 0.25, i.e. the variance of Y is 4 times greater than the variance of X, (with say df1 = df2 = 10), then p = 0.98 based on 1 - CDF. So, I would go for less than 0.05 OR greater than 0.95 to test for significantly different variances. Apr 24, 2018 at 2:19
  • 1
    I would personally advised against using Bartlett's test for non-normal data as the statistic is is not reliable with "moderate" departures from normality. Better off with Levene's. Aug 1, 2019 at 10:24
13

For anyone who came here searching for an ANOVA F-test or to compare between models for feature selection

10

To do a one way anova you can use

import scipy.stats as stats

stats.f_oneway(a,b)

One way Anova checks if the variance between the groups is greater then the variance within groups, and computes the probability of observing this variance ratio using F-distribution. A good tutorial can be found here:

https://www.khanacademy.org/math/probability/statistics-inferential/anova/v/anova-1-calculating-sst-total-sum-of-squares

1
  • 1
    I assume the OP asked about comparing the variances of two samples, for which the test statistic is the ratio of the sample variances. If you run ANOVA on two groups, the omnibus F test statistic value will be the ratio of the between-group variance vs within-group variance, which is not the same thing. Try it! :-) Aug 4, 2020 at 16:27
4

if you need a two-tailed test, you can proceed as follow, i choosed 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)

if you want to proceed to an ANOVA like test where only large values can cause rejection, you can proceed to right-tail test, you need to pay attention to the order of variances (fstatistics = var1/var2 or 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)

The left-tailed can be done as follow :

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)
4

Here is a simple function to calculate the one-sided or two-sided F test with Python and SciPy. The results have been checked against the output of the var.test() function in R. Please keep in mind the warnings mentioned in the other answers concerning the sensitivity of the F-test to non-normality.

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
2
  • In standard statistical practice, ddof=1 provides an unbiased estimator of the variance of a hypothetical infinite population. ddof=0 provides a maximum likelihood estimate of the variance for normally distributed variables.
    – JustDoIt
    Feb 17, 2022 at 9:36
  • 1
    @JustDoIt Thanks but my goal with this Python function was to reproduce exactly the output of the var.test() function in R. Hence I'd rather not mess with the ddof settings. Feb 17, 2022 at 11:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.