Python scipy.stats.kurtosis() Examples

The following are code examples for showing how to use scipy.stats.kurtosis(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: LaserTOF   Author: kyleuckert   File: test_stats.py    MIT License 6 votes vote down vote up
def test_describe_axis_none(self):
        x = np.vstack((np.ones((3, 4)), 2 * np.ones((2, 4))))

        # expected values
        e_nobs, e_minmax = (20, (1.0, 2.0))
        e_mean = 1.3999999999999999
        e_var = 0.25263157894736848
        e_skew = 0.4082482904638634
        e_kurt = -1.8333333333333333

        # actual values
        a = stats.describe(x, axis=None)

        assert_equal(a.nobs, e_nobs)
        assert_almost_equal(a.minmax, e_minmax)
        assert_almost_equal(a.mean, e_mean)
        assert_almost_equal(a.variance, e_var)
        assert_array_almost_equal(a.skewness, e_skew, decimal=13)
        assert_array_almost_equal(a.kurtosis, e_kurt, decimal=13) 
Example 2
Project: LaserTOF   Author: kyleuckert   File: test_distributions.py    MIT License 6 votes vote down vote up
def test_tukeylambda_stats_ticket_1545():
    # Some test for the variance and kurtosis of the Tukey Lambda distr.
    # See test_tukeylamdba_stats.py for more tests.

    mv = stats.tukeylambda.stats(0, moments='mvsk')
    # Known exact values:
    expected = [0, np.pi**2/3, 0, 1.2]
    assert_almost_equal(mv, expected, decimal=10)

    mv = stats.tukeylambda.stats(3.13, moments='mvsk')
    # 'expected' computed with mpmath.
    expected = [0, 0.0269220858861465102, 0, -0.898062386219224104]
    assert_almost_equal(mv, expected, decimal=10)

    mv = stats.tukeylambda.stats(0.14, moments='mvsk')
    # 'expected' computed with mpmath.
    expected = [0, 2.11029702221450250, 0, -0.02708377353223019456]
    assert_almost_equal(mv, expected, decimal=10) 
Example 3
Project: genimpro   Author: bastustrump   File: rp_extract.py    MIT License 6 votes vote down vote up
def calc_statistical_features(matrix):

    result = np.zeros((matrix.shape[0],7))
    
    result[:,0] = np.mean(matrix, axis=1)
    result[:,1] = np.var(matrix, axis=1, dtype=np.float64) # the values for variance differ between MATLAB and Numpy!
    result[:,2] = stats.skew(matrix, axis=1)
    result[:,3] = stats.kurtosis(matrix, axis=1, fisher=False) # Matlab calculates Pearson's Kurtosis
    result[:,4] = np.median(matrix, axis=1)
    result[:,5] = np.min(matrix, axis=1)
    result[:,6] = np.max(matrix, axis=1)
    
    result[np.where(np.isnan(result))] = 0
    
    return result


# PSYCHO-ACOUSTIC TRANSFORMS as individual functions


# Transform 2 Mel Scale: NOT USED by rp_extract, but included for testing purposes or for import into other programs 
Example 4
Project: QuantStudio   Author: Scorpi000   File: RiskMeasureFun.py    GNU General Public License v3.0 6 votes vote down vote up
def estimateVaR(x,alpha,method):
    N = x.shape[0]
    Avg = x.mean()
    Std = x.std()
    if method=='历史模拟':
        VaR = np.percentile(x,alpha)
        CVaR = x[x<VaR].mean()
    elif method=='正态分布':
        VaR = norm.ppf(alpha,loc=Avg,scale=Std)
        CVaR = Avg-Std/alpha*norm.pdf((VaR-Avg)/Std)
    elif method=='Cornish-Fisher':
        x = (x-Avg)/Std
        S = skew(x)
        K = kurtosis(x)-3
        q = norm.ppf(alpha)
        VaR = Avg+Std*(q+1/6*(q**2-1)*S+1/24*(q**3-3*q)*K-1/36*(2*q**3-5*q)*S**2)
        m1 = quad(lambda x:x*1/(2*np.pi)**0.5*np.exp(-x**2/2),-np.inf,q)[0]/alpha
        m2 = quad(lambda x:x**2*1/(2*np.pi)**0.5*np.exp(-x**2/2),-np.inf,q)[0]/alpha
        m3 = quad(lambda x:x**3*1/(2*np.pi)**0.5*np.exp(-x**2/2),-np.inf,q)[0]/alpha
        CVaR = Avg+Std*(m1+1/6*(m2-1)*S+1/24*(m3-3*m1)*K-1/36*(2*m3-5*m1)*S**2)
    return (VaR,CVaR) 
Example 5
Project: recruit   Author: Frank-qlu   File: test_analytics.py    Apache License 2.0 6 votes vote down vote up
def test_kurt(self, float_frame_with_na, float_frame, float_string_frame):
        from scipy.stats import kurtosis

        def alt(x):
            if len(x) < 4:
                return np.nan
            return kurtosis(x, bias=False)

        assert_stat_op_calc('kurt', alt, float_frame_with_na)
        assert_stat_op_api('kurt', float_frame, float_string_frame)

        index = MultiIndex(levels=[['bar'], ['one', 'two', 'three'], [0, 1]],
                           codes=[[0, 0, 0, 0, 0, 0],
                                  [0, 1, 2, 0, 1, 2],
                                  [0, 1, 0, 1, 0, 1]])
        df = DataFrame(np.random.randn(6, 3), index=index)

        kurt = df.kurt()
        kurt2 = df.kurt(level=0).xs('bar')
        tm.assert_series_equal(kurt, kurt2, check_names=False)
        assert kurt.name is None
        assert kurt2.name == 'bar' 
Example 6
Project: QFiPy   Author: kouzapo   File: equities.py    MIT License 6 votes vote down vote up
def descriptive_stats(self):
		closeDF = pd.read_csv('data/historical_data/' + self.quote + '.dat')['Adj Close']
		log_returns = np.log(closeDF / closeDF.shift(1)).dropna()

		desc = log_returns.describe()
		skewness = self.calc_skewness()
		kurtosis = self.calc_kurtosis()

		print('-----Descriptive Statistics for ' + self.quote + '-----')
		print('count\t', desc['count'])
		print('mean\t', round(desc['mean'], 6))
		print('std\t', round(desc['std'], 6))
		print('skew\t', round(skewness, 6))
		print('kurt\t', round(kurtosis, 6))
		print('min\t', round(desc['min'], 6))
		print('max\t', round(desc['max'], 6))
		print('25%\t', round(desc['25%'], 6))
		print('50%\t', round(desc['50%'], 6))
		print('75%\t', round(desc['75%'], 6)) 
Example 7
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_stats.py    GNU General Public License v3.0 6 votes vote down vote up
def test_kurtosis(self):
        #   sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4
        #   sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5
        #   Set flags for axis = 0 and
        #   fisher=0 (Pearson's defn of kurtosis for compatiability with Matlab)
        y = stats.kurtosis(self.testmathworks,0,fisher=0,bias=1)
        assert_approx_equal(y, 2.1658856802973,10)

        # Note that MATLAB has confusing docs for the following case
        #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
        #  kurtosis(x)  gives a biased estimate of Fisher's skewness (Pearson-3)
        #  The MATLAB docs imply that both should give Fisher's
        y = stats.kurtosis(self.testmathworks,fisher=0,bias=0)
        assert_approx_equal(y, 3.663542721189047,10)
        y = stats.kurtosis(self.testcase,0,0)
        assert_approx_equal(y,1.64) 
Example 8
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 6 votes vote down vote up
def symbols_kurtosis(self):
        from scipy.stats import kurtosis
        symbol_counts_dict = self._get_symbols_per_category()
        ## None is for checking empty, no categorical columns
        if not symbol_counts_dict:
            # return np.nan
            return 0

        symbol_counts = list(symbol_counts_dict.values())

        return kurtosis(symbol_counts, bias = False)  

    


    ##todo: Note we can evaluate symbol probabilities too.

    #----------------------------------------------------------------------
    # Kustosis related - For all non-categorical columns 
Example 9
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_analytics.py    MIT License 6 votes vote down vote up
def test_kurt(self):
        from scipy.stats import kurtosis
        alt = lambda x: kurtosis(x, bias=False)
        self._check_stat_op('kurt', alt)

        index = MultiIndex(levels=[['bar'], ['one', 'two', 'three'], [0, 1]],
                           labels=[[0, 0, 0, 0, 0, 0], [0, 1, 2, 0, 1, 2],
                                   [0, 1, 0, 1, 0, 1]])
        s = Series(np.random.randn(6), index=index)
        tm.assert_almost_equal(s.kurt(), s.kurt(level=0)['bar'])

        # test corner cases, kurt() returns NaN unless there's at least 4
        # values
        min_N = 4
        for i in range(1, min_N + 1):
            s = Series(np.ones(i))
            df = DataFrame(np.ones((i, i)))
            if i < min_N:
                assert np.isnan(s.kurt())
                assert np.isnan(df.kurt()).all()
            else:
                assert 0 == s.kurt()
                assert (df.kurt() == 0).all() 
Example 10
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_analytics.py    MIT License 6 votes vote down vote up
def test_kurt(self):
        from scipy.stats import kurtosis

        def alt(x):
            if len(x) < 4:
                return np.nan
            return kurtosis(x, bias=False)

        self._check_stat_op('kurt', alt)

        index = MultiIndex(levels=[['bar'], ['one', 'two', 'three'], [0, 1]],
                           labels=[[0, 0, 0, 0, 0, 0],
                                   [0, 1, 2, 0, 1, 2],
                                   [0, 1, 0, 1, 0, 1]])
        df = DataFrame(np.random.randn(6, 3), index=index)

        kurt = df.kurt()
        kurt2 = df.kurt(level=0).xs('bar')
        tm.assert_series_equal(kurt, kurt2, check_names=False)
        assert kurt.name is None
        assert kurt2.name == 'bar' 
Example 11
Project: vnpy_crypto   Author: birforce   File: check_moments.py    MIT License 6 votes vote down vote up
def nct_kurt_bug():
    '''test for incorrect kurtosis of nct

    D. Hogben, R. S. Pinkham, M. B. Wilk: The Moments of the Non-Central
    t-DistributionAuthor(s): Biometrika, Vol. 48, No. 3/4 (Dec., 1961),
    pp. 465-468
    '''
    from numpy.testing import assert_almost_equal
    mvsk_10_1 = (1.08372, 1.325546, 0.39993, 1.2499424941142943)
    assert_almost_equal(stats.nct.stats(10, 1, moments='mvsk'), mvsk_10_1, decimal=6)
    c1=np.array([1.08372])
    c2=np.array([.0755460, 1.25000])
    c3 = np.array([.0297802, .580566])
    c4 = np.array([0.0425458, 1.17491, 6.25])

    #calculation for df=10, for arbitrary nc
    nc = 1
    mc1 = c1.item()
    mc2 = (c2*nc**np.array([2,0])).sum()
    mc3 = (c3*nc**np.array([3,1])).sum()
    mc4 = c4=np.array([0.0425458, 1.17491, 6.25])
    mvsk_nc = mc2mvsk((mc1,mc2,mc3,mc4)) 
Example 12
Project: vnpy_crypto   Author: birforce   File: tmodel.py    MIT License 6 votes vote down vote up
def _set_start_params(self, start_params=None, use_kurtosis=False):
        if start_params is not None:
            self.start_params = start_params
        else:
            from statsmodels.regression.linear_model import OLS
            res_ols = OLS(self.endog, self.exog).fit()
            start_params = 0.1*np.ones(self.k_params)
            start_params[:self.k_vars] = res_ols.params

            if self.fix_df is False:

                if use_kurtosis:
                    kurt = stats.kurtosis(res_ols.resid)
                    df = 6./kurt + 4
                else:
                    df = 5

                start_params[-2] = df
                #TODO adjust scale for df
                start_params[-1] = np.sqrt(res_ols.scale)

            self.start_params = start_params 
Example 13
Project: vnpy_crypto   Author: birforce   File: test_analytics.py    MIT License 6 votes vote down vote up
def test_kurt(self):
        from scipy.stats import kurtosis
        alt = lambda x: kurtosis(x, bias=False)
        self._check_stat_op('kurt', alt)

        index = MultiIndex(levels=[['bar'], ['one', 'two', 'three'], [0, 1]],
                           labels=[[0, 0, 0, 0, 0, 0], [0, 1, 2, 0, 1, 2],
                                   [0, 1, 0, 1, 0, 1]])
        s = Series(np.random.randn(6), index=index)
        tm.assert_almost_equal(s.kurt(), s.kurt(level=0)['bar'])

        # test corner cases, kurt() returns NaN unless there's at least 4
        # values
        min_N = 4
        for i in range(1, min_N + 1):
            s = Series(np.ones(i))
            df = DataFrame(np.ones((i, i)))
            if i < min_N:
                assert np.isnan(s.kurt())
                assert np.isnan(df.kurt()).all()
            else:
                assert 0 == s.kurt()
                assert (df.kurt() == 0).all() 
Example 14
Project: vnpy_crypto   Author: birforce   File: test_analytics.py    MIT License 6 votes vote down vote up
def test_kurt(self):
        from scipy.stats import kurtosis

        def alt(x):
            if len(x) < 4:
                return np.nan
            return kurtosis(x, bias=False)

        self._check_stat_op('kurt', alt)

        index = MultiIndex(levels=[['bar'], ['one', 'two', 'three'], [0, 1]],
                           labels=[[0, 0, 0, 0, 0, 0],
                                   [0, 1, 2, 0, 1, 2],
                                   [0, 1, 0, 1, 0, 1]])
        df = DataFrame(np.random.randn(6, 3), index=index)

        kurt = df.kurt()
        kurt2 = df.kurt(level=0).xs('bar')
        tm.assert_series_equal(kurt, kurt2, check_names=False)
        assert kurt.name is None
        assert kurt2.name == 'bar' 
Example 15
Project: ble5-nrf52-mac   Author: tomasero   File: test_stats.py    MIT License 6 votes vote down vote up
def test_describe_axis_none(self):
        x = np.vstack((np.ones((3, 4)), 2 * np.ones((2, 4))))

        # expected values
        e_nobs, e_minmax = (20, (1.0, 2.0))
        e_mean = 1.3999999999999999
        e_var = 0.25263157894736848
        e_skew = 0.4082482904638634
        e_kurt = -1.8333333333333333

        # actual values
        a = stats.describe(x, axis=None)

        assert_equal(a.nobs, e_nobs)
        assert_almost_equal(a.minmax, e_minmax)
        assert_almost_equal(a.mean, e_mean)
        assert_almost_equal(a.variance, e_var)
        assert_array_almost_equal(a.skewness, e_skew, decimal=13)
        assert_array_almost_equal(a.kurtosis, e_kurt, decimal=13) 
Example 16
Project: ble5-nrf52-mac   Author: tomasero   File: test_distributions.py    MIT License 6 votes vote down vote up
def test_tukeylambda_stats_ticket_1545():
    # Some test for the variance and kurtosis of the Tukey Lambda distr.
    # See test_tukeylamdba_stats.py for more tests.

    mv = stats.tukeylambda.stats(0, moments='mvsk')
    # Known exact values:
    expected = [0, np.pi**2/3, 0, 1.2]
    assert_almost_equal(mv, expected, decimal=10)

    mv = stats.tukeylambda.stats(3.13, moments='mvsk')
    # 'expected' computed with mpmath.
    expected = [0, 0.0269220858861465102, 0, -0.898062386219224104]
    assert_almost_equal(mv, expected, decimal=10)

    mv = stats.tukeylambda.stats(0.14, moments='mvsk')
    # 'expected' computed with mpmath.
    expected = [0, 2.11029702221450250, 0, -0.02708377353223019456]
    assert_almost_equal(mv, expected, decimal=10) 
Example 17
Project: phoneticSimilarity   Author: ronggong   File: regression_train_predict.py    GNU Affero General Public License v3.0 5 votes vote down vote up
def phnLevelStatistics(x):
    # feature = [np.mean(x), np.std(x), np.min(x), np.max(x), np.median(x), skew(x), kurtosis(x)]
    feature = [np.mean(x), np.std(x), skew(x), kurtosis(x)]

    return feature 
Example 18
Project: LaserTOF   Author: kyleuckert   File: test_mstats_basic.py    MIT License 5 votes vote down vote up
def test_kurtosis(self):
        # Set flags for axis = 0 and fisher=0 (Pearson's definition of kurtosis
        # for compatibility with Matlab)
        y = mstats.kurtosis(self.testmathworks,0,fisher=0,bias=1)
        assert_almost_equal(y, 2.1658856802973,10)
        # Note that MATLAB has confusing docs for the following case
        #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
        #  kurtosis(x)  gives a biased estimate of Fisher's skewness (Pearson-3)
        #  The MATLAB docs imply that both should give Fisher's
        y = mstats.kurtosis(self.testmathworks,fisher=0, bias=0)
        assert_almost_equal(y, 3.663542721189047,10)
        y = mstats.kurtosis(self.testcase,0,0)
        assert_almost_equal(y,1.64)

        # test that kurtosis works on multidimensional masked arrays
        correct_2d = ma.array(np.array([-1.5, -3., -1.47247052385, 0.,
                                        -1.26979517952]),
                              mask=np.array([False, False, False, True,
                                             False], dtype=bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1),
                                  correct_2d)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row), correct_2d[i])

        correct_2d_bias_corrected = ma.array(
            np.array([-1.5, -3., -1.88988209538, 0., -0.5234638463918877]),
            mask=np.array([False, False, False, True, False], dtype=bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1,
                                                  bias=False),
                                  correct_2d_bias_corrected)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row, bias=False),
                                correct_2d_bias_corrected[i])

        # Check consistency between stats and mstats implementations
        assert_array_almost_equal_nulp(mstats.kurtosis(self.testcase_2d[2, :]),
                                       stats.kurtosis(self.testcase_2d[2, :]),
                                       nulp=4) 
