【问题标题】:How to correctly use scipy's skew and kurtosis functions?如何正确使用 scipy 的 skew 和 kurtosis 函数?
【发布时间】:2018-01-11 00:57:27
【问题描述】:

skewness 是衡量数据集对称性的参数,kurtosis 是衡量其尾部与正态分布相比的重度,例如参见 @ 987654321@.

scipy.stats 提供了一种计算这两个量的简单方法,请参阅scipy.stats.kurtosisscipy.stats.skew

据我了解,使用刚才提到的函数,normal distribution 的偏度和峰度都应该为 0。但是,我的代码并非如此:

import numpy as np
from scipy.stats import kurtosis
from scipy.stats import skew

x = np.linspace( -5, 5, 1000 )
y = 1./(np.sqrt(2.*np.pi)) * np.exp( -.5*(x)**2  )  # normal distribution

print( 'excess kurtosis of normal distribution (should be 0): {}'.format( kurtosis(y) ))
print( 'skewness of normal distribution (should be 0): {}'.format( skew(y) ))

输出是:

正态分布的超峰度(应为0):-0.307393087742

正态分布偏度(应为0):1.11082371392

我做错了什么?

我使用的版本是

python: 2.7.6
scipy : 0.17.1
numpy : 1.12.1

【问题讨论】:

    标签: python numpy scipy statistics


    【解决方案1】:

    这些函数计算probability density distribution 的矩(这就是为什么它只需要一个参数)并且不关心值的“函数形式”。

    这些适用于“随机数据集”(将它们视为均值、标准差、方差等度量):

    import numpy as np
    from scipy.stats import kurtosis, skew
    
    x = np.random.normal(0, 2, 10000)   # create random values based on a normal distribution
    
    print( 'excess kurtosis of normal distribution (should be 0): {}'.format( kurtosis(x) ))
    print( 'skewness of normal distribution (should be 0): {}'.format( skew(x) ))
    

    给出:

    excess kurtosis of normal distribution (should be 0): -0.024291887786943356
    skewness of normal distribution (should be 0): 0.009666157036010928
    

    改变随机值的数量可以提高准确度:

    x = np.random.normal(0, 2, 10000000)
    

    导致:

    excess kurtosis of normal distribution (should be 0): -0.00010309478605163847
    skewness of normal distribution (should be 0): -0.0006751744848755031
    

    在您的情况下,函数“假定”每个值具有相同的“概率”(因为这些值是均匀分布的,并且每个值只出现一次)所以从 skewkurtosis 的角度来看,它正在处理具有非高斯概率密度(不确定这到底是什么),这解释了为什么结果值甚至不接近 0

    import numpy as np
    from scipy.stats import kurtosis, skew
    
    x_random = np.random.normal(0, 2, 10000)
    
    x = np.linspace( -5, 5, 10000 )
    y = 1./(np.sqrt(2.*np.pi)) * np.exp( -.5*(x)**2  )  # normal distribution
    
    import matplotlib.pyplot as plt
    
    f, (ax1, ax2) = plt.subplots(1, 2)
    ax1.hist(x_random, bins='auto')
    ax1.set_title('probability density (random)')
    ax2.hist(y, bins='auto')
    ax2.set_title('(your dataset)')
    plt.tight_layout()
    

    【讨论】:

    • 很好的解释,谢谢!我在这里“误用”了峰度和偏度函数。比较显示两个直方图就很清楚了。 (我最初的动机是以某种方式量化曲线的重尾;看来我必须考虑一些不同的东西)。
    • 如果要可视化尾部,请使用正常的分位数 - 分位数 (q-q) 图而不是直方图。即使尾部很重,尾部密度也接近于零,因此在直方图中不容易看到尾部。但在 q-q 图中,重尾很明显。 q-q 图的视觉外观与实际峰度统计量之间也存在简单的数学联系。
    【解决方案2】:

    您正在使用密度函数的“形状”作为数据。 这些函数旨在与从分布中采样的数据一起使用。 如果您从分布中抽样,您将获得随着样本量增加而接近正确值的样本统计量。 要绘制数据,我建议使用直方图。

    %matplotlib inline
    import numpy as np
    import pandas as pd
    from scipy.stats import kurtosis
    from scipy.stats import skew
    
    import matplotlib.pyplot as plt
    
    plt.style.use('ggplot')
    
    data = np.random.normal(0, 1, 10000000)
    np.var(data)
    
    plt.hist(data, bins=60)
    
    print("mean : ", np.mean(data))
    print("var  : ", np.var(data))
    print("skew : ",skew(data))
    print("kurt : ",kurtosis(data))
    

    输出:

    mean :  0.000410213500847
    var  :  0.999827716979
    skew :  0.00012294118186476907
    kurt :  0.0033554829466604374
    

    除非您处理的是解析表达式,否则在使用数据时获得零的可能性极小。

    【讨论】:

    • 直方图的要点。但我会使用更多的垃圾箱,我个人最喜欢plt.hits(data, bins='auto') - 但这需要更新的 NumPy 版本,因此可能不适合所有读者。
    • 谢谢,我现在看到我“误用”了函数峰度和偏度(我想量化曲线/信号的“重尾”)。
    • 也许您应该打开另一个问题,提供有关信号的更多信息
    • @purpleTentacle 我正在考虑这一点,但首先我试图在文献中找到一些东西(而且,它可能不属于这里,但更多的是物理学或统计学,我猜)
    猜你喜欢
    • 2016-01-11
    • 2017-05-16
    • 2020-08-09
    • 2019-12-31
    • 2023-04-02
    • 1970-01-01
    • 2019-02-04
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多