# 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
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
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
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
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)
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
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
def descriptive_stats(self):
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
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
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
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
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
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
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
start_params[-1] = np.sqrt(res_ols.scale)

self.start_params = start_params 
Example 13
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
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
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
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
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
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]),
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
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
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
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
def test_kurtosis_array_scalar(self):
assert_equal(type(stats.kurtosis([1,2,3])), float) 
Example 23
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
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
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
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
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
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
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
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
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)
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
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
def test_rolling_kurt(self):
from scipy.stats import kurtosis
self._check_moment_func(lambda x: kurtosis(x, bias=False),
name='kurt') 
Example 34
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
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
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
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
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]),
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
def test_rolling_kurt(self):
from scipy.stats import kurtosis
self._check_moment_func(lambda x: kurtosis(x, bias=False),
name='kurt') 
Example 54
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
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
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
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
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
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
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
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
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
def test_rolling_kurt(self):
from scipy.stats import kurtosis
self._check_moment_func(lambda x: kurtosis(x, bias=False),
name='kurt') 
Example 64
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
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
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]),
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
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
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
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
def test_kurtosis_array_scalar(self):
assert_equal(type(stats.kurtosis([1,2,3])), float) 
Example 71
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
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
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
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
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]),
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
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
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
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
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])
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
distex.append(distname)
#print 'expect', me, ve, se, ke,
#print mnc1, mnc2, mnc3, mnc4

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

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