Example 19
Project: LaserTOF   Author: kyleuckert   File: test_mstats_basic.py    MIT License 5 votes vote down vote up
def test_kurtosis(self):
        for n in self.get_n():
            x, y, xm, ym = self.generate_xy_sample(n)
            r = stats.kurtosis(x)
            rm = stats.mstats.kurtosis(xm)
            assert_almost_equal(r, rm, 10)

            r = stats.kurtosis(y)
            rm = stats.mstats.kurtosis(ym)
            assert_almost_equal(r, rm, 10) 
Example 20
Project: LaserTOF   Author: kyleuckert   File: test_mstats_basic.py    MIT License 5 votes vote down vote up
def test_describe_result_attributes(self):
        actual = mstats.describe(np.arange(5))
        attributes = ('nobs', 'minmax', 'mean', 'variance', 'skewness',
                      'kurtosis')
        check_named_results(actual, attributes, ma=True) 
Example 21
Project: LaserTOF   Author: kyleuckert   File: test_stats.py    MIT License 5 votes vote down vote up
def test_kurtosis(self):
        # Scalar test case
        y = stats.kurtosis(self.scalar_testcase)
        assert_approx_equal(y, -3.0)
        #   sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4
        #   sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5
        #   Set flags for axis = 0 and
        #   fisher=0 (Pearson's defn of kurtosis for compatiability with Matlab)
        y = stats.kurtosis(self.testmathworks, 0, fisher=0, bias=1)
        assert_approx_equal(y, 2.1658856802973, 10)

        # Note that MATLAB has confusing docs for the following case
        #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
        #  kurtosis(x)  gives a biased estimate of Fisher's skewness (Pearson-3)
        #  The MATLAB docs imply that both should give Fisher's
        y = stats.kurtosis(self.testmathworks, fisher=0, bias=0)
        assert_approx_equal(y, 3.663542721189047, 10)
        y = stats.kurtosis(self.testcase, 0, 0)
        assert_approx_equal(y, 1.64)

        x = np.arange(10.)
        x[9] = np.nan
        assert_equal(stats.kurtosis(x), np.nan)
        assert_almost_equal(stats.kurtosis(x, nan_policy='omit'), -1.230000)
        assert_raises(ValueError, stats.kurtosis, x, nan_policy='raise')
        assert_raises(ValueError, stats.kurtosis, x, nan_policy='foobar') 
