<s> """<STR_LIT>""" <EOL> from __future__ import division , print_function , absolute_import <EOL> __all__ = [ '<STR_LIT>' , <EOL> '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' , ] <EOL> import numpy as np <EOL> from numpy import float_ , int_ , ndarray <EOL> import numpy . ma as ma <EOL> from numpy . ma import MaskedArray <EOL> from . import mstats_basic as mstats <EOL> from scipy . stats . distributions import norm , beta , t , binom <EOL> def hdquantiles ( data , prob = list ( [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] ) , axis = None , var = False , ) : <EOL> """<STR_LIT>""" <EOL> def _hd_1D ( data , prob , var ) : <EOL> "<STR_LIT>" <EOL> xsorted = np . squeeze ( np . sort ( data . compressed ( ) . view ( ndarray ) ) ) <EOL> n = xsorted . size <EOL> hd = np . empty ( ( <NUM_LIT:2> , len ( prob ) ) , float_ ) <EOL> if n < <NUM_LIT:2> : <EOL> hd . flat = np . nan <EOL> if var : <EOL> return hd <EOL> return hd [ <NUM_LIT:0> ] <EOL> v = np . arange ( n + <NUM_LIT:1> ) / float ( n ) <EOL> betacdf = beta . cdf <EOL> for ( i , p ) in enumerate ( prob ) : <EOL> _w = betacdf ( v , ( n + <NUM_LIT:1> ) * p , ( n + <NUM_LIT:1> ) * ( <NUM_LIT:1> - p ) ) <EOL> w = _w [ <NUM_LIT:1> : ] - _w [ : - <NUM_LIT:1> ] <EOL> hd_mean = np . dot ( w , xsorted ) <EOL> hd [ <NUM_LIT:0> , i ] = hd_mean <EOL> hd [ <NUM_LIT:1> , i ] = np . dot ( w , ( xsorted - hd_mean ) ** <NUM_LIT:2> ) <EOL> hd [ <NUM_LIT:0> , prob == <NUM_LIT:0> ] = xsorted [ <NUM_LIT:0> ] <EOL> hd [ <NUM_LIT:0> , prob == <NUM_LIT:1> ] = xsorted [ - <NUM_LIT:1> ] <EOL> if var : <EOL> hd [ <NUM_LIT:1> , prob == <NUM_LIT:0> ] = hd [ <NUM_LIT:1> , prob == <NUM_LIT:1> ] = np . nan <EOL> return hd <EOL> return hd [ <NUM_LIT:0> ] <EOL> data = ma . array ( data , copy = False , dtype = float_ ) <EOL> p = np . array ( prob , copy = False , ndmin = <NUM_LIT:1> ) <EOL> if ( axis is None ) or ( data . ndim == <NUM_LIT:1> ) : <EOL> result = _hd_1D ( data , p , var ) <EOL> else : <EOL> if data . ndim > <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % data . ndim ) <EOL> result = ma . apply_along_axis ( _hd_1D , axis , data , p , var ) <EOL> return ma . fix_invalid ( result , copy = False ) <EOL> def hdmedian ( data , axis = - <NUM_LIT:1> , var = False ) : <EOL> """<STR_LIT>""" <EOL> result = hdquantiles ( data , [ <NUM_LIT:0.5> ] , axis = axis , var = var ) <EOL> return result . squeeze ( ) <EOL> def hdquantiles_sd ( data , prob = list ( [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] ) , axis = None ) : <EOL> """<STR_LIT>""" <EOL> def _hdsd_1D ( data , prob ) : <EOL> "<STR_LIT>" <EOL> xsorted = np . sort ( data . compressed ( ) ) <EOL> n = len ( xsorted ) <EOL> hdsd = np . empty ( len ( prob ) , float_ ) <EOL> if n < <NUM_LIT:2> : <EOL> hdsd . flat = np . nan <EOL> vv = np . arange ( n ) / float ( n - <NUM_LIT:1> ) <EOL> betacdf = beta . cdf <EOL> for ( i , p ) in enumerate ( prob ) : <EOL> _w = betacdf ( vv , ( n + <NUM_LIT:1> ) * p , ( n + <NUM_LIT:1> ) * ( <NUM_LIT:1> - p ) ) <EOL> w = _w [ <NUM_LIT:1> : ] - _w [ : - <NUM_LIT:1> ] <EOL> mx_ = np . fromiter ( [ np . dot ( w , xsorted [ np . r_ [ list ( range ( <NUM_LIT:0> , k ) ) , <EOL> list ( range ( k + <NUM_LIT:1> , n ) ) ] . astype ( int_ ) ] ) <EOL> for k in range ( n ) ] , dtype = float_ ) <EOL> mx_var = np . array ( mx_ . var ( ) , copy = False , ndmin = <NUM_LIT:1> ) * n / float ( n - <NUM_LIT:1> ) <EOL> hdsd [ i ] = float ( n - <NUM_LIT:1> ) * np . sqrt ( np . diag ( mx_var ) . diagonal ( ) / float ( n ) ) <EOL> return hdsd <EOL> data = ma . array ( data , copy = False , dtype = float_ ) <EOL> p = np . array ( prob , copy = False , ndmin = <NUM_LIT:1> ) <EOL> if ( axis is None ) : <EOL> result = _hdsd_1D ( data , p ) <EOL> else : <EOL> if data . ndim > <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % data . ndim ) <EOL> result = ma . apply_along_axis ( _hdsd_1D , axis , data , p ) <EOL> return ma . fix_invalid ( result , copy = False ) . ravel ( ) <EOL> def trimmed_mean_ci ( data , limits = ( <NUM_LIT> , <NUM_LIT> ) , inclusive = ( True , True ) , <EOL> alpha = <NUM_LIT> , axis = None ) : <EOL> """<STR_LIT>""" <EOL> data = ma . array ( data , copy = False ) <EOL> trimmed = mstats . trimr ( data , limits = limits , inclusive = inclusive , axis = axis ) <EOL> tmean = trimmed . mean ( axis ) <EOL> tstde = mstats . trimmed_stde ( data , limits = limits , inclusive = inclusive , axis = axis ) <EOL> df = trimmed . count ( axis ) - <NUM_LIT:1> <EOL> tppf = t . ppf ( <NUM_LIT:1> - alpha / <NUM_LIT> , df ) <EOL> return np . array ( ( tmean - tppf * tstde , tmean + tppf * tstde ) ) <EOL> def mjci ( data , prob = [ <NUM_LIT> , <NUM_LIT:0.5> , <NUM_LIT> ] , axis = None ) : <EOL> """<STR_LIT>""" <EOL> def _mjci_1D ( data , p ) : <EOL> data = np . sort ( data . compressed ( ) ) <EOL> n = data . size <EOL> prob = ( np . array ( p ) * n + <NUM_LIT:0.5> ) . astype ( int_ ) <EOL> betacdf = beta . cdf <EOL> mj = np . empty ( len ( prob ) , float_ ) <EOL> x = np . arange ( <NUM_LIT:1> , n + <NUM_LIT:1> , dtype = float_ ) / n <EOL> y = x - <NUM_LIT:1.> / n <EOL> for ( i , m ) in enumerate ( prob ) : <EOL> W = betacdf ( x , m - <NUM_LIT:1> , n - m ) - betacdf ( y , m - <NUM_LIT:1> , n - m ) <EOL> C1 = np . dot ( W , data ) <EOL> C2 = np . dot ( W , data ** <NUM_LIT:2> ) <EOL> mj [ i ] = np . sqrt ( C2 - C1 ** <NUM_LIT:2> ) <EOL> return mj <EOL> data = ma . array ( data , copy = False ) <EOL> if data . ndim > <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % data . ndim ) <EOL> p = np . array ( prob , copy = False , ndmin = <NUM_LIT:1> ) <EOL> if ( axis is None ) : <EOL> return _mjci_1D ( data , p ) <EOL> else : <EOL> return ma . apply_along_axis ( _mjci_1D , axis , data , p ) <EOL> def mquantiles_cimj ( data , prob = [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] , alpha = <NUM_LIT> , axis = None ) : <EOL> """<STR_LIT>""" <EOL> alpha = min ( alpha , <NUM_LIT:1> - alpha ) <EOL> z = norm . ppf ( <NUM_LIT:1> - alpha / <NUM_LIT> ) <EOL> xq = mstats . mquantiles ( data , prob , alphap = <NUM_LIT:0> , betap = <NUM_LIT:0> , axis = axis ) <EOL> smj = mjci ( data , prob , axis = axis ) <EOL> return ( xq - z * smj , xq + z * smj ) <EOL> def median_cihs ( data , alpha = <NUM_LIT> , axis = None ) : <EOL> """<STR_LIT>""" <EOL> def _cihs_1D ( data , alpha ) : <EOL> data = np . sort ( data . compressed ( ) ) <EOL> n = len ( data ) <EOL> alpha = min ( alpha , <NUM_LIT:1> - alpha ) <EOL> k = int ( binom . _ppf ( alpha / <NUM_LIT> , n , <NUM_LIT:0.5> ) ) <EOL> gk = binom . cdf ( n - k , n , <NUM_LIT:0.5> ) - binom . cdf ( k - <NUM_LIT:1> , n , <NUM_LIT:0.5> ) <EOL> if gk < <NUM_LIT:1> - alpha : <EOL> k -= <NUM_LIT:1> <EOL> gk = binom . cdf ( n - k , n , <NUM_LIT:0.5> ) - binom . cdf ( k - <NUM_LIT:1> , n , <NUM_LIT:0.5> ) <EOL> gkk = binom . cdf ( n - k - <NUM_LIT:1> , n , <NUM_LIT:0.5> ) - binom . cdf ( k , n , <NUM_LIT:0.5> ) <EOL> I = ( gk - <NUM_LIT:1> + alpha ) / ( gk - gkk ) <EOL> lambd = ( n - k ) * I / float ( k + ( n - <NUM_LIT:2> * k ) * I ) <EOL> lims = ( lambd * data [ k ] + ( <NUM_LIT:1> - lambd ) * data [ k - <NUM_LIT:1> ] , <EOL> lambd * data [ n - k - <NUM_LIT:1> ] + ( <NUM_LIT:1> - lambd ) * data [ n - k ] ) <EOL> return lims <EOL> data = ma . rray ( data , copy = False ) <EOL> if ( axis is None ) : <EOL> result = _cihs_1D ( data . compressed ( ) , alpha ) <EOL> else : <EOL> if data . ndim > <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % data . ndim ) <EOL> result = ma . apply_along_axis ( _cihs_1D , axis , data , alpha ) <EOL> return result <EOL> def compare_medians_ms ( group_1 , group_2 , axis = None ) : <EOL> """<STR_LIT>""" <EOL> ( med_1 , med_2 ) = ( ma . median ( group_1 , axis = axis ) , ma . median ( group_2 , axis = axis ) ) <EOL> ( std_1 , std_2 ) = ( mstats . stde_median ( group_1 , axis = axis ) , <EOL> mstats . stde_median ( group_2 , axis = axis ) ) <EOL> W = np . abs ( med_1 - med_2 ) / ma . sqrt ( std_1 ** <NUM_LIT:2> + std_2 ** <NUM_LIT:2> ) <EOL> return <NUM_LIT:1> - norm . cdf ( W ) <EOL> def idealfourths ( data , axis = None ) : <EOL> """<STR_LIT>""" <EOL> def _idf ( data ) : <EOL> x = data . compressed ( ) <EOL> n = len ( x ) <EOL> if n < <NUM_LIT:3> : <EOL> return [ np . nan , np . nan ] <EOL> ( j , h ) = divmod ( n / <NUM_LIT> + <NUM_LIT:5> / <NUM_LIT> , <NUM_LIT:1> ) <EOL> j = int ( j ) <EOL> qlo = ( <NUM_LIT:1> - h ) * x [ j - <NUM_LIT:1> ] + h * x [ j ] <EOL> k = n - j <EOL> qup = ( <NUM_LIT:1> - h ) * x [ k ] + h * x [ k - <NUM_LIT:1> ] <EOL> return [ qlo , qup ] <EOL> data = ma . sort ( data , axis = axis ) . view ( MaskedArray ) <EOL> if ( axis is None ) : <EOL> return _idf ( data ) <EOL> else : <EOL> return ma . apply_along_axis ( _idf , axis , data ) <EOL> def rsh ( data , points = None ) : <EOL> """<STR_LIT>""" <EOL> data = ma . array ( data , copy = False ) <EOL> if points is None : <EOL> points = data <EOL> else : <EOL> points = np . array ( points , copy = False , ndmin = <NUM_LIT:1> ) <EOL> if data . ndim != <NUM_LIT:1> : <EOL> raise AttributeError ( "<STR_LIT>" ) <EOL> n = data . count ( ) <EOL> r = idealfourths ( data , axis = None ) <EOL> h = <NUM_LIT> * ( r [ - <NUM_LIT:1> ] - r [ <NUM_LIT:0> ] ) / n ** ( <NUM_LIT:1.> / <NUM_LIT:5> ) <EOL> nhi = ( data [ : , None ] <= points [ None , : ] + h ) . sum ( <NUM_LIT:0> ) <EOL> nlo = ( data [ : , None ] < points [ None , : ] - h ) . sum ( <NUM_LIT:0> ) <EOL> return ( nhi - nlo ) / ( <NUM_LIT> * n * h ) </s>
