Python: ttest

Description

NCL’s ttest function reprogrammed in Python. It returns an estimate of the statistical significance.

This function is the Python equivalent of the NCL function “ttest”, yyr/ncl

NCL documentation ttest: https://www.ncl.ucar.edu/Document/Functions/Built-in/ttest.shtml

Software requirements

• Python 3

• numpy

• scipy

Example script

ttest.py

def ttest(ave1, var1, s1, ave2, var2, s2, iflag=True):
    '''
    NCL's ttest function reprogrammed in Python. 
    It returns an estimate of the statistical significance.

    This function is the Python equivalent of the NCL function "ttest", 
        https://github.com/yyr/ncl/blob/master/ni/src/lib/nfpfort/ttest.f
    NCL documentation ttest: 
        https://www.ncl.ucar.edu/Document/Functions/Built-in/ttest.shtml

	The returned probability is two-tailed. The ttest uses the incomplete beta 
	function (betainc) to calculate the probability. Example 2 at betainc 
	illustrates how to get the one-tailed probability.

    Parameter:
        ave1, ave2  average of sample 1, average of sample 2
        var1, var2  variance of sample 1, variance  of sample 2
        s1, s2      size  of sample 1, size  of sample 2
        iflag       False: The two samples are assumed to have same population variance
                    True:  The two samples are assumed to have different population 
                           variances                    
        Ideally, s1 and s2 >=30

    Returns: 
        Level of significance probability that a random variable following the 
        students t distribution will exceed tval in absolute value.

    Fortran subroutine name: DTTEST(AVE1,VAR1,S1,AVE2,VAR2,S2,IFLAG,ALPHA,TVAL,IER)
    NCL function name:       ttest (ave1,var1,s1,ave2,var2,s2,iflag,tval_opt)
    '''
    import numpy as np
    from scipy import stats, special
    
    ier = 0
    if s1 < 2. or s2 < 2.: ier = 1
    if var1.any() <= 0 or var2.any() <= 0 : ier = ier + 10
    if ier != 0: 
        print(f'ERROR: ttest returns error {ier}')

    if iflag:
        df = s1 + s2 - 2.
        varp = ((s1-1.)*var1 + (s2-1.)*var2)/df
        tval = (ave1-ave2)/np.sqrt(varp* (1./s1+1./s2))
    else:
        df   = (var1/s1 + var2/s2)**2 / ((var1/s1)**2 / (s1-1.)+(var2/s2)**2 / (s2-1.))
        tval = (ave1-ave2)/np.sqrt(var1/s1+var2/s2)
    
    x = df/(df+tval**2)  # x:    Upper limit of integration. X must be in (0,1) inclusive.
    pin = 0.5*df         # pin:  First beta distribution parameter.  PIN must be > 0.0.
    qin = 0.5            # qin:  Second beta distribution parameter. QIN must be > 0.0.

    return special.betainc(pin, qin, x)