Example 22
Project: LaserTOF   Author: kyleuckert   File: test_stats.py    MIT License 5 votes vote down vote up
def test_kurtosis_array_scalar(self):
        assert_equal(type(stats.kurtosis([1,2,3])), float) 
Example 23
Project: LaserTOF   Author: kyleuckert   File: test_stats.py    MIT License 5 votes vote down vote up
def test_describe_result_attributes(self):
        actual = stats.describe(np.arange(5))
        attributes = ('nobs', 'minmax', 'mean', 'variance', 'skewness',
                      'kurtosis')
        check_named_results(actual, attributes) 
Example 24
Project: LaserTOF   Author: kyleuckert   File: test_distributions.py    MIT License 5 votes vote down vote up
def test_moments(self):
        X = stats.skewnorm.rvs(a=4, size=int(1e6), loc=5, scale=2)
        assert_array_almost_equal([np.mean(X), np.var(X), stats.skew(X), stats.kurtosis(X)],
                                   stats.skewnorm.stats(a=4, loc=5, scale=2, moments='mvsk'),
                                   decimal=2)

        X = stats.skewnorm.rvs(a=-4, size=int(1e6), loc=5, scale=2)
        assert_array_almost_equal([np.mean(X), np.var(X), stats.skew(X), stats.kurtosis(X)],
                                   stats.skewnorm.stats(a=-4, loc=5, scale=2, moments='mvsk'),
                                   decimal=2) 
Example 25
Project: LaserTOF   Author: kyleuckert   File: test_distributions.py    MIT License 5 votes vote down vote up
def test_powerlaw_stats():
    """Test the powerlaw stats function.

    This unit test is also a regression test for ticket 1548.

    The exact values are:
    mean:
        mu = a / (a + 1)
    variance:
        sigma**2 = a / ((a + 2) * (a + 1) ** 2)
    skewness:
        One formula (see http://en.wikipedia.org/wiki/Skewness) is
            gamma_1 = (E[X**3] - 3*mu*E[X**2] + 2*mu**3) / sigma**3
        A short calculation shows that E[X**k] is a / (a + k), so gamma_1
        can be implemented as
            n = a/(a+3) - 3*(a/(a+1))*a/(a+2) + 2*(a/(a+1))**3
            d = sqrt(a/((a+2)*(a+1)**2)) ** 3
            gamma_1 = n/d
        Either by simplifying, or by a direct calculation of mu_3 / sigma**3,
        one gets the more concise formula:
            gamma_1 = -2.0 * ((a - 1) / (a + 3)) * sqrt((a + 2) / a)
    kurtosis: (See http://en.wikipedia.org/wiki/Kurtosis)
        The excess kurtosis is
            gamma_2 = mu_4 / sigma**4 - 3
        A bit of calculus and algebra (sympy helps) shows that
            mu_4 = 3*a*(3*a**2 - a + 2) / ((a+1)**4 * (a+2) * (a+3) * (a+4))
        so
            gamma_2 = 3*(3*a**2 - a + 2) * (a+2) / (a*(a+3)*(a+4)) - 3
        which can be rearranged to
            gamma_2 = 6 * (a**3 - a**2 - 6*a + 2) / (a*(a+3)*(a+4))
    """
    cases = [(1.0, (0.5, 1./12, 0.0, -1.2)),
             (2.0, (2./3, 2./36, -0.56568542494924734, -0.6))]
    for a, exact_mvsk in cases:
        mvsk = stats.powerlaw.stats(a, moments="mvsk")
        assert_array_almost_equal(mvsk, exact_mvsk) 
Example 26
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_nanops.py    MIT License 5 votes vote down vote up
def test_nankurt(self):
        from scipy.stats import kurtosis

        func1 = partial(kurtosis, fisher=True)
        func = partial(self._skew_kurt_wrap, func=func1)
        with np.errstate(invalid="ignore"):
            self.check_funs(
                nanops.nankurt,
                func,
                allow_complex=False,
                allow_str=False,
                allow_date=False,
                allow_tdelta=False,
            ) 
Example 27
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_nanops.py    MIT License 5 votes vote down vote up
def setup_method(self, method):
        # Test data + kurtosis value (computed with scipy.stats.kurtosis)
        self.samples = np.sin(np.linspace(0, 1, 200))
        self.actual_kurt = -1.2058303433799713 
Example 28
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_moments.py    MIT License 5 votes vote down vote up
def test_rolling_kurt(self, raw):
        from scipy.stats import kurtosis

        self._check_moment_func(lambda x: kurtosis(x, bias=False), name="kurt", raw=raw) 
Example 29
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_stat_reductions.py    MIT License 5 votes vote down vote up
def test_kurt(self):
        from scipy.stats import kurtosis

        string_series = tm.makeStringSeries().rename("series")

        alt = lambda x: kurtosis(x, bias=False)
        self._check_stat_op("kurt", alt, string_series)

        index = pd.MultiIndex(
            levels=[["bar"], ["one", "two", "three"], [0, 1]],
            codes=[[0, 0, 0, 0, 0, 0], [0, 1, 2, 0, 1, 2], [0, 1, 0, 1, 0, 1]],
        )
        s = Series(np.random.randn(6), index=index)
        tm.assert_almost_equal(s.kurt(), s.kurt(level=0)["bar"])

        # test corner cases, kurt() returns NaN unless there's at least 4
        # values
        min_N = 4
        for i in range(1, min_N + 1):
            s = Series(np.ones(i))
            df = DataFrame(np.ones((i, i)))
            if i < min_N:
                assert np.isnan(s.kurt())
                assert np.isnan(df.kurt()).all()
            else:
                assert 0 == s.kurt()
                assert (df.kurt() == 0).all() 
