<s> '''<STR_LIT>''' <EOL> from __future__ import print_function <EOL> from statsmodels . compat . python import lzip , range , lrange , zip <EOL> import scipy . stats <EOL> import numpy <EOL> import numpy as np <EOL> import math <EOL> import copy <EOL> from scipy import stats <EOL> from statsmodels . iolib . table import SimpleTable <EOL> from numpy . testing import assert_almost_equal , assert_equal <EOL> from statsmodels . stats . multitest import multipletests , _ecdf as ecdf , fdrcorrection as fdrcorrection0 , fdrcorrection_twostage <EOL> from statsmodels . graphics import utils <EOL> qcrit = '''<STR_LIT>''' <EOL> res = [ line . split ( ) for line in qcrit . replace ( '<STR_LIT>' , '<STR_LIT>' ) . split ( '<STR_LIT:\n>' ) ] <EOL> c = np . array ( res [ <NUM_LIT:2> : - <NUM_LIT:1> ] ) . astype ( float ) <EOL> ccols = np . arange ( <NUM_LIT:2> , <NUM_LIT:11> ) <EOL> crows = c [ : , <NUM_LIT:0> ] <EOL> cv005 = c [ : , <NUM_LIT:1> : : <NUM_LIT:2> ] <EOL> cv001 = c [ : , <NUM_LIT:2> : : <NUM_LIT:2> ] <EOL> from scipy import interpolate <EOL> def get_tukeyQcrit ( k , df , alpha = <NUM_LIT> ) : <EOL> '''<STR_LIT>''' <EOL> if alpha == <NUM_LIT> : <EOL> intp = interpolate . interp1d ( crows , cv005 [ : , k - <NUM_LIT:2> ] ) <EOL> elif alpha == <NUM_LIT> : <EOL> intp = interpolate . interp1d ( crows , cv001 [ : , k - <NUM_LIT:2> ] ) <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> return intp ( df ) <EOL> def get_tukeyQcrit2 ( k , df , alpha = <NUM_LIT> ) : <EOL> '''<STR_LIT>''' <EOL> from statsmodels . stats . libqsturng import qsturng <EOL> return qsturng ( <NUM_LIT:1> - alpha , k , df ) <EOL> def Tukeythreegene ( first , second , third ) : <EOL> firstmean = numpy . mean ( first ) <EOL> secondmean = numpy . mean ( second ) <EOL> thirdmean = numpy . mean ( third ) <EOL> firststd = numpy . std ( first ) <EOL> secondstd = numpy . std ( second ) <EOL> thirdstd = numpy . std ( third ) <EOL> firsts2 = math . pow ( firststd , <NUM_LIT:2> ) <EOL> seconds2 = math . pow ( secondstd , <NUM_LIT:2> ) <EOL> thirds2 = math . pow ( thirdstd , <NUM_LIT:2> ) <EOL> mserrornum = firsts2 * <NUM_LIT:2> + seconds2 * <NUM_LIT:2> + thirds2 * <NUM_LIT:2> <EOL> mserrorden = ( len ( first ) + len ( second ) + len ( third ) ) - <NUM_LIT:3> <EOL> mserror = mserrornum / mserrorden <EOL> standarderror = math . sqrt ( mserror / len ( first ) ) <EOL> dftotal = len ( first ) + len ( second ) + len ( third ) - <NUM_LIT:1> <EOL> dfgroups = <NUM_LIT:2> <EOL> dferror = dftotal - dfgroups <EOL> qcrit = <NUM_LIT:0.5> <EOL> qcrit = get_tukeyQcrit ( <NUM_LIT:3> , dftotal , alpha = <NUM_LIT> ) <EOL> qtest3to1 = ( math . fabs ( thirdmean - firstmean ) ) / standarderror <EOL> qtest3to2 = ( math . fabs ( thirdmean - secondmean ) ) / standarderror <EOL> qtest2to1 = ( math . fabs ( secondmean - firstmean ) ) / standarderror <EOL> conclusion = [ ] <EOL> print ( qtest3to1 ) <EOL> print ( qtest3to2 ) <EOL> print ( qtest2to1 ) <EOL> if ( qtest3to1 > qcrit ) : <EOL> conclusion . append ( '<STR_LIT>' ) <EOL> else : <EOL> conclusion . append ( '<STR_LIT>' ) <EOL> if ( qtest3to2 > qcrit ) : <EOL> conclusion . append ( '<STR_LIT>' ) <EOL> else : <EOL> conclusion . append ( '<STR_LIT>' ) <EOL> if ( qtest2to1 > qcrit ) : <EOL> conclusion . append ( '<STR_LIT>' ) <EOL> else : <EOL> conclusion . append ( '<STR_LIT>' ) <EOL> return conclusion <EOL> def Tukeythreegene2 ( genes ) : <EOL> """<STR_LIT>""" <EOL> means = [ ] <EOL> stds = [ ] <EOL> for gene in genes : <EOL> means . append ( numpy . mean ( gene ) ) <EOL> std . append ( numpy . std ( gene ) ) <EOL> stds2 = [ ] <EOL> for std in stds : <EOL> stds2 . append ( math . pow ( std , <NUM_LIT:2> ) ) <EOL> mserrornum = sum ( stds2 ) * <NUM_LIT:2> <EOL> mserrorden = ( len ( genes [ <NUM_LIT:0> ] ) + len ( genes [ <NUM_LIT:1> ] ) + len ( genes [ <NUM_LIT:2> ] ) ) - <NUM_LIT:3> <EOL> mserror = mserrornum / mserrorden <EOL> def catstack ( args ) : <EOL> x = np . hstack ( args ) <EOL> labels = np . hstack ( [ k * np . ones ( len ( arr ) ) for k , arr in enumerate ( args ) ] ) <EOL> return x , labels <EOL> def maxzero ( x ) : <EOL> '''<STR_LIT>''' <EOL> x = np . asarray ( x ) <EOL> cond1 = x [ : - <NUM_LIT:1> ] < <NUM_LIT:0> <EOL> cond2 = x [ <NUM_LIT:1> : ] > <NUM_LIT:0> <EOL> allzeros = np . nonzero ( ( cond1 & cond2 ) | ( x [ <NUM_LIT:1> : ] == <NUM_LIT:0> ) ) [ <NUM_LIT:0> ] + <NUM_LIT:1> <EOL> if x [ - <NUM_LIT:1> ] >= <NUM_LIT:0> : <EOL> maxz = max ( allzeros ) <EOL> else : <EOL> maxz = None <EOL> return maxz , allzeros <EOL> def maxzerodown ( x ) : <EOL> '''<STR_LIT>''' <EOL> x = np . asarray ( x ) <EOL> cond1 = x [ : - <NUM_LIT:1> ] > <NUM_LIT:0> <EOL> cond2 = x [ <NUM_LIT:1> : ] < <NUM_LIT:0> <EOL> allzeros = np . nonzero ( ( cond1 & cond2 ) | ( x [ <NUM_LIT:1> : ] == <NUM_LIT:0> ) ) [ <NUM_LIT:0> ] + <NUM_LIT:1> <EOL> if x [ - <NUM_LIT:1> ] <= <NUM_LIT:0> : <EOL> maxz = max ( allzeros ) <EOL> else : <EOL> maxz = None <EOL> return maxz , allzeros <EOL> def rejectionline ( n , alpha = <NUM_LIT:0.5> ) : <EOL> '''<STR_LIT>''' <EOL> t = np . arange ( n ) / float ( n ) <EOL> frej = t / ( t * ( <NUM_LIT:1> - alpha ) + alpha ) <EOL> return frej <EOL> def fdrcorrection_bak ( pvals , alpha = <NUM_LIT> , method = '<STR_LIT>' ) : <EOL> '''<STR_LIT>''' <EOL> pvals = np . asarray ( pvals ) <EOL> pvals_sortind = np . argsort ( pvals ) <EOL> pvals_sorted = pvals [ pvals_sortind ] <EOL> pecdf = ecdf ( pvals_sorted ) <EOL> if method in [ '<STR_LIT:i>' , '<STR_LIT>' , '<STR_LIT:p>' , '<STR_LIT>' ] : <EOL> rline = pvals_sorted / alpha <EOL> elif method in [ '<STR_LIT:n>' , '<STR_LIT>' ] : <EOL> cm = np . sum ( <NUM_LIT:1.> / np . arange ( <NUM_LIT:1> , len ( pvals ) ) ) <EOL> rline = pvals_sorted / alpha * cm <EOL> elif method in [ '<STR_LIT:g>' , '<STR_LIT>' ] : <EOL> rline = pvals_sorted / ( pvals_sorted * ( <NUM_LIT:1> - alpha ) + alpha ) <EOL> elif method in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> cm = np . sum ( np . arange ( len ( pvals ) ) ) <EOL> rline = pvals_sorted / alpha / cm <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> reject = pecdf >= rline <EOL> if reject . any ( ) : <EOL> rejectmax = max ( np . nonzero ( reject ) [ <NUM_LIT:0> ] ) <EOL> else : <EOL> rejectmax = <NUM_LIT:0> <EOL> reject [ : rejectmax ] = True <EOL> return reject [ pvals_sortind . argsort ( ) ] <EOL> def mcfdr ( nrepl = <NUM_LIT:100> , nobs = <NUM_LIT:50> , ntests = <NUM_LIT:10> , ntrue = <NUM_LIT:6> , mu = <NUM_LIT:0.5> , alpha = <NUM_LIT> , rho = <NUM_LIT:0.> ) : <EOL> '''<STR_LIT>''' <EOL> nfalse = ntests - ntrue <EOL> locs = np . array ( [ <NUM_LIT:0.> ] * ntrue + [ mu ] * ( ntests - ntrue ) ) <EOL> results = [ ] <EOL> for i in range ( nrepl ) : <EOL> rvs = locs + randmvn ( rho , size = ( nobs , ntests ) ) <EOL> tt , tpval = stats . ttest_1samp ( rvs , <NUM_LIT:0> ) <EOL> res = fdrcorrection_bak ( np . abs ( tpval ) , alpha = alpha , method = '<STR_LIT:i>' ) <EOL> res0 = fdrcorrection0 ( np . abs ( tpval ) , alpha = alpha ) <EOL> results . append ( [ np . sum ( res [ : ntrue ] ) , np . sum ( res [ ntrue : ] ) ] + <EOL> [ np . sum ( res0 [ : ntrue ] ) , np . sum ( res0 [ ntrue : ] ) ] + <EOL> res . tolist ( ) + <EOL> np . sort ( tpval ) . tolist ( ) + <EOL> [ np . sum ( tpval [ : ntrue ] < alpha ) , <EOL> np . sum ( tpval [ ntrue : ] < alpha ) ] + <EOL> [ np . sum ( tpval [ : ntrue ] < alpha / ntests ) , <EOL> np . sum ( tpval [ ntrue : ] < alpha / ntests ) ] ) <EOL> return np . array ( results ) <EOL> def randmvn ( rho , size = ( <NUM_LIT:1> , <NUM_LIT:2> ) , standardize = False ) : <EOL> '''<STR_LIT>''' <EOL> nobs , nvars = size <EOL> if <NUM_LIT:0> < rho and rho < <NUM_LIT:1> : <EOL> rvs = np . random . randn ( nobs , nvars + <NUM_LIT:1> ) <EOL> rvs2 = rvs [ : , : - <NUM_LIT:1> ] * np . sqrt ( ( <NUM_LIT:1> - rho ) ) + rvs [ : , - <NUM_LIT:1> : ] * np . sqrt ( rho ) <EOL> elif rho == <NUM_LIT:0> : <EOL> rvs2 = np . random . randn ( nobs , nvars ) <EOL> elif rho < <NUM_LIT:0> : <EOL> if rho < - <NUM_LIT:1.> / ( nvars - <NUM_LIT:1> ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> elif rho == - <NUM_LIT:1.> / ( nvars - <NUM_LIT:1> ) : <EOL> rho = - <NUM_LIT:1.> / ( nvars - <NUM_LIT:1> + <NUM_LIT> ) <EOL> A = rho * np . ones ( ( nvars , nvars ) ) + ( <NUM_LIT:1> - rho ) * np . eye ( nvars ) <EOL> rvs2 = np . dot ( np . random . randn ( nobs , nvars ) , np . linalg . cholesky ( A ) . T ) <EOL> if standardize : <EOL> rvs2 = stats . zscore ( rvs2 ) <EOL> return rvs2 <EOL> def tiecorrect ( xranks ) : <EOL> '''<STR_LIT>''' <EOL> rankbincount = np . bincount ( np . asarray ( xranks , dtype = int ) ) <EOL> nties = rankbincount [ rankbincount > <NUM_LIT:1> ] <EOL> ntot = float ( len ( xranks ) ) ; <EOL> tiecorrection = <NUM_LIT:1> - ( nties ** <NUM_LIT:3> - nties ) . sum ( ) / ( ntot ** <NUM_LIT:3> - ntot ) <EOL> return tiecorrection <EOL> class GroupsStats ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __init__ ( self , x , useranks = False , uni = None , intlab = None ) : <EOL> '''<STR_LIT>''' <EOL> self . x = np . asarray ( x ) <EOL> if intlab is None : <EOL> uni , intlab = np . unique ( x [ : , <NUM_LIT:1> ] , return_inverse = True ) <EOL> elif uni is None : <EOL> uni = np . unique ( x [ : , <NUM_LIT:1> ] ) <EOL> self . useranks = useranks <EOL> self . uni = uni <EOL> self . intlab = intlab <EOL> self . groupnobs = groupnobs = np . bincount ( intlab ) <EOL> self . runbasic ( useranks = useranks ) <EOL> def runbasic_old ( self , useranks = False ) : <EOL> x = self . x <EOL> if useranks : <EOL> self . xx = x [ : , <NUM_LIT:1> ] . argsort ( ) . argsort ( ) + <NUM_LIT:1> <EOL> else : <EOL> self . xx = x [ : , <NUM_LIT:0> ] <EOL> self . groupsum = groupranksum = np . bincount ( self . intlab , weights = self . xx ) <EOL> self . groupmean = grouprankmean = groupranksum * <NUM_LIT:1.0> / self . groupnobs <EOL> self . groupmeanfilter = grouprankmean [ self . intlab ] <EOL> def runbasic ( self , useranks = False ) : <EOL> x = self . x <EOL> if useranks : <EOL> xuni , xintlab = np . unique ( x [ : , <NUM_LIT:0> ] , return_inverse = True ) <EOL> ranksraw = x [ : , <NUM_LIT:0> ] . argsort ( ) . argsort ( ) + <NUM_LIT:1> <EOL> self . xx = GroupsStats ( np . column_stack ( [ ranksraw , xintlab ] ) , <EOL> useranks = False ) . groupmeanfilter <EOL> else : <EOL> self . xx = x [ : , <NUM_LIT:0> ] <EOL> self . groupsum = groupranksum = np . bincount ( self . intlab , weights = self . xx ) <EOL> self . groupmean = grouprankmean = groupranksum * <NUM_LIT:1.0> / self . groupnobs <EOL> self . groupmeanfilter = grouprankmean [ self . intlab ] <EOL> def groupdemean ( self ) : <EOL> return self . xx - self . groupmeanfilter <EOL> def groupsswithin ( self ) : <EOL> xtmp = self . groupdemean ( ) <EOL> return np . bincount ( self . intlab , weights = xtmp ** <NUM_LIT:2> ) <EOL> def groupvarwithin ( self ) : <EOL> return self . groupsswithin ( ) / ( self . groupnobs - <NUM_LIT:1> ) <EOL> class TukeyHSDResults ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , mc_object , results_table , q_crit , reject = None , <EOL> meandiffs = None , std_pairs = None , confint = None , df_total = None , <EOL> reject2 = None , variance = None ) : <EOL> self . _multicomp = mc_object <EOL> self . _results_table = results_table <EOL> self . q_crit = q_crit <EOL> self . reject = reject <EOL> self . meandiffs = meandiffs <EOL> self . std_pairs = std_pairs <EOL> self . confint = confint <EOL> self . df_total = df_total <EOL> self . reject2 = reject2 <EOL> self . variance = variance <EOL> self . data = self . _multicomp . data <EOL> self . groups = self . _multicomp . groups <EOL> self . groupsunique = self . _multicomp . groupsunique <EOL> def __str__ ( self ) : <EOL> return str ( self . _results_table ) <EOL> def summary ( self ) : <EOL> '''<STR_LIT>''' <EOL> return self . _results_table <EOL> def _simultaneous_ci ( self ) : <EOL> """<STR_LIT>""" <EOL> self . halfwidths = simultaneous_ci ( self . q_crit , self . variance , <EOL> self . _multicomp . groupstats . groupnobs , <EOL> self . _multicomp . pairindices ) <EOL> def plot_simultaneous ( self , comparison_name = None , ax = None , figsize = ( <NUM_LIT:10> , <NUM_LIT:6> ) , <EOL> xlabel = None , ylabel = None ) : <EOL> """<STR_LIT>""" <EOL> fig , ax1 = utils . create_mpl_ax ( ax ) <EOL> if figsize is not None : <EOL> fig . set_size_inches ( figsize ) <EOL> if getattr ( self , '<STR_LIT>' , None ) is None : <EOL> self . _simultaneous_ci ( ) <EOL> means = self . _multicomp . groupstats . groupmean <EOL> sigidx = [ ] <EOL> nsigidx = [ ] <EOL> minrange = [ means [ i ] - self . halfwidths [ i ] for i in range ( len ( means ) ) ] <EOL> maxrange = [ means [ i ] + self . halfwidths [ i ] for i in range ( len ( means ) ) ] <EOL> if comparison_name is None : <EOL> ax1 . errorbar ( means , lrange ( len ( means ) ) , xerr = self . halfwidths , <EOL> marker = '<STR_LIT:o>' , linestyle = '<STR_LIT:None>' , color = '<STR_LIT:k>' , ecolor = '<STR_LIT:k>' ) <EOL> else : <EOL> if comparison_name not in self . groupsunique : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> midx = np . where ( self . groupsunique == comparison_name ) [ <NUM_LIT:0> ] <EOL> for i in range ( len ( means ) ) : <EOL> if self . groupsunique [ i ] == comparison_name : <EOL> continue <EOL> if ( min ( maxrange [ i ] , maxrange [ midx ] ) - <EOL> max ( minrange [ i ] , minrange [ midx ] ) < <NUM_LIT:0> ) : <EOL> sigidx . append ( i ) <EOL> else : <EOL> nsigidx . append ( i ) <EOL> ax1 . errorbar ( means [ midx ] , midx , xerr = self . halfwidths [ midx ] , <EOL> marker = '<STR_LIT:o>' , linestyle = '<STR_LIT:None>' , color = '<STR_LIT:b>' , ecolor = '<STR_LIT:b>' ) <EOL> ax1 . plot ( [ minrange [ midx ] ] * <NUM_LIT:2> , [ - <NUM_LIT:1> , self . _multicomp . ngroups ] , <EOL> linestyle = '<STR_LIT>' , color = '<STR_LIT>' ) <EOL> ax1 . plot ( [ maxrange [ midx ] ] * <NUM_LIT:2> , [ - <NUM_LIT:1> , self . _multicomp . ngroups ] , <EOL> linestyle = '<STR_LIT>' , color = '<STR_LIT>' ) <EOL> if len ( sigidx ) > <NUM_LIT:0> : <EOL> ax1 . errorbar ( means [ sigidx ] , sigidx , <EOL> xerr = self . halfwidths [ sigidx ] , marker = '<STR_LIT:o>' , <EOL> linestyle = '<STR_LIT:None>' , color = '<STR_LIT:r>' , ecolor = '<STR_LIT:r>' ) <EOL> if len ( nsigidx ) > <NUM_LIT:0> : <EOL> ax1 . errorbar ( means [ nsigidx ] , nsigidx , <EOL> xerr = self . halfwidths [ nsigidx ] , marker = '<STR_LIT:o>' , <EOL> linestyle = '<STR_LIT:None>' , color = '<STR_LIT>' , ecolor = '<STR_LIT>' ) <EOL> ax1 . set_title ( '<STR_LIT>' ) <EOL> r = np . max ( maxrange ) - np . min ( minrange ) <EOL> ax1 . set_ylim ( [ - <NUM_LIT:1> , self . _multicomp . ngroups ] ) <EOL> ax1 . set_xlim ( [ np . min ( minrange ) - r / <NUM_LIT> , np . max ( maxrange ) + r / <NUM_LIT> ] ) <EOL> ax1 . set_yticklabels ( np . insert ( self . groupsunique . astype ( str ) , <NUM_LIT:0> , '<STR_LIT>' ) ) <EOL> ax1 . set_yticks ( np . arange ( - <NUM_LIT:1> , len ( means ) + <NUM_LIT:1> ) ) <EOL> ax1 . set_xlabel ( xlabel if xlabel is not None else '<STR_LIT>' ) <EOL> ax1 . set_ylabel ( ylabel if ylabel is not None else '<STR_LIT>' ) <EOL> return fig <EOL> class MultiComparison ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __init__ ( self , data , groups , group_order = None ) : <EOL> if len ( data ) != len ( groups ) : <EOL> raise ValueError ( '<STR_LIT>' % ( len ( data ) , len ( groups ) ) ) <EOL> self . data = np . asarray ( data ) <EOL> self . groups = groups = np . asarray ( groups ) <EOL> if group_order is None : <EOL> self . groupsunique , self . groupintlab = np . unique ( groups , <EOL> return_inverse = True ) <EOL> else : <EOL> for grp in group_order : <EOL> if grp not in groups : <EOL> raise ValueError ( <EOL> "<STR_LIT>" % grp ) <EOL> self . groupsunique = np . array ( group_order ) <EOL> self . groupintlab = np . empty ( len ( data ) , int ) <EOL> self . groupintlab . fill ( - <NUM_LIT> ) <EOL> count = <NUM_LIT:0> <EOL> for name in self . groupsunique : <EOL> idx = np . where ( self . groups == name ) [ <NUM_LIT:0> ] <EOL> count += len ( idx ) <EOL> self . groupintlab [ idx ] = np . where ( self . groupsunique == name ) [ <NUM_LIT:0> ] <EOL> if count != data . shape [ <NUM_LIT:0> ] : <EOL> import warnings <EOL> warnings . warn ( '<STR_LIT>' + <EOL> '<STR_LIT>' ) <EOL> mask_keep = self . groupintlab != - <NUM_LIT> <EOL> self . groupintlab = self . groupintlab [ mask_keep ] <EOL> self . data = self . data [ mask_keep ] <EOL> self . groups = self . groups [ mask_keep ] <EOL> if len ( self . groupsunique ) < <NUM_LIT:2> : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> self . datali = [ data [ self . groups == k ] for k in self . groupsunique ] <EOL> self . pairindices = np . triu_indices ( len ( self . groupsunique ) , <NUM_LIT:1> ) <EOL> self . nobs = self . data . shape [ <NUM_LIT:0> ] <EOL> self . ngroups = len ( self . groupsunique ) <EOL> def getranks ( self ) : <EOL> '''<STR_LIT>''' <EOL> self . ranks = GroupsStats ( np . column_stack ( [ self . data , self . groupintlab ] ) , <EOL> useranks = True ) <EOL> self . rankdata = self . ranks . groupmeanfilter <EOL> def kruskal ( self , pairs = None , multimethod = '<STR_LIT:T>' ) : <EOL> '''<STR_LIT>''' <EOL> self . getranks ( ) <EOL> tot = self . nobs <EOL> meanranks = self . ranks . groupmean <EOL> groupnobs = self . ranks . groupnobs <EOL> f = ( tot * ( tot + <NUM_LIT:1.> ) / <NUM_LIT> ) / stats . tiecorrect ( self . rankdata ) <EOL> print ( '<STR_LIT>' ) <EOL> for i , j in zip ( * self . pairindices ) : <EOL> pdiff = np . abs ( meanranks [ i ] - meanranks [ j ] ) <EOL> se = np . sqrt ( f * np . sum ( <NUM_LIT:1.> / groupnobs [ [ i , j ] ] ) ) <EOL> Q = pdiff / se <EOL> print ( i , j , pdiff , se , pdiff / se , pdiff / se > <NUM_LIT> ) <EOL> print ( stats . norm . sf ( Q ) * <NUM_LIT:2> ) <EOL> return stats . norm . sf ( Q ) * <NUM_LIT:2> <EOL> def allpairtest ( self , testfunc , alpha = <NUM_LIT> , method = '<STR_LIT>' , pvalidx = <NUM_LIT:1> ) : <EOL> '''<STR_LIT>''' <EOL> res = [ ] <EOL> for i , j in zip ( * self . pairindices ) : <EOL> res . append ( testfunc ( self . datali [ i ] , self . datali [ j ] ) ) <EOL> res = np . array ( res ) <EOL> reject , pvals_corrected , alphacSidak , alphacBonf = multipletests ( res [ : , pvalidx ] , alpha = <NUM_LIT> , method = method ) <EOL> i1 , i2 = self . pairindices <EOL> if pvals_corrected is None : <EOL> resarr = np . array ( lzip ( self . groupsunique [ i1 ] , self . groupsunique [ i2 ] , <EOL> np . round ( res [ : , <NUM_LIT:0> ] , <NUM_LIT:4> ) , <EOL> np . round ( res [ : , <NUM_LIT:1> ] , <NUM_LIT:4> ) , <EOL> reject ) , <EOL> dtype = [ ( '<STR_LIT>' , object ) , <EOL> ( '<STR_LIT>' , object ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , np . bool8 ) ] ) <EOL> else : <EOL> resarr = np . array ( lzip ( self . groupsunique [ i1 ] , self . groupsunique [ i2 ] , <EOL> np . round ( res [ : , <NUM_LIT:0> ] , <NUM_LIT:4> ) , <EOL> np . round ( res [ : , <NUM_LIT:1> ] , <NUM_LIT:4> ) , <EOL> np . round ( pvals_corrected , <NUM_LIT:4> ) , <EOL> reject ) , <EOL> dtype = [ ( '<STR_LIT>' , object ) , <EOL> ( '<STR_LIT>' , object ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , np . bool8 ) ] ) <EOL> from statsmodels . iolib . table import SimpleTable <EOL> results_table = SimpleTable ( resarr , headers = resarr . dtype . names ) <EOL> results_table . title = ( <EOL> '<STR_LIT>' <EOL> % ( testfunc . __name__ , '<STR_LIT>' , alpha , method ) + <EOL> '<STR_LIT>' <EOL> % ( alphacSidak , alphacBonf ) ) <EOL> return results_table , ( res , reject , pvals_corrected , <EOL> alphacSidak , alphacBonf ) , resarr <EOL> def tukeyhsd ( self , alpha = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> self . groupstats = GroupsStats ( <EOL> np . column_stack ( [ self . data , self . groupintlab ] ) , <EOL> useranks = False ) <EOL> gmeans = self . groupstats . groupmean <EOL> gnobs = self . groupstats . groupnobs <EOL> var_ = np . var ( self . groupstats . groupdemean ( ) , ddof = len ( gmeans ) ) <EOL> res = tukeyhsd ( gmeans , gnobs , var_ , df = None , alpha = alpha , q_crit = None ) <EOL> resarr = np . array ( lzip ( self . groupsunique [ res [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] ] , self . groupsunique [ res [ <NUM_LIT:0> ] [ <NUM_LIT:1> ] ] , <EOL> np . round ( res [ <NUM_LIT:2> ] , <NUM_LIT:4> ) , <EOL> np . round ( res [ <NUM_LIT:4> ] [ : , <NUM_LIT:0> ] , <NUM_LIT:4> ) , <EOL> np . round ( res [ <NUM_LIT:4> ] [ : , <NUM_LIT:1> ] , <NUM_LIT:4> ) , <EOL> res [ <NUM_LIT:1> ] ) , <EOL> dtype = [ ( '<STR_LIT>' , object ) , <EOL> ( '<STR_LIT>' , object ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , np . bool8 ) ] ) <EOL> results_table = SimpleTable ( resarr , headers = resarr . dtype . names ) <EOL> results_table . title = '<STR_LIT>' + '<STR_LIT>' % alpha <EOL> return TukeyHSDResults ( self , results_table , res [ <NUM_LIT:5> ] , res [ <NUM_LIT:1> ] , res [ <NUM_LIT:2> ] , <EOL> res [ <NUM_LIT:3> ] , res [ <NUM_LIT:4> ] , res [ <NUM_LIT:6> ] , res [ <NUM_LIT:7> ] , var_ ) <EOL> def rankdata ( x ) : <EOL> '''<STR_LIT>''' <EOL> uni , intlab = np . unique ( x [ : , <NUM_LIT:0> ] , return_inverse = True ) <EOL> groupnobs = np . bincount ( intlab ) <EOL> groupxsum = np . bincount ( intlab , weights = X [ : , <NUM_LIT:0> ] ) <EOL> groupxmean = groupxsum * <NUM_LIT:1.0> / groupnobs <EOL> rankraw = x [ : , <NUM_LIT:0> ] . argsort ( ) . argsort ( ) <EOL> groupranksum = np . bincount ( intlab , weights = rankraw ) <EOL> grouprankmean = groupranksum * <NUM_LIT:1.0> / groupnobs + <NUM_LIT:1> <EOL> return grouprankmean [ intlab ] <EOL> def compare_ordered ( vals , alpha ) : <EOL> '''<STR_LIT>''' <EOL> vals = np . asarray ( vals ) <EOL> alphaf = alpha <EOL> sortind = np . argsort ( vals ) <EOL> pvals = vals [ sortind ] <EOL> sortrevind = sortind . argsort ( ) <EOL> ntests = len ( vals ) <EOL> v1 , v2 = np . triu_indices ( ntests , <NUM_LIT:1> ) <EOL> for i in range ( <NUM_LIT:4> ) : <EOL> for j in range ( <NUM_LIT:4> , i , - <NUM_LIT:1> ) : <EOL> print ( i , j ) <EOL> def varcorrection_unbalanced ( nobs_all , srange = False ) : <EOL> '''<STR_LIT>''' <EOL> nobs_all = np . asarray ( nobs_all ) <EOL> if not srange : <EOL> return ( <NUM_LIT:1.> / nobs_all ) . sum ( ) <EOL> else : <EOL> return ( <NUM_LIT:1.> / nobs_all ) . sum ( ) / len ( nobs_all ) <EOL> def varcorrection_pairs_unbalanced ( nobs_all , srange = False ) : <EOL> '''<STR_LIT>''' <EOL> n1 , n2 = np . meshgrid ( nobs_all , nobs_all ) <EOL> if not srange : <EOL> return ( <NUM_LIT:1.> / n1 + <NUM_LIT:1.> / n2 ) <EOL> else : <EOL> return ( <NUM_LIT:1.> / n1 + <NUM_LIT:1.> / n2 ) / <NUM_LIT> <EOL> def varcorrection_unequal ( var_all , nobs_all , df_all ) : <EOL> '''<STR_LIT>''' <EOL> var_all = np . asarray ( var_all ) <EOL> var_over_n = var_all * <NUM_LIT:1.> / nobs_all <EOL> varjoint = var_over_n . sum ( ) <EOL> dfjoint = varjoint ** <NUM_LIT:2> / ( var_over_n ** <NUM_LIT:2> * df_all ) . sum ( ) <EOL> return varjoint , dfjoint <EOL> def varcorrection_pairs_unequal ( var_all , nobs_all , df_all ) : <EOL> '''<STR_LIT>''' <EOL> v1 , v2 = np . meshgrid ( var_all , var_all ) <EOL> n1 , n2 = np . meshgrid ( nobs_all , nobs_all ) <EOL> df1 , df2 = np . meshgrid ( df_all , df_all ) <EOL> varjoint = v1 / n1 + v2 / n2 <EOL> dfjoint = varjoint ** <NUM_LIT:2> / ( df1 * ( v1 / n1 ) ** <NUM_LIT:2> + df2 * ( v2 / n2 ) ** <NUM_LIT:2> ) <EOL> return varjoint , dfjoint <EOL> def tukeyhsd ( mean_all , nobs_all , var_all , df = None , alpha = <NUM_LIT> , q_crit = None ) : <EOL> '''<STR_LIT>''' <EOL> mean_all = np . asarray ( mean_all ) <EOL> n_means = len ( mean_all ) <EOL> if df is None : <EOL> df = nobs_all - <NUM_LIT:1> <EOL> if np . size ( df ) == <NUM_LIT:1> : <EOL> df_total = n_means * df <EOL> df = np . ones ( n_means ) * df <EOL> else : <EOL> df_total = np . sum ( df ) <EOL> if ( np . size ( nobs_all ) == <NUM_LIT:1> ) and ( np . size ( var_all ) == <NUM_LIT:1> ) : <EOL> var_pairs = <NUM_LIT:1.> * var_all / nobs_all * np . ones ( ( n_means , n_means ) ) <EOL> elif np . size ( var_all ) == <NUM_LIT:1> : <EOL> var_pairs = var_all * varcorrection_pairs_unbalanced ( nobs_all , <EOL> srange = True ) <EOL> elif np . size ( var_all ) > <NUM_LIT:1> : <EOL> var_pairs , df_sum = varcorrection_pairs_unequal ( nobs_all , var_all , df ) <EOL> var_pairs /= <NUM_LIT> <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> meandiffs_ = mean_all - mean_all [ : , None ] <EOL> std_pairs_ = np . sqrt ( var_pairs ) <EOL> idx1 , idx2 = np . triu_indices ( n_means , <NUM_LIT:1> ) <EOL> meandiffs = meandiffs_ [ idx1 , idx2 ] <EOL> std_pairs = std_pairs_ [ idx1 , idx2 ] <EOL> st_range = np . abs ( meandiffs ) / std_pairs <EOL> df_total_ = max ( df_total , <NUM_LIT:5> ) <EOL> if q_crit is None : <EOL> q_crit = get_tukeyQcrit2 ( n_means , df_total , alpha = alpha ) <EOL> reject = st_range > q_crit <EOL> crit_int = std_pairs * q_crit <EOL> reject2 = np . abs ( meandiffs ) > crit_int <EOL> confint = np . column_stack ( ( meandiffs - crit_int , meandiffs + crit_int ) ) <EOL> return ( idx1 , idx2 ) , reject , meandiffs , std_pairs , confint , q_crit , df_total , reject2 <EOL> def simultaneous_ci ( q_crit , var , groupnobs , pairindices = None ) : <EOL> """<STR_LIT>""" <EOL> ng = len ( groupnobs ) <EOL> if pairindices is None : <EOL> pairindices = np . triu_indices ( ng , <NUM_LIT:1> ) <EOL> gvar = var / groupnobs <EOL> d12 = np . sqrt ( gvar [ pairindices [ <NUM_LIT:0> ] ] + gvar [ pairindices [ <NUM_LIT:1> ] ] ) <EOL> d = np . zeros ( ( ng , ng ) ) <EOL> d [ pairindices ] = d12 <EOL> d = d + d . conj ( ) . T <EOL> sum1 = np . sum ( d12 ) <EOL> sum2 = np . sum ( d , axis = <NUM_LIT:0> ) <EOL> if ( ng > <NUM_LIT:2> ) : <EOL> w = ( ( ng - <NUM_LIT:1.> ) * sum2 - sum1 ) / ( ( ng - <NUM_LIT:1.> ) * ( ng - <NUM_LIT> ) ) <EOL> else : <EOL> w = sum1 * np . ones ( <NUM_LIT:2> , <NUM_LIT:1> ) / <NUM_LIT> <EOL> return ( q_crit / np . sqrt ( <NUM_LIT:2> ) ) * w <EOL> def distance_st_range ( mean_all , nobs_all , var_all , df = None , triu = False ) : <EOL> '''<STR_LIT>''' <EOL> mean_all = np . asarray ( mean_all ) <EOL> n_means = len ( mean_all ) <EOL> if df is None : <EOL> df = nobs_all - <NUM_LIT:1> <EOL> if np . size ( df ) == <NUM_LIT:1> : <EOL> df_total = n_means * df <EOL> else : <EOL> df_total = np . sum ( df ) <EOL> if ( np . size ( nobs_all ) == <NUM_LIT:1> ) and ( np . size ( var_all ) == <NUM_LIT:1> ) : <EOL> var_pairs = <NUM_LIT:1.> * var_all / nobs_all * np . ones ( ( n_means , n_means ) ) <EOL> elif np . size ( var_all ) == <NUM_LIT:1> : <EOL> var_pairs = var_all * varcorrection_pairs_unbalanced ( nobs_all , <EOL> srange = True ) <EOL> elif np . size ( var_all ) > <NUM_LIT:1> : <EOL> var_pairs , df_sum = varcorrection_pairs_unequal ( nobs_all , var_all , df ) <EOL> var_pairs /= <NUM_LIT> <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> meandiffs = mean_all - mean_all [ : , None ] <EOL> std_pairs = np . sqrt ( var_pairs ) <EOL> idx1 , idx2 = np . triu_indices ( n_means , <NUM_LIT:1> ) <EOL> if triu : <EOL> meandiffs = meandiffs_ [ idx1 , idx2 ] <EOL> std_pairs = std_pairs_ [ idx1 , idx2 ] <EOL> st_range = np . abs ( meandiffs ) / std_pairs <EOL> return st_range , meandiffs , std_pairs , ( idx1 , idx2 ) <EOL> def contrast_allpairs ( nm ) : <EOL> '''<STR_LIT>''' <EOL> contr = [ ] <EOL> for i in range ( nm ) : <EOL> for j in range ( i + <NUM_LIT:1> , nm ) : <EOL> contr_row = np . zeros ( nm ) <EOL> contr_row [ i ] = <NUM_LIT:1> <EOL> contr_row [ j ] = - <NUM_LIT:1> <EOL> contr . append ( contr_row ) <EOL> return np . array ( contr ) <EOL> def contrast_all_one ( nm ) : <EOL> '''<STR_LIT>''' <EOL> contr = np . column_stack ( ( np . ones ( nm - <NUM_LIT:1> ) , - np . eye ( nm - <NUM_LIT:1> ) ) ) <EOL> return contr <EOL> def contrast_diff_mean ( nm ) : <EOL> '''<STR_LIT>''' <EOL> return np . eye ( nm ) - np . ones ( ( nm , nm ) ) / nm <EOL> def tukey_pvalues ( std_range , nm , df ) : <EOL> from statsmodels . sandbox . distributions . multivariate import mvstdtprob <EOL> contr = contrast_allpairs ( nm ) <EOL> corr = np . dot ( contr , contr . T ) / <NUM_LIT> <EOL> tstat = std_range / np . sqrt ( <NUM_LIT:2> ) * np . ones ( corr . shape [ <NUM_LIT:0> ] ) <EOL> return multicontrast_pvalues ( tstat , corr , df = df ) <EOL> def test_tukey_pvalues ( ) : <EOL> res = tukey_pvalues ( <NUM_LIT> , <NUM_LIT:3> , <NUM_LIT:16> ) <EOL> assert_almost_equal ( <NUM_LIT> , res [ <NUM_LIT:0> ] , <NUM_LIT:3> ) <EOL> assert_almost_equal ( <NUM_LIT> * np . ones ( <NUM_LIT:3> ) , res [ <NUM_LIT:1> ] , <NUM_LIT:3> ) <EOL> def multicontrast_pvalues ( tstat , tcorr , df = None , dist = '<STR_LIT:t>' , alternative = '<STR_LIT>' ) : <EOL> '''<STR_LIT>''' <EOL> from statsmodels . sandbox . distributions . multivariate import mvstdtprob <EOL> if ( df is None ) and ( dist == '<STR_LIT:t>' ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> tstat = np . asarray ( tstat ) <EOL> ntests = len ( tstat ) <EOL> cc = np . abs ( tstat ) <EOL> pval_global = <NUM_LIT:1> - mvstdtprob ( - cc , cc , tcorr , df ) <EOL> pvals = [ ] <EOL> for ti in cc : <EOL> limits = ti * np . ones ( ntests ) <EOL> pvals . append ( <NUM_LIT:1> - mvstdtprob ( - cc , cc , tcorr , df ) ) <EOL> return pval_global , np . asarray ( pvals ) <EOL> class StepDown ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __init__ ( self , vals , nobs_all , var_all , df = None ) : <EOL> self . vals = vals <EOL> self . n_vals = len ( vals ) <EOL> self . nobs_all = nobs_all <EOL> self . var_all = var_all <EOL> self . df = df <EOL> def get_crit ( self , alpha ) : <EOL> q_crit = get_tukeyQcrit ( self . n_vals , self . df , alpha = alpha ) <EOL> return q_crit * np . ones ( self . n_vals ) <EOL> def get_distance_matrix ( self ) : <EOL> '''<STR_LIT>''' <EOL> dres = distance_st_range ( self . vals , self . nobs_all , self . var_all , df = self . df ) <EOL> self . distance_matrix = dres [ <NUM_LIT:0> ] <EOL> def iter_subsets ( self , indices ) : <EOL> for ii in range ( len ( indices ) ) : <EOL> idxsub = copy . copy ( indices ) <EOL> idxsub . pop ( ii ) <EOL> yield idxsub <EOL> def check_set ( self , indices ) : <EOL> '''<STR_LIT>''' <EOL> indtup = tuple ( indices ) <EOL> if indtup in self . cache_result : <EOL> return self . cache_result [ indtup ] <EOL> else : <EOL> set_distance_matrix = self . distance_matrix [ np . asarray ( indices ) [ : , None ] , indices ] <EOL> n_elements = len ( indices ) <EOL> if np . any ( set_distance_matrix > self . crit [ n_elements - <NUM_LIT:1> ] ) : <EOL> res = True <EOL> else : <EOL> res = False <EOL> self . cache_result [ indtup ] = res <EOL> return res <EOL> def stepdown ( self , indices ) : <EOL> print ( indices ) <EOL> if self . check_set ( indices ) : <EOL> if ( len ( indices ) > <NUM_LIT:2> ) : <EOL> for subs in self . iter_subsets ( indices ) : <EOL> self . stepdown ( subs ) <EOL> else : <EOL> self . rejected . append ( tuple ( indices ) ) <EOL> else : <EOL> self . accepted . append ( tuple ( indices ) ) <EOL> return indices <EOL> def run ( self , alpha ) : <EOL> '''<STR_LIT>''' <EOL> self . cache_result = { } <EOL> self . crit = self . get_crit ( alpha ) <EOL> self . accepted = [ ] <EOL> self . rejected = [ ] <EOL> self . get_distance_matrix ( ) <EOL> self . stepdown ( lrange ( self . n_vals ) ) <EOL> return list ( set ( self . accepted ) ) , list ( set ( sd . rejected ) ) <EOL> def homogeneous_subsets ( vals , dcrit ) : <EOL> '''<STR_LIT>''' <EOL> nvals = len ( vals ) <EOL> indices_ = lrange ( nvals ) <EOL> rejected = [ ] <EOL> subsetsli = [ ] <EOL> if np . size ( dcrit ) == <NUM_LIT:1> : <EOL> dcrit = dcrit * np . ones ( ( nvals , nvals ) ) <EOL> def subsets ( vals , indices_ ) : <EOL> '''<STR_LIT>''' <EOL> i , j = ( indices_ [ <NUM_LIT:0> ] , indices_ [ - <NUM_LIT:1> ] ) <EOL> if vals [ - <NUM_LIT:1> ] - vals [ <NUM_LIT:0> ] > dcrit [ i , j ] : <EOL> rejected . append ( ( indices_ [ <NUM_LIT:0> ] , indices_ [ - <NUM_LIT:1> ] ) ) <EOL> return [ subsets ( vals [ : - <NUM_LIT:1> ] , indices_ [ : - <NUM_LIT:1> ] ) , <EOL> subsets ( vals [ <NUM_LIT:1> : ] , indices_ [ <NUM_LIT:1> : ] ) , <EOL> ( indices_ [ <NUM_LIT:0> ] , indices_ [ - <NUM_LIT:1> ] ) ] <EOL> else : <EOL> subsetsli . append ( tuple ( indices_ ) ) <EOL> return indices_ <EOL> res = subsets ( vals , indices_ ) <EOL> all_pairs = [ ( i , j ) for i in range ( nvals ) for j in range ( nvals - <NUM_LIT:1> , i , - <NUM_LIT:1> ) ] <EOL> rejs = set ( rejected ) <EOL> not_rejected = list ( set ( all_pairs ) - rejs ) <EOL> return list ( rejs ) , not_rejected , list ( set ( subsetsli ) ) , res <EOL> def set_partition ( ssli ) : <EOL> '''<STR_LIT>''' <EOL> part = [ ] <EOL> for s in sorted ( list ( set ( ssli ) ) , key = len ) [ : : - <NUM_LIT:1> ] : <EOL> s_ = set ( s ) . copy ( ) <EOL> if not any ( set ( s_ ) . intersection ( set ( t ) ) for t in part ) : <EOL> part . append ( s ) <EOL> missing = list ( set ( i for ll in ssli for i in ll ) <EOL> - set ( i for ll in part for i in ll ) ) <EOL> return part , missing <EOL> def set_remove_subs ( ssli ) : <EOL> '''<STR_LIT>''' <EOL> part = [ ] <EOL> for s in sorted ( list ( set ( ssli ) ) , key = lambda x : len ( set ( x ) ) ) [ : : - <NUM_LIT:1> ] : <EOL> if not any ( set ( s ) . issubset ( set ( t ) ) for t in part ) : <EOL> part . append ( s ) <EOL> return part <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> examples = [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , <EOL> '<STR_LIT>' , '<STR_LIT:None>' ] <EOL> if '<STR_LIT>' in examples : <EOL> x = np . array ( [ [ <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:1> ] ] ) . T + np . random . randn ( <NUM_LIT:3> , <NUM_LIT:20> ) <EOL> print ( Tukeythreegene ( * x ) ) <EOL> if ( '<STR_LIT>' in examples ) or ( '<STR_LIT>' in examples ) : <EOL> x1 = [ <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:0> , - <NUM_LIT:1> , - <NUM_LIT:1> , - <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:1> , - <NUM_LIT:1> , <NUM_LIT:1> ] <EOL> print ( lzip ( np . arange ( len ( x1 ) ) , x1 ) ) <EOL> print ( maxzero ( x1 ) ) <EOL> print ( maxzerodown ( - np . array ( x1 ) ) ) <EOL> locs = np . linspace ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:10> ) <EOL> locs = np . array ( [ <NUM_LIT:0.> ] * <NUM_LIT:6> + [ <NUM_LIT> ] * <NUM_LIT:4> ) <EOL> rvs = locs + stats . norm . rvs ( size = ( <NUM_LIT:20> , <NUM_LIT:10> ) ) <EOL> tt , tpval = stats . ttest_1samp ( rvs , <NUM_LIT:0> ) <EOL> tpval_sortind = np . argsort ( tpval ) <EOL> tpval_sorted = tpval [ tpval_sortind ] <EOL> reject = tpval_sorted < ecdf ( tpval_sorted ) * <NUM_LIT> <EOL> reject2 = max ( np . nonzero ( reject ) ) <EOL> print ( reject ) <EOL> res = np . array ( lzip ( np . round ( rvs . mean ( <NUM_LIT:0> ) , <NUM_LIT:4> ) , np . round ( tpval , <NUM_LIT:4> ) , <EOL> reject [ tpval_sortind . argsort ( ) ] ) , <EOL> dtype = [ ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , float ) , <EOL> ( '<STR_LIT>' , np . bool8 ) ] ) <EOL> print ( SimpleTable ( res , headers = res . dtype . names ) ) <EOL> print ( fdrcorrection_bak ( tpval , alpha = <NUM_LIT> ) ) <EOL> print ( reject ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( '<STR_LIT>' , multipletests ( tpval , alpha = <NUM_LIT> , method = '<STR_LIT>' ) ) <EOL> print ( '<STR_LIT>' , multipletests ( tpval , alpha = <NUM_LIT> , method = '<STR_LIT>' ) ) <EOL> print ( '<STR_LIT>' , multipletests ( tpval , alpha = <NUM_LIT> , method = '<STR_LIT>' ) ) <EOL> print ( '<STR_LIT>' , multipletests ( tpval , alpha = <NUM_LIT> , method = '<STR_LIT>' ) ) <EOL> pvals = np . array ( '<STR_LIT>' <EOL> '<STR_LIT>' . split ( ) , float ) <EOL> print ( '<STR_LIT>' ) <EOL> for meth in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] : <EOL> print ( meth ) <EOL> print ( multipletests ( pvals , alpha = <NUM_LIT> , method = meth ) ) <EOL> if '<STR_LIT>' in examples : <EOL> mcres = mcfdr ( nobs = <NUM_LIT:100> , nrepl = <NUM_LIT:1000> , ntests = <NUM_LIT:30> , ntrue = <NUM_LIT:30> , mu = <NUM_LIT:0.1> , alpha = <NUM_LIT> , rho = <NUM_LIT> ) <EOL> mcmeans = np . array ( mcres ) . mean ( <NUM_LIT:0> ) <EOL> print ( mcmeans ) <EOL> print ( mcmeans [ <NUM_LIT:0> ] / <NUM_LIT> , <NUM_LIT:1> - mcmeans [ <NUM_LIT:1> ] / <NUM_LIT> ) <EOL> print ( mcmeans [ : <NUM_LIT:4> ] , mcmeans [ - <NUM_LIT:4> : ] ) <EOL> if '<STR_LIT>' in examples : <EOL> rvsmvn = randmvn ( <NUM_LIT> , ( <NUM_LIT> , <NUM_LIT:5> ) ) <EOL> print ( np . corrcoef ( rvsmvn , rowvar = <NUM_LIT:0> ) ) <EOL> print ( rvsmvn . var ( <NUM_LIT:0> ) ) <EOL> if '<STR_LIT>' in examples : <EOL> print ( get_tukeyQcrit ( <NUM_LIT:8> , <NUM_LIT:8> , alpha = <NUM_LIT> ) , <NUM_LIT> ) <EOL> print ( get_tukeyQcrit ( <NUM_LIT:8> , <NUM_LIT:8> , alpha = <NUM_LIT> ) , <NUM_LIT> ) <EOL> if '<STR_LIT>' in examples : <EOL> X = np . array ( [ [ <NUM_LIT> , <NUM_LIT:1> ] , [ <NUM_LIT> , <NUM_LIT:1> ] , [ <NUM_LIT> , <NUM_LIT:1> ] , [ <NUM_LIT> , <NUM_LIT:1> ] , [ <NUM_LIT> , <NUM_LIT:1> ] , <EOL> [ <NUM_LIT> , <NUM_LIT:1> ] , [ <NUM_LIT> , <NUM_LIT:1> ] , [ <NUM_LIT> , <NUM_LIT:1> ] , [ <NUM_LIT> , <NUM_LIT:2> ] , [ <NUM_LIT> , <NUM_LIT:2> ] , <EOL> [ <NUM_LIT> , <NUM_LIT:2> ] , [ <NUM_LIT> , <NUM_LIT:2> ] , [ <NUM_LIT> , <NUM_LIT:2> ] , [ <NUM_LIT> , <NUM_LIT:2> ] , [ <NUM_LIT> , <NUM_LIT:2> ] , <EOL> [ <NUM_LIT> , <NUM_LIT:2> ] , [ <NUM_LIT> , <NUM_LIT:3> ] , [ <NUM_LIT> , <NUM_LIT:3> ] , [ <NUM_LIT> , <NUM_LIT:3> ] , [ <NUM_LIT> , <NUM_LIT:3> ] , <EOL> [ <NUM_LIT> , <NUM_LIT:3> ] , [ <NUM_LIT> , <NUM_LIT:3> ] , [ <NUM_LIT> , <NUM_LIT:3> ] , [ <NUM_LIT> , <NUM_LIT:4> ] , [ <NUM_LIT> , <NUM_LIT:4> ] , <EOL> [ <NUM_LIT> , <NUM_LIT:4> ] , [ <NUM_LIT> , <NUM_LIT:4> ] , [ <NUM_LIT> , <NUM_LIT:4> ] , [ <NUM_LIT> , <NUM_LIT:4> ] , [ <NUM_LIT> , <NUM_LIT:4> ] , <EOL> [ <NUM_LIT> , <NUM_LIT:4> ] ] ) <EOL> xli = [ X [ X [ : , <NUM_LIT:1> ] == k , <NUM_LIT:0> ] for k in range ( <NUM_LIT:1> , <NUM_LIT:5> ) ] <EOL> xranks = stats . rankdata ( X [ : , <NUM_LIT:0> ] ) <EOL> xranksli = [ xranks [ X [ : , <NUM_LIT:1> ] == k ] for k in range ( <NUM_LIT:1> , <NUM_LIT:5> ) ] <EOL> xnobs = np . array ( [ len ( x ) for x in xli ] ) <EOL> meanranks = [ item . mean ( ) for item in xranksli ] <EOL> sumranks = [ item . sum ( ) for item in xranksli ] <EOL> stats . norm . sf ( <NUM_LIT> ) <EOL> stats . norm . isf ( <NUM_LIT> ) <EOL> mrs = np . sort ( meanranks ) <EOL> v1 , v2 = np . triu_indices ( <NUM_LIT:4> , <NUM_LIT:1> ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( mrs [ v2 ] - mrs [ v1 ] ) <EOL> diffidx = np . argsort ( mrs [ v2 ] - mrs [ v1 ] ) [ : : - <NUM_LIT:1> ] <EOL> mrs [ v2 [ diffidx ] ] - mrs [ v1 [ diffidx ] ] <EOL> print ( '<STR_LIT>' ) <EOL> for i , j in zip ( v2 [ diffidx ] , v1 [ diffidx ] ) : <EOL> print ( i , j , stats . kruskal ( xli [ i ] , xli [ j ] ) ) <EOL> mwu , mwupval = stats . mannwhitneyu ( xli [ i ] , xli [ j ] , use_continuity = False ) <EOL> print ( mwu , mwupval * <NUM_LIT:2> , mwupval * <NUM_LIT:2> < <NUM_LIT> / <NUM_LIT> , mwupval * <NUM_LIT:2> < <NUM_LIT:0.1> / <NUM_LIT> ) <EOL> uni , intlab = np . unique ( X [ : , <NUM_LIT:0> ] , return_inverse = True ) <EOL> groupnobs = np . bincount ( intlab ) <EOL> groupxsum = np . bincount ( intlab , weights = X [ : , <NUM_LIT:0> ] ) <EOL> groupxmean = groupxsum * <NUM_LIT:1.0> / groupnobs <EOL> rankraw = X [ : , <NUM_LIT:0> ] . argsort ( ) . argsort ( ) <EOL> groupranksum = np . bincount ( intlab , weights = rankraw ) <EOL> grouprankmean = groupranksum * <NUM_LIT:1.0> / groupnobs + <NUM_LIT:1> <EOL> assert_almost_equal ( grouprankmean [ intlab ] , stats . rankdata ( X [ : , <NUM_LIT:0> ] ) , <NUM_LIT:15> ) <EOL> gs = GroupsStats ( X , useranks = True ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( gs . groupmeanfilter ) <EOL> print ( grouprankmean [ intlab ] ) <EOL> xuni , xintlab = np . unique ( X [ : , <NUM_LIT:0> ] , return_inverse = True ) <EOL> gs2 = GroupsStats ( np . column_stack ( [ X [ : , <NUM_LIT:0> ] , xintlab ] ) , useranks = True ) <EOL> rankbincount = np . bincount ( xranks . astype ( int ) ) <EOL> nties = rankbincount [ rankbincount > <NUM_LIT:1> ] <EOL> ntot = float ( len ( xranks ) ) ; <EOL> tiecorrection = <NUM_LIT:1> - ( nties ** <NUM_LIT:3> - nties ) . sum ( ) / ( ntot ** <NUM_LIT:3> - ntot ) <EOL> assert_almost_equal ( tiecorrection , stats . tiecorrect ( xranks ) , <NUM_LIT:15> ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( tiecorrection ) <EOL> print ( tiecorrect ( xranks ) ) <EOL> tot = X . shape [ <NUM_LIT:0> ] <EOL> t = <NUM_LIT> <EOL> f = ( tot * ( tot + <NUM_LIT:1.> ) / <NUM_LIT> ) - ( t / ( <NUM_LIT> * ( tot - <NUM_LIT:1.> ) ) ) <EOL> f = ( tot * ( tot + <NUM_LIT:1.> ) / <NUM_LIT> ) / stats . tiecorrect ( xranks ) <EOL> print ( '<STR_LIT>' ) <EOL> for i , j in zip ( v2 [ diffidx ] , v1 [ diffidx ] ) : <EOL> pdiff = np . abs ( meanranks [ i ] - meanranks [ j ] ) <EOL> se = np . sqrt ( f * np . sum ( <NUM_LIT:1.> / xnobs [ [ i , j ] ] ) ) <EOL> print ( i , j , pdiff , se , pdiff / se , pdiff / se > <NUM_LIT> ) <EOL> multicomp = MultiComparison ( * X . T ) <EOL> multicomp . kruskal ( ) <EOL> gsr = GroupsStats ( X , useranks = True ) <EOL> print ( '<STR_LIT>' ) <EOL> for i in range ( <NUM_LIT:10> ) : <EOL> x1 , x2 = ( np . random . randn ( <NUM_LIT:30> , <NUM_LIT:2> ) + np . array ( [ <NUM_LIT:0> , <NUM_LIT:0.5> ] ) ) . T <EOL> skw = stats . kruskal ( x1 , x2 ) <EOL> mc2 = MultiComparison ( np . r_ [ x1 , x2 ] , np . r_ [ np . zeros ( len ( x1 ) ) , np . ones ( len ( x2 ) ) ] ) <EOL> newskw = mc2 . kruskal ( ) <EOL> print ( skw , np . sqrt ( skw [ <NUM_LIT:0> ] ) , skw [ <NUM_LIT:1> ] - newskw , ( newskw / skw [ <NUM_LIT:1> ] - <NUM_LIT:1> ) * <NUM_LIT:100> ) <EOL> tablett , restt , arrtt = multicomp . allpairtest ( stats . ttest_ind ) <EOL> tablemw , resmw , arrmw = multicomp . allpairtest ( stats . mannwhitneyu ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( tablett ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( tablemw ) <EOL> tablemwhs , resmw , arrmw = multicomp . allpairtest ( stats . mannwhitneyu , method = '<STR_LIT>' ) <EOL> print ( '<STR_LIT>' ) <EOL> print ( tablemwhs ) <EOL> if '<STR_LIT>' in examples : <EOL> xli = ( np . random . randn ( <NUM_LIT> , <NUM_LIT:4> ) + np . array ( [ <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0.5> , <NUM_LIT:0.5> ] ) ) . T <EOL> xrvs , xrvsgr = catstack ( xli ) <EOL> multicompr = MultiComparison ( xrvs , xrvsgr ) <EOL> tablett , restt , arrtt = multicompr . allpairtest ( stats . ttest_ind ) <EOL> print ( tablett ) <EOL> xli = [ [ <NUM_LIT:8> , <NUM_LIT:10> , <NUM_LIT:9> , <NUM_LIT:10> , <NUM_LIT:9> ] , [ <NUM_LIT:7> , <NUM_LIT:8> , <NUM_LIT:5> , <NUM_LIT:8> , <NUM_LIT:5> ] , [ <NUM_LIT:4> , <NUM_LIT:8> , <NUM_LIT:7> , <NUM_LIT:5> , <NUM_LIT:7> ] ] <EOL> x , l = catstack ( xli ) <EOL> gs4 = GroupsStats ( np . column_stack ( [ x , l ] ) ) <EOL> print ( gs4 . groupvarwithin ( ) ) <EOL> gmeans = np . array ( [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] ) <EOL> gnobs = np . array ( [ <NUM_LIT:8> , <NUM_LIT:8> , <NUM_LIT:7> , <NUM_LIT:8> ] ) <EOL> sd = StepDown ( gmeans , gnobs , <NUM_LIT> , [ <NUM_LIT> ] ) <EOL> pvals = [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <EOL> <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] <EOL> print ( fdrcorrection0 ( pvals , alpha = <NUM_LIT> , method = '<STR_LIT>' ) ) <EOL> print ( fdrcorrection_twostage ( pvals , alpha = <NUM_LIT> , iter = False ) ) <EOL> res_tst = fdrcorrection_twostage ( pvals , alpha = <NUM_LIT> , iter = False ) <EOL> assert_almost_equal ( [ <NUM_LIT> , <NUM_LIT> ] , res_tst [ - <NUM_LIT:1> ] [ : <NUM_LIT:2> ] , <NUM_LIT:3> ) <EOL> assert_equal ( <NUM_LIT:8> , res_tst [ <NUM_LIT:0> ] . sum ( ) ) <EOL> print ( fdrcorrection_twostage ( pvals , alpha = <NUM_LIT> , iter = True ) ) <EOL> print ( '<STR_LIT>' , multipletests ( pvals , alpha = <NUM_LIT> , method = '<STR_LIT>' ) ) <EOL> test_tukey_pvalues ( ) <EOL> tukey_pvalues ( <NUM_LIT> , <NUM_LIT:3> , <NUM_LIT:16> ) </s>