Example 30
Project: QuantStudio   Author: Scorpi000   File: StrategyTestFun.py    GNU General Public License v3.0 5 votes vote down vote up
def calcAdjustedSharpeRatio(wealth_seq, risk_free_rate=0.0, expected_return=None):
    SR = calcSharpeRatio(wealth_seq, risk_free_rate=risk_free_rate, expected_return=expected_return)
    YieldSeq = calcYieldSeq(wealth_seq)
    Skewness = skew(YieldSeq,axis=0,nan_policy='omit')
    Kurtosis = kurtosis(YieldSeq,axis=0,nan_policy='omit')
    return SR*(1+Skewness/6*SR-(Kurtosis-3)/24*SR**2)
# 计算信息比率, wealth_seq: 净值序列, array, benchmark_wealth_seq: 基准净值序列, array 
Example 31
Project: QuantStudio   Author: Scorpi000   File: StrategyTestFun.py    GNU General Public License v3.0 5 votes vote down vote up
def calcVaR(wealth_seq, alpha=0.05, method="Historical"):
    YieldSeq = calcYieldSeq(wealth_seq)
    if method=='Historical':
        VaR = np.nanpercentile(YieldSeq,alpha)
        CVaR = np.nanmean(YieldSeq[YieldSeq<VaR])
        return (VaR,CVaR)
    Avg = np.nanmean(YieldSeq)
    Std = np.nanstd(YieldSeq)
    if method=='Norm':
        VaR = norm.ppf(alpha,loc=Avg,scale=Std)
        CVaR = Avg-Std/alpha*norm.pdf((VaR-Avg)/Std)
    elif method=='Cornish-Fisher':
        x = (YieldSeq-Avg) / Std
        S = skew(x)
        K = kurtosis(x)-3
        q = norm.ppf(alpha)
        VaR = Avg+Std*(q+1/6*(q**2-1)*S+1/24*(q**3-3*q)*K-1/36*(2*q**3-5*q)*S**2)
        m1 = quad(lambda x:x*1/(2*np.pi)**0.5*np.exp(-x**2/2),-np.inf,q)[0]/alpha
        m2 = quad(lambda x:x**2*1/(2*np.pi)**0.5*np.exp(-x**2/2),-np.inf,q)[0]/alpha
        m3 = quad(lambda x:x**3*1/(2*np.pi)**0.5*np.exp(-x**2/2),-np.inf,q)[0]/alpha
        CVaR = Avg+Std*(m1+1/6*(m2-1)*S+1/24*(m3-3*m1)*K-1/36*(2*m3-5*m1)*S**2)
    else:
        VaR = np.nan
        CVaR = np.nan
    return (VaR,CVaR)
# 计算财富上升波段, wealth_seq: 净值序列, array 
Example 32
Project: QuantStudio   Author: Scorpi000   File: RiskMeasureFun.py    GNU General Public License v3.0 5 votes vote down vote up
def estimate_u(x):
    x = x[~np.isnan(x)]
    Avg = x.mean()
    K = kurtosis(x)
    while K>=3:
        Pos = np.abs(x-Avg).argmax()
        x = np.delete(x,Pos,0)
        Avg = x.mean()
        K = kurtosis(x)
    return np.abs(x).max() 
Example 33
Project: recruit   Author: Frank-qlu   File: test_window.py    Apache License 2.0 5 votes vote down vote up
def test_rolling_kurt(self):
        from scipy.stats import kurtosis
        self._check_moment_func(lambda x: kurtosis(x, bias=False),
                                name='kurt') 
Example 34
Project: recruit   Author: Frank-qlu   File: test_nanops.py    Apache License 2.0 5 votes vote down vote up
def test_nankurt(self):
        from scipy.stats import kurtosis
        func1 = partial(kurtosis, fisher=True)
        func = partial(self._skew_kurt_wrap, func=func1)
        with np.errstate(invalid='ignore'):
            self.check_funs(nanops.nankurt, func, allow_complex=False,
                            allow_str=False, allow_date=False,
                            allow_tdelta=False) 
Example 35
Project: recruit   Author: Frank-qlu   File: test_nanops.py    Apache License 2.0 5 votes vote down vote up
def setup_method(self, method):
        # Test data + kurtosis value (computed with scipy.stats.kurtosis)
        self.samples = np.sin(np.linspace(0, 1, 200))
        self.actual_kurt = -1.2058303433799713 
Example 36
Project: recruit   Author: Frank-qlu   File: test_stat_reductions.py    Apache License 2.0 5 votes vote down vote up
def test_kurt(self):
        from scipy.stats import kurtosis

        string_series = tm.makeStringSeries().rename('series')

        alt = lambda x: kurtosis(x, bias=False)
        self._check_stat_op('kurt', alt, string_series)

        index = pd.MultiIndex(
            levels=[['bar'], ['one', 'two', 'three'], [0, 1]],
            codes=[[0, 0, 0, 0, 0, 0], [0, 1, 2, 0, 1, 2], [0, 1, 0, 1, 0, 1]]
        )
        s = Series(np.random.randn(6), index=index)
        tm.assert_almost_equal(s.kurt(), s.kurt(level=0)['bar'])

        # test corner cases, kurt() returns NaN unless there's at least 4
        # values
        min_N = 4
        for i in range(1, min_N + 1):
            s = Series(np.ones(i))
            df = DataFrame(np.ones((i, i)))
            if i < min_N:
                assert np.isnan(s.kurt())
                assert np.isnan(df.kurt()).all()
            else:
                assert 0 == s.kurt()
                assert (df.kurt() == 0).all() 
Example 37
Project: QFiPy   Author: kouzapo   File: equities.py    MIT License 5 votes vote down vote up
def calc_kurtosis(self):
		"""
		This method calculates the kurtosis of the asset based on the historical returns.

		Returns:
		-------
			kurtosis: float, the kurtosis.
		"""

		kurtosis = stats.kurtosis(self.calc_log_returns())

		return kurtosis 
Example 38
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_mstats_basic.py    GNU General Public License v3.0 5 votes vote down vote up
def test_kurtosis(self):
        # Set flags for axis = 0 and fisher=0 (Pearson's definition of kurtosis
        # for compatibility with Matlab)
        y = mstats.kurtosis(self.testmathworks,0,fisher=0,bias=1)
        assert_almost_equal(y, 2.1658856802973,10)
        # Note that MATLAB has confusing docs for the following case
        #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
        #  kurtosis(x)  gives a biased estimate of Fisher's skewness (Pearson-3)
        #  The MATLAB docs imply that both should give Fisher's
        y = mstats.kurtosis(self.testmathworks,fisher=0, bias=0)
        assert_almost_equal(y, 3.663542721189047,10)
        y = mstats.kurtosis(self.testcase,0,0)
        assert_almost_equal(y,1.64)

        # test that kurtosis works on multidimensional masked arrays
        correct_2d = ma.array(np.array([-1.5, -3., -1.47247052385, 0.,
                                        -1.26979517952]),
                              mask=np.array([False, False, False, True,
                                             False], dtype=np.bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1),
                                  correct_2d)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row), correct_2d[i])

        correct_2d_bias_corrected = ma.array(
            np.array([-1.5, -3., -1.88988209538, 0., -0.5234638463918877]),
            mask=np.array([False, False, False, True, False], dtype=np.bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1,
                                                  bias=False),
                                  correct_2d_bias_corrected)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row, bias=False),
                                correct_2d_bias_corrected[i])

        # Check consistency between stats and mstats implementations
        assert_array_almost_equal_nulp(mstats.kurtosis(self.testcase_2d[2, :]),
                                       stats.kurtosis(self.testcase_2d[2, :])) 
Example 39
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_mstats_basic.py    GNU General Public License v3.0 5 votes vote down vote up
def test_kurtosis(self):
        for n in self.get_n():
            x, y, xm, ym = self.generate_xy_sample(n)
            r = stats.kurtosis(x)
            rm = stats.mstats.kurtosis(xm)
            assert_almost_equal(r, rm, 10)

            r = stats.kurtosis(y)
            rm = stats.mstats.kurtosis(ym)
            assert_almost_equal(r, rm, 10) 
Example 40
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_stats.py    GNU General Public License v3.0 5 votes vote down vote up
def test_kurtosis_array_scalar(self):
        assert_equal(type(stats.kurtosis([1,2,3])), float) 
Example 41
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_stats.py    GNU General Public License v3.0 5 votes vote down vote up
def test_describe_result_attributes(self):
        actual = stats.describe(np.arange(5))
        attributes = ('nobs', 'minmax', 'mean', 'variance', 'skewness',
                      'kurtosis')
        for i, attr in enumerate(attributes):
            assert_equal(actual[i], getattr(actual, attr)) 
Example 42
Project: sonidero   Author: ivanvladimir   File: stats.py    GNU General Public License v2.0 5 votes vote down vote up
def calculate(self,data,axis=0):
        values=[]
        for stat in self.stats:
            if stat.startswith('mean'):
                ndata=np.mean(data,axis=axis)
            elif stat.startswith('std'):
                ndata=np.std(data,axis=axis)
            elif stat.startswith('median'):
                ndata=np.median(data,axis=axis)
            elif stat.startswith('q1'):
                ndata=[scoreatpercentile(data[:,col],25)
                        for col in range(data.shape[1]) ]
            elif stat.startswith('q2'):
                ndata=[scoreatpercentile(data[:,col],50)
                        for col in range(data.shape[1]) ]
            elif stat.startswith('q3'):
                ndata=[scoreatpercentile(data[:,col],75)
                        for col in range(data.shape[1]) ]
            elif stat.startswith('skew'):
                ndata=[skew(data[:,col])
                        for col in range(data.shape[1]) ]
            elif stat.startswith('kurtosis'):
                ndata=[kurtosis(data[:,col])
                        for col in range(data.shape[1]) ]
            elif stat.startswith('max'):
                ndata=np.amax(data,axis=axis)
            elif stat.startswith('min'):
                ndata=np.amin(data,axis=axis)
            values.append(ndata)
        return np.hstack(values) 
Example 43
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def corr_with_dependent_abs_75p(self):
        """ 75p absolute pearson correlation with dependent variable
        returns np.nan for classificaiton problems. Uses df_encoded
        ie dataframe with categorical columns encoded automatically.
        """
        if self.prediction_type == 'classification':
            return np.nan
        else:
            abs_corr_with_dependent = self._get_corr_with_dependent().abs()
            return np.nanpercentile(abs_corr_with_dependent, 75)

    #todo: try kurtosis and skew for correl values without abs. 
Example 44
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def corr_with_dependent_abs_kurtosis(self):
        """ kurtosis of absolute pearson correlation with dependent variable
        returns np.nan for classificaiton problems. Uses df_encoded
        ie dataframe with categorical columns encoded automatically.
        """
        from scipy.stats import kurtosis
        if self.prediction_type == 'classification':
            return np.nan
        else:
            abs_corr_with_dependent = self._get_corr_with_dependent().abs()
            return kurtosis(abs_corr_with_dependent, bias = False) 
Example 45
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def class_prob_median(self):
        if self.prediction_type=='regression':
            return np.nan
        else:
            class_probablities = self._get_class_probablity()
            return class_probablities.median()

    #todo: add kurtosis and skew here too. Classes will be usually less, so 
    #may not make sense.


    #----------------------------------------------------------------------
    # Symbols related - All the categorical columns 
Example 46
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def _get_kurtosis_per_num_column(self):
        """Sets an dictionary with kurtosis per numerical column"""

        if self.kurtosis_dict == None:
            self.kurtosis_dict = {}
            numerical_cols = list(set(self.independent_cols) - set(self.categorical_cols)) 
            for column in numerical_cols:
                self.kurtosis_dict[column] = kurtosis(self.df[column].dropna(), bias = False)
            return self.kurtosis_dict
        else:
            return self.kurtosis_dict 
Example 47
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def kurtosis_median(self):
        """ Median kurtosis per columns """

        kurtosis_dict = self._get_kurtosis_per_num_column()
        ## None is for checking empty, no categorical columns
        
        if not kurtosis_dict:
            # return np.nan
            return 0
        
        kurtosisses = list(kurtosis_dict.values())

        return np.nanmedian(kurtosisses) 
Example 48
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def kurtosis_min(self):
        """ Min kurtosis per columns """

        kurtosis_dict = self._get_kurtosis_per_num_column()
        ## None is for checking empty, no categorical columns
        
        if not kurtosis_dict:
            # return np.nan
            return 0
        
        kurtosisses = list(kurtosis_dict.values())

        return np.min(kurtosisses) 
Example 49
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def kurtosis_max(self):
        """ Max kurtosis per columns """

        kurtosis_dict = self._get_kurtosis_per_num_column()
        ## None is for checking empty, no categorical columns
        
        if not kurtosis_dict:
            # return np.nan
            return 0
        
        kurtosisses = list(kurtosis_dict.values())

        return np.max(kurtosisses) 
Example 50
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def kurtosis_std(self):
        """ STD of kurtosis per columns """

        kurtosis_dict = self._get_kurtosis_per_num_column()
        ## None is for checking empty, no categorical columns
        
        if not kurtosis_dict:
            # return np.nan
            return 0
        
        kurtosisses = list(kurtosis_dict.values())

        return np.nanstd(kurtosisses) 
Example 51
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def kurtosis_kurtosis(self):
        """ Kurtosis of kurtosis per columns """

        kurtosis_dict = self._get_kurtosis_per_num_column()
        ## None is for checking empty, no categorical columns
        
        if not kurtosis_dict:
            # return np.nan
            return 0
        
        kurtosisses = list(kurtosis_dict.values())

        return kurtosis(kurtosisses, bias = False) 
Example 52
Project: pennai   Author: EpistasisLab   File: dataset_describe.py    GNU General Public License v3.0 5 votes vote down vote up
def skew_kurtosis(self):
        """ kurtosis of skew in all numerical columns """

        skew_dict = self._get_skew_per_num_column()
        ## None is for checking empty, no categorical columns
        
        if not skew_dict:
            # return np.nan
            return 0
        
        skews = list(skew_dict.values())

        return kurtosis(skews, bias = False) 
Example 53
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_window.py    MIT License 5 votes vote down vote up
def test_rolling_kurt(self):
        from scipy.stats import kurtosis
        self._check_moment_func(lambda x: kurtosis(x, bias=False),
                                name='kurt') 
Example 54
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_nanops.py    MIT License 5 votes vote down vote up
def test_nankurt(self):
        from scipy.stats import kurtosis
        func1 = partial(kurtosis, fisher=True)
        func = partial(self._skew_kurt_wrap, func=func1)
        with np.errstate(invalid='ignore'):
            self.check_funs(nanops.nankurt, func, allow_complex=False,
                            allow_str=False, allow_date=False,
                            allow_tdelta=False) 
Example 55
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_nanops.py    MIT License 5 votes vote down vote up
def setup_method(self, method):
        # Test data + kurtosis value (computed with scipy.stats.kurtosis)
        self.samples = np.sin(np.linspace(0, 1, 200))
        self.actual_kurt = -1.2058303433799713 
Example 56
Project: vnpy_crypto   Author: birforce   File: descriptivestats.py    MIT License 5 votes vote down vote up
def _kurtosis(a):
    '''wrapper for scipy.stats.kurtosis that returns nan instead of raising Error

    missing options
    '''
    try:
        res = stats.kurtosis(a)
    except ValueError:
        res = np.nan
    return res 
Example 57
Project: vnpy_crypto   Author: birforce   File: stattools.py    MIT License 5 votes vote down vote up
def _kr3(y, alpha=5.0, beta=50.0):
    """
    KR3 estimator from Kim & White

    Parameters
    ----------
    y : array-like, 1-d
    alpha : float, optional
        Lower cut-off for measuring expectation in tail.
    beta :  float, optional
        Lower cut-off for measuring expectation in center.

    Returns
    -------
    kr3 : float
        Robust kurtosis estimator based on standardized lower- and upper-tail
        expected values

    Notes
    -----
    .. [*] Tae-Hwan Kim and Halbert White, "On more robust estimation of
       skewness and kurtosis," Finance Research Letters, vol. 1, pp. 56-73,
       March 2004.
    """
    perc = (alpha, 100.0 - alpha, beta, 100.0 - beta)
    lower_alpha, upper_alpha, lower_beta, upper_beta = np.percentile(y, perc)
    l_alpha = np.mean(y[y < lower_alpha])
    u_alpha = np.mean(y[y > upper_alpha])

    l_beta = np.mean(y[y < lower_beta])
    u_beta = np.mean(y[y > upper_beta])

    return (u_alpha - l_alpha) / (u_beta - l_beta) 
Example 58
Project: vnpy_crypto   Author: birforce   File: check_moments.py    MIT License 5 votes vote down vote up
def mnc2mvsk(args):
    '''convert central moments to mean, variance, skew, kurtosis
    '''
    #convert four non-central moments to central moments
    mnc, mnc2, mnc3, mnc4 = args
    mc = mnc
    mc2 = mnc2 - mnc*mnc
    mc3 = mnc3 - (3*mc*mc2+mc**3) # 3rd central moment
    mc4 = mnc4 - (4*mc*mc3+6*mc*mc*mc2+mc**4)
    return mc2mvsk((mc, mc2, mc3, mc4)) 
Example 59
Project: vnpy_crypto   Author: birforce   File: descriptive.py    MIT License 5 votes vote down vote up
def _opt_kurt(self, nuis_params):
        """
        Called by test_kurt.  This function is optimized over
        nuisance parameters mu and sigma

        Parameters
        ----------
        nuis_params : 1darray
            An array with a nuisance mean and variance parameter

        Returns
        -------
        llr : float
            The log likelihood ratio of a pre-speified kurtosis holding the
            nuisance parameters constant
        """
        endog = self.endog
        nobs = self.nobs
        mu_data = endog - nuis_params[0]
        sig_data = ((endog - nuis_params[0]) ** 2) - nuis_params[1]
        kurt_data = (((((endog - nuis_params[0]) ** 4) / \
                    (nuis_params[1] ** 2))) - 3) - self.kurt0
        est_vect = np.column_stack((mu_data, sig_data, kurt_data))
        eta_star = self._modif_newton(np.array([1. / nobs,
                                               1. / nobs,
                                               1. / nobs]), est_vect,
                                               np.ones(nobs) * (1. / nobs))
        denom = 1 + np.dot(eta_star, est_vect.T)
        self.new_weights = 1. / nobs * 1. / denom
        llr = np.sum(np.log(nobs * self.new_weights))
        return -2 * llr 
Example 60
Project: vnpy_crypto   Author: birforce   File: descriptive.py    MIT License 5 votes vote down vote up
def _opt_skew_kurt(self, nuis_params):
        """
        Called by test_joint_skew_kurt.  This function is optimized over
        nuisance parameters mu and sigma

        Parameters
        -----------
        nuis_params : 1darray
            An array with a nuisance mean and variance parameter

        Returns
        ------
        llr : float
            The log likelihood ratio of a pre-speified skewness and
            kurtosis holding the nuisance parameters constant.
        """
        endog = self.endog
        nobs = self.nobs
        mu_data = endog - nuis_params[0]
        sig_data = ((endog - nuis_params[0]) ** 2) - nuis_params[1]
        skew_data = ((((endog - nuis_params[0]) ** 3) / \
                    (nuis_params[1] ** 1.5))) - self.skew0
        kurt_data = (((((endog - nuis_params[0]) ** 4) / \
                    (nuis_params[1] ** 2))) - 3) - self.kurt0
        est_vect = np.column_stack((mu_data, sig_data, skew_data, kurt_data))
        eta_star = self._modif_newton(np.array([1. / nobs,
                                               1. / nobs,
                                               1. / nobs,
                                               1. / nobs]), est_vect,
                                               np.ones(nobs) * (1. / nobs))
        denom = 1. + np.dot(eta_star, est_vect.T)
        self.new_weights = 1. / nobs * 1. / denom
        llr = np.sum(np.log(nobs * self.new_weights))
        return -2 * llr 
Example 61
Project: vnpy_crypto   Author: birforce   File: descriptive.py    MIT License 5 votes vote down vote up
def test_kurt(self, kurt0, return_weights=False):
        """
        Returns -2 x log-likelihood and the p-value for the hypothesized
        kurtosis.

        Parameters
        ----------
        kurt0 : float
            Kurtosis value to be tested

        return_weights : bool
            If True, function also returns the weights that
            maximize the likelihood ratio. Default is False.

        Returns
        -------
        test_results : tuple
            The log-likelihood ratio and p-value of kurt0
        """
        self.kurt0 = kurt0
        start_nuisance = np.array([self.endog.mean(),
                                       self.endog.var()])

        llr = optimize.fmin_powell(self._opt_kurt, start_nuisance,
                                     full_output=1, disp=0)[1]
        p_val = chi2.sf(llr, 1)
        if return_weights:
            return  llr, p_val, self.new_weights.T
        return llr, p_val 
Example 62
Project: vnpy_crypto   Author: birforce   File: descriptive.py    MIT License 5 votes vote down vote up
def test_joint_skew_kurt(self, skew0, kurt0, return_weights=False):
        """
        Returns - 2 x log-likelihood and the p-value for the joint
        hypothesis test for skewness and kurtosis

        Parameters
        ----------
        skew0 : float
            Skewness value to be tested
        kurt0 : float
            Kurtosis value to be tested

        return_weights : bool
            If True, function also returns the weights that
            maximize the likelihood ratio. Default is False.

        Returns
        -------
        test_results : tuple
            The log-likelihood ratio and p-value  of the joint hypothesis test.
        """
        self.skew0 = skew0
        self.kurt0 = kurt0
        start_nuisance = np.array([self.endog.mean(),
                                       self.endog.var()])

        llr = optimize.fmin_powell(self._opt_skew_kurt, start_nuisance,
                                     full_output=1, disp=0)[1]
        p_val = chi2.sf(llr, 2)
        if return_weights:
            return llr, p_val, self.new_weights.T
        return llr, p_val 
Example 63
Project: vnpy_crypto   Author: birforce   File: test_window.py    MIT License 5 votes vote down vote up
def test_rolling_kurt(self):
        from scipy.stats import kurtosis
        self._check_moment_func(lambda x: kurtosis(x, bias=False),
                                name='kurt') 
Example 64
Project: vnpy_crypto   Author: birforce   File: test_nanops.py    MIT License 5 votes vote down vote up
def setup_method(self, method):
        # Test data + kurtosis value (computed with scipy.stats.kurtosis)
        self.samples = np.sin(np.linspace(0, 1, 200))
        self.actual_kurt = -1.2058303433799713 
Example 65
Project: ironn   Author: mfqqio   File: get_roughness.py    GNU General Public License v3.0 5 votes vote down vote up
def get_roughness_feat(self, aoi):
        """
        Retrieve roughness features (skewness, kurtosis, average pixel) of each given aoi

        Parameters:
        -----------
        aoi: numpy.ndarray (dimension: height * width * channel)
            3D arrays that represent each subimage (i.e. each rock type)

        Returns:
        --------
        pixels: list
            a list that contains pixels on 3 channels (Blue, Green, Red)
        feats: list
            a list that contains 3 roughness features on 3 channels with the following order:
            ["SkewnessBlue","KurtosisBlue","MeanPixelBlue",
            "SkewnessGreen","KurtosisGreen","MeanPixelGreen",
             "SkewnessRed","KurtosisRed","MeanPixelRed"]
        """
        feats = []
        pixels = []
        for col in range(3): # 3 colour channels
            img_flatten = aoi[:,:,col].reshape(aoi.shape[0]*aoi.shape[1],1)
            img_nonzero = img_flatten[img_flatten.nonzero()]
            feats.append(skew(img_nonzero))
            feats.append(kurtosis(img_nonzero))
            feats.append(np.mean(img_nonzero))
            pixels.append(img_nonzero)
        return pixels, feats 
Example 66
Project: ble5-nrf52-mac   Author: tomasero   File: test_mstats_basic.py    MIT License 5 votes vote down vote up
def test_kurtosis(self):
        # Set flags for axis = 0 and fisher=0 (Pearson's definition of kurtosis
        # for compatibility with Matlab)
        y = mstats.kurtosis(self.testmathworks, 0, fisher=0, bias=1)
        assert_almost_equal(y, 2.1658856802973, 10)
        # Note that MATLAB has confusing docs for the following case
        #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
        #  kurtosis(x) gives a biased estimate of Fisher's skewness (Pearson-3)
        #  The MATLAB docs imply that both should give Fisher's
        y = mstats.kurtosis(self.testmathworks, fisher=0, bias=0)
        assert_almost_equal(y, 3.663542721189047, 10)
        y = mstats.kurtosis(self.testcase, 0, 0)
        assert_almost_equal(y, 1.64)

        # test that kurtosis works on multidimensional masked arrays
        correct_2d = ma.array(np.array([-1.5, -3., -1.47247052385, 0.,
                                        -1.26979517952]),
                              mask=np.array([False, False, False, True,
                                             False], dtype=bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1),
                                  correct_2d)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row), correct_2d[i])

        correct_2d_bias_corrected = ma.array(
            np.array([-1.5, -3., -1.88988209538, 0., -0.5234638463918877]),
            mask=np.array([False, False, False, True, False], dtype=bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1,
                                                  bias=False),
                                  correct_2d_bias_corrected)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row, bias=False),
                                correct_2d_bias_corrected[i])

        # Check consistency between stats and mstats implementations
        assert_array_almost_equal_nulp(mstats.kurtosis(self.testcase_2d[2, :]),
                                       stats.kurtosis(self.testcase_2d[2, :]),
                                       nulp=4) 
Example 67
Project: ble5-nrf52-mac   Author: tomasero   File: test_mstats_basic.py    MIT License 5 votes vote down vote up
def test_kurtosis(self):
        for n in self.get_n():
            x, y, xm, ym = self.generate_xy_sample(n)
            r = stats.kurtosis(x)
            rm = stats.mstats.kurtosis(xm)
            assert_almost_equal(r, rm, 10)

            r = stats.kurtosis(y)
            rm = stats.mstats.kurtosis(ym)
            assert_almost_equal(r, rm, 10) 
Example 68
Project: ble5-nrf52-mac   Author: tomasero   File: test_mstats_basic.py    MIT License 5 votes vote down vote up
def test_describe_result_attributes(self):
        actual = mstats.describe(np.arange(5))
        attributes = ('nobs', 'minmax', 'mean', 'variance', 'skewness',
                      'kurtosis')
        check_named_results(actual, attributes, ma=True) 
Example 69
Project: ble5-nrf52-mac   Author: tomasero   File: test_stats.py    MIT License 5 votes vote down vote up
def test_kurtosis(self):
        # Scalar test case
        y = stats.kurtosis(self.scalar_testcase)
        assert_approx_equal(y, -3.0)
        #   sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4
        #   sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5
        #   Set flags for axis = 0 and
        #   fisher=0 (Pearson's defn of kurtosis for compatibility with Matlab)
        y = stats.kurtosis(self.testmathworks, 0, fisher=0, bias=1)
        assert_approx_equal(y, 2.1658856802973, 10)

        # Note that MATLAB has confusing docs for the following case
        #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
        #  kurtosis(x)  gives a biased estimate of Fisher's skewness (Pearson-3)
        #  The MATLAB docs imply that both should give Fisher's
        y = stats.kurtosis(self.testmathworks, fisher=0, bias=0)
        assert_approx_equal(y, 3.663542721189047, 10)
        y = stats.kurtosis(self.testcase, 0, 0)
        assert_approx_equal(y, 1.64)

        x = np.arange(10.)
        x[9] = np.nan
        assert_equal(stats.kurtosis(x), np.nan)
        assert_almost_equal(stats.kurtosis(x, nan_policy='omit'), -1.230000)
        assert_raises(ValueError, stats.kurtosis, x, nan_policy='raise')
        assert_raises(ValueError, stats.kurtosis, x, nan_policy='foobar') 
Example 70
Project: ble5-nrf52-mac   Author: tomasero   File: test_stats.py    MIT License 5 votes vote down vote up
def test_kurtosis_array_scalar(self):
        assert_equal(type(stats.kurtosis([1,2,3])), float) 
Example 71
Project: ble5-nrf52-mac   Author: tomasero   File: test_stats.py    MIT License 5 votes vote down vote up
def test_describe_result_attributes(self):
        actual = stats.describe(np.arange(5))
        attributes = ('nobs', 'minmax', 'mean', 'variance', 'skewness',
                      'kurtosis')
        check_named_results(actual, attributes) 
Example 72
Project: ble5-nrf52-mac   Author: tomasero   File: test_distributions.py    MIT License 5 votes vote down vote up
def test_moments(self):
        X = stats.skewnorm.rvs(a=4, size=int(1e6), loc=5, scale=2)
        expected = [np.mean(X), np.var(X), stats.skew(X), stats.kurtosis(X)]
        computed = stats.skewnorm.stats(a=4, loc=5, scale=2, moments='mvsk')
        assert_array_almost_equal(computed, expected, decimal=2)

        X = stats.skewnorm.rvs(a=-4, size=int(1e6), loc=5, scale=2)
        expected = [np.mean(X), np.var(X), stats.skew(X), stats.kurtosis(X)]
        computed = stats.skewnorm.stats(a=-4, loc=5, scale=2, moments='mvsk')
        assert_array_almost_equal(computed, expected, decimal=2) 
Example 73
Project: ble5-nrf52-mac   Author: tomasero   File: test_distributions.py    MIT License 5 votes vote down vote up
def test_powerlaw_stats():
    """Test the powerlaw stats function.

    This unit test is also a regression test for ticket 1548.

    The exact values are:
    mean:
        mu = a / (a + 1)
    variance:
        sigma**2 = a / ((a + 2) * (a + 1) ** 2)
    skewness:
        One formula (see https://en.wikipedia.org/wiki/Skewness) is
            gamma_1 = (E[X**3] - 3*mu*E[X**2] + 2*mu**3) / sigma**3
        A short calculation shows that E[X**k] is a / (a + k), so gamma_1
        can be implemented as
            n = a/(a+3) - 3*(a/(a+1))*a/(a+2) + 2*(a/(a+1))**3
            d = sqrt(a/((a+2)*(a+1)**2)) ** 3
            gamma_1 = n/d
        Either by simplifying, or by a direct calculation of mu_3 / sigma**3,
        one gets the more concise formula:
            gamma_1 = -2.0 * ((a - 1) / (a + 3)) * sqrt((a + 2) / a)
    kurtosis: (See https://en.wikipedia.org/wiki/Kurtosis)
        The excess kurtosis is
            gamma_2 = mu_4 / sigma**4 - 3
        A bit of calculus and algebra (sympy helps) shows that
            mu_4 = 3*a*(3*a**2 - a + 2) / ((a+1)**4 * (a+2) * (a+3) * (a+4))
        so
            gamma_2 = 3*(3*a**2 - a + 2) * (a+2) / (a*(a+3)*(a+4)) - 3
        which can be rearranged to
            gamma_2 = 6 * (a**3 - a**2 - 6*a + 2) / (a*(a+3)*(a+4))
    """
    cases = [(1.0, (0.5, 1./12, 0.0, -1.2)),
             (2.0, (2./3, 2./36, -0.56568542494924734, -0.6))]
    for a, exact_mvsk in cases:
        mvsk = stats.powerlaw.stats(a, moments="mvsk")
        assert_array_almost_equal(mvsk, exact_mvsk) 
Example 74
Project: Computable   Author: ktraunmueller   File: test_moments.py    MIT License 5 votes vote down vote up
def test_rolling_kurt(self):
        try:
            from scipy.stats import kurtosis
        except ImportError:
            raise nose.SkipTest('no scipy')
        self._check_moment_func(mom.rolling_kurt,
                                lambda x: kurtosis(x, bias=False)) 
Example 75
Project: Computable   Author: ktraunmueller   File: test_mstats_basic.py    MIT License 5 votes vote down vote up
def test_kurtosis(self):
        #    sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4
        #    sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5
        #    Set flags for axis = 0 and
        #    fisher=0 (Pearson's definition of kurtosis for compatibility with Matlab)
        y = mstats.kurtosis(self.testmathworks,0,fisher=0,bias=1)
        assert_almost_equal(y, 2.1658856802973,10)
        # Note that MATLAB has confusing docs for the following case
        #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
        #  kurtosis(x)  gives a biased estimate of Fisher's skewness (Pearson-3)
        #  The MATLAB docs imply that both should give Fisher's
        y = mstats.kurtosis(self.testmathworks,fisher=0, bias=0)
        assert_almost_equal(y, 3.663542721189047,10)
        y = mstats.kurtosis(self.testcase,0,0)
        assert_almost_equal(y,1.64)

        # test that kurtosis works on multidimensional masked arrays
        correct_2d = ma.array(np.array([-1.5, -3., -1.47247052385, 0.,
                                        -1.26979517952]),
                              mask=np.array([False, False, False, True,
                                             False], dtype=np.bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1),
                                  correct_2d)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row), correct_2d[i])

        correct_2d_bias_corrected = ma.array(
            np.array([-1.5, -3., -1.88988209538, 0., -0.5234638463918877]),
            mask=np.array([False, False, False, True, False], dtype=np.bool))
        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1,
                                                  bias=False),
                                  correct_2d_bias_corrected)
        for i, row in enumerate(self.testcase_2d):
            assert_almost_equal(mstats.kurtosis(row, bias=False),
                                correct_2d_bias_corrected[i])

        # Check consistency between stats and mstats implementations
        assert_array_almost_equal_nulp(mstats.kurtosis(self.testcase_2d[2, :]),
                                       stats.kurtosis(self.testcase_2d[2, :])) 
Example 76
Project: vnpy_crypto   Author: birforce   File: descriptivestats.py    MIT License 4 votes vote down vote up
def __init__(self, dataset):
        self.dataset = dataset

        #better if this is initially a list to define order, or use an
        # ordered dict. First position is the function
        # Second position is the tuple/list of column names/numbers
        # third is are the results in order of the columns
        self.univariate = dict(
            obs = [len, None, None],
            mean = [np.mean, None, None],
            std = [np.std, None, None],
            min = [np.min, None, None],
            max = [np.max, None, None],
            ptp = [np.ptp, None, None],
            var = [np.var, None, None],
            mode_val = [self._mode_val, None, None],
            mode_bin = [self._mode_bin, None, None],
            median = [np.median, None, None],
            skew = [stats.skew, None, None],
            uss = [lambda x: np.sum(np.asarray(x)**2, axis=0), None, None],
            kurtosis = [stats.kurtosis, None, None],
            percentiles = [self._percentiles, None, None],
            #BUG: not single value
            #sign_test_M = [self.sign_test_m, None, None],
            #sign_test_P = [self.sign_test_p, None, None]
        )

        #TODO: Basic stats for strings
        #self.strings = dict(
            #unique = [np.unique, None, None],
            #number_uniq = [len(
            #most = [
            #least = [

        #TODO: Multivariate
        #self.multivariate = dict(
            #corrcoef(x[, y, rowvar, bias]),
            #cov(m[, y, rowvar, bias]),
            #histogram2d(x, y[, bins, range, normed, weights])
            #)
        self._arraytype = None
        self._columns_list = None 
Example 77
Project: vnpy_crypto   Author: birforce   File: stattools.py    MIT License 4 votes vote down vote up
def jarque_bera(resids, axis=0):
    r"""
    Calculates the Jarque-Bera test for normality

    Parameters
    -----------
    data : array-like
        Data to test for normality
    axis : int, optional
        Axis to use if data has more than 1 dimension. Default is 0

    Returns
    -------
    JB : float or array
        The Jarque-Bera test statistic
    JBpv : float or array
        The pvalue of the test statistic
    skew : float or array
        Estimated skewness of the data
    kurtosis : float or array
        Estimated kurtosis of the data

    Notes
    -----
    Each output returned has 1 dimension fewer than data


    The Jarque-Bera test statistic tests the null that the data is normally
    distributed against an alternative that the data follow some other
    distribution. The test statistic is based on two moments of the data,
    the skewness, and the kurtosis, and has an asymptotic :math:`\chi^2_2`
    distribution.

    The test statistic is defined

    .. math:: JB = n(S^2/6+(K-3)^2/24)

    where n is the number of data points, S is the sample skewness, and K is
    the sample kurtosis of the data.
    """
    resids = np.asarray(resids)
    # Calculate residual skewness and kurtosis
    skew = stats.skew(resids, axis=axis)
    kurtosis = 3 + stats.kurtosis(resids, axis=axis)

    # Calculate the Jarque-Bera test for normality
    n = resids.shape[axis]
    jb = (n / 6.) * (skew ** 2 + (1 / 4.) * (kurtosis - 3) ** 2)
    jb_pv = stats.chi2.sf(jb, 2)

    return jb, jb_pv, skew, kurtosis 
Example 78
Project: vnpy_crypto   Author: birforce   File: stattools.py    MIT License 4 votes vote down vote up
def expected_robust_kurtosis(ab=(5.0, 50.0), dg=(2.5, 25.0)):
    """
    Calculates the expected value of the robust kurtosis measures in Kim and
    White assuming the data are normally distributed.

    Parameters
    ----------
    ab: iterable, optional
        Contains 100*(alpha, beta) in the kr3 measure where alpha is the tail
        quantile cut-off for measuring the extreme tail and beta is the central
        quantile cutoff for the standardization of the measure
    db: iterable, optional
        Contains 100*(delta, gamma) in the kr4 measure where delta is the tail
        quantile for measuring extreme values and gamma is the central quantile
        used in the the standardization of the measure

    Returns
    -------
    ekr: array, 4-element
        Contains the expected values of the 4 robust kurtosis measures

    Notes
    -----
    See `robust_kurtosis` for definitions of the robust kurtosis measures
    """

    alpha, beta = ab
    delta, gamma = dg
    expected_value = np.zeros(4)
    ppf = stats.norm.ppf
    pdf = stats.norm.pdf
    q1, q2, q3, q5, q6, q7 = ppf(np.array((1.0, 2.0, 3.0, 5.0, 6.0, 7.0)) / 8)
    expected_value[0] = 3

    expected_value[1] = ((q7 - q5) + (q3 - q1)) / (q6 - q2)

    q_alpha, q_beta = ppf(np.array((alpha / 100.0, beta / 100.0)))
    expected_value[2] = (2 * pdf(q_alpha) / alpha) / (2 * pdf(q_beta) / beta)

    q_delta, q_gamma = ppf(np.array((delta / 100.0, gamma / 100.0)))
    expected_value[3] = (-2.0 * q_delta) / (-2.0 * q_gamma)

    return expected_value 
Example 79
Project: vnpy_crypto   Author: birforce   File: check_moments.py    MIT License 4 votes vote down vote up
def check_cont_basic():
    #results saved in module global variable

    for distname, distargs in distcont[:]:
        #if distname not in distex_0: continue
        distfn = getattr(stats, distname)
##        np.random.seed(765456)
##        sn = 1000
##        rvs = distfn.rvs(size=sn,*arg)
##        sm = rvs.mean()
##        sv = rvs.var()
##        skurt = stats.kurtosis(rvs)
##        sskew = stats.skew(rvs)
        m,v,s,k = distfn.stats(*distargs, **dict(moments='mvsk'))
        st = np.array([m,v,s,k])
        mask = np.isfinite(st)
        if mask.sum() < 4:
            distnonfinite.append(distname)
        print(distname)
        #print 'stats ', m,v,s,k
        expect = distfn.expect
        expect = lambda *args, **kwds : expect_v2(distfn, *args, **kwds)

        special_kwds = specialcases.get(distname, {})
        mnc0 = expect(mom_nc0, args=distargs, **special_kwds)
        mnc1 = expect(args=distargs, **special_kwds)
        mnc2 = expect(mom_nc2, args=distargs, **special_kwds)
        mnc3 = expect(mom_nc3, args=distargs, **special_kwds)
        mnc4 = expect(mom_nc4, args=distargs, **special_kwds)

        mnc1_lc = expect(args=distargs, loc=1, scale=2, **special_kwds)
        #print mnc1, mnc2, mnc3, mnc4
        try:
            me, ve, se, ke = mnc2mvsk((mnc1, mnc2, mnc3, mnc4))
        except:
            print('exception', mnc1, mnc2, mnc3, mnc4, st)
            me, ve, se, ke = [np.nan]*4
            if mask.size > 0:
                distex.append(distname)
        #print 'expect', me, ve, se, ke,
        #print mnc1, mnc2, mnc3, mnc4

        em = np.array([me, ve, se, ke])

        diff = st[mask] - em[mask]
        print(diff, mnc1_lc - (1 + 2*mnc1))
        if np.size(diff)>0 and np.max(np.abs(diff)) > 1e-3:
            distlow.append(distname)
        else:
            distok.append(distname)

        res[distname] = [mnc0, st, em, diff, mnc1_lc] 
Example 80
Project: vnpy_crypto   Author: birforce   File: descriptive.py    MIT License 4 votes vote down vote up
def ci_kurt(self, sig=.05, upper_bound=None, lower_bound=None):
        """
        Returns the confidence interval for kurtosis.

        Parameters
        ----------

        sig : float
            The significance level.  Default is .05

        upper_bound : float
            Maximum value of kurtosis the upper limit can be.
            Default is .99 confidence limit assuming normality.

        lower_bound : float
            Minimum value of kurtosis the lower limit can be.
            Default is .99 confidence limit assuming normality.

        Returns
        --------
        Interval : tuple
            Lower and upper confidence limit

        Notes
        -----
        For small n, upper_bound and lower_bound may have to be
        provided by the user.  Consider using test_kurt to find
        values close to the desired significance level.

        If function returns f(a) and f(b) must have different signs, consider
        expanding the bounds.
        """
        endog = self.endog
        nobs = self.nobs
        if upper_bound is None:
            upper_bound = kurtosis(endog) + \
            (2.5 * (2. * ((6. * nobs * (nobs - 1.)) / \
              ((nobs - 2.) * (nobs + 1.) * \
               (nobs + 3.))) ** .5) * \
               (((nobs ** 2.) - 1.) / ((nobs - 3.) *\
                 (nobs + 5.))) ** .5)
        if lower_bound is None:
            lower_bound = kurtosis(endog) - \
            (2.5 * (2. * ((6. * nobs * (nobs - 1.)) / \
              ((nobs - 2.) * (nobs + 1.) * \
               (nobs + 3.))) ** .5) * \
               (((nobs ** 2.) - 1.) / ((nobs - 3.) *\
                 (nobs + 5.))) ** .5)
        self.r0 = chi2.ppf(1 - sig, 1)
        llim = optimize.brentq(self._ci_limits_kurt, lower_bound, \
                             kurtosis(endog))
        ulim = optimize.brentq(self._ci_limits_kurt, kurtosis(endog), \
                             upper_bound)
        return   llim, ulim