偏度和峰度的置信区间及标准误差

Question

请告诉我如何计算偏度和峰度，以及它们各自的标准误差和置信区间(即偏度SE和峰度S.E )我找到了两个软件包

1) package:'measure'只能计算偏度和峰度。

2) package:'rela'可以同时计算偏度和峰度，但在计算过程中，默认情况下使用引导，没有命令将其关闭。

Answer 1

我只是在复制和粘贴由Howard在这里中发布的代码

# Skewness and kurtosis and their standard errors as implement by SPSS
#
# Reference: pp 451-452 of
# http://support.spss.com/ProductsExt/SPSS/Documentation/Manuals/16.0/SPSS 16.0 Algorithms.pdf
# 
# See also: Suggestion for Using Powerful and Informative Tests of Normality,
# Ralph B. D'Agostino, Albert Belanger, Ralph B. D'Agostino, Jr.,
# The American Statistician, Vol. 44, No. 4 (Nov., 1990), pp. 316-321

spssSkewKurtosis=function(x) {
  w=length(x)
  m1=mean(x)
  m2=sum((x-m1)^2)
  m3=sum((x-m1)^3)
  m4=sum((x-m1)^4)
  s1=sd(x)
  skew=w*m3/(w-1)/(w-2)/s1^3
  sdskew=sqrt( 6*w*(w-1) / ((w-2)*(w+1)*(w+3)) )
  kurtosis=(w*(w+1)*m4 - 3*m2^2*(w-1)) / ((w-1)*(w-2)*(w-3)*s1^4)
  sdkurtosis=sqrt( 4*(w^2-1) * sdskew^2 / ((w-3)*(w+5)) )
  mat=matrix(c(skew,kurtosis, sdskew,sdkurtosis), 2,
        dimnames=list(c("skew","kurtosis"), c("estimate","se")))
  return(mat)
}

要获得变量的偏度和峰度以及它们的标准错误，只需运行以下函数：

x <- rnorm(100)
spssSkewKurtosis(x)

##             estimate    se
##    skew       -0.684 0.241
##    kurtosis    0.273 0.478

Answer 2

标准错误对正态分布有效，但对其他分布无效。为了了解原因，您可以运行以下代码(它使用了上面所示的spssSkewKurtosis函数)来估计通过取峰度估计值正负1.96标准误差而得到的区间的真实置信水平：

set.seed(12345)
Nsim = 10000
Correct = numeric(Nsim)
b1.ols = numeric(Nsim)
b1.alt = numeric(Nsim)
for (i in 1:Nsim) {
 Data = rnorm(1000)  
 Kurt = spssSkewKurtosis(Data)[2,1]
 seKurt =  spssSkewKurtosis(Data)[2,2]
  LowerLimit = Kurt -1.96*seKurt
  UpperLimit = Kurt +1.96*seKurt
  Correct[i] = LowerLimit <= 0 & 0 <= UpperLimit  
 }

TrueConfLevel = mean(Correct)
TrueConfLevel

这给你0.9496，可以接受地接近预期的95%，所以当数据来自正态分布时，标准错误会像预期的那样工作。但是，如果将Data = rnorm(1000)改为Data = runif(1000)，则假设数据来自一个均匀分布，其理论(超额)峰度为-1.2。从Correct[i] = LowerLimit <= 0 & 0 <= UpperLimit到Correct[i] = LowerLimit <= -1.2 & -1.2 <= UpperLimit进行相应的更改会得到结果1.0，这意味着95%的间隔总是正确的，而不是95%的样本。因此，对于(轻尾)均匀分布，标准误差似乎被高估(太大).

如果将Data = rnorm(1000)改为Data = rexp(1000)，则假设数据来自指数分布，其理论(超额)峰度为6.0。从Correct[i] = LowerLimit <= 0 & 0 <= UpperLimit到Correct[i] = LowerLimit <= 6.0 & 6.0 <= UpperLimit的相应变化给出了结果0.1007，这意味着95%的区间仅对10.07%的样本是正确的，而对95%的样本则是正确的。因此，对于(重尾)指数分布，标准误差似乎被低估(太小).

这些标准错误对于非正态分布是非常不正确的，正如上面的模拟所示.因此，这些标准误差的唯一用途是将估计的峰度与预期的理论正常值(0.0)进行比较；例如，使用假设检验。它们不能用来构造真实峰度的置信区间。

Answer 3

@HBat 是对的：如果您的样本数据是高斯的，您可以使用来自维基百科的方程计算标准误差

n = len(sample)
se_skew = ((6*n*(n-1))/((n-2)*(n+1)*(n+3)))**0.5

但是，@BigBendRegion 也是对的：如果您的数据不是高斯的，这是行不通的。那么你可能需要引导。

R有一个DescTools包，它可以引导信任区间来进行倾斜(除其他事情外)。可以使用rpy2将其包括在python中，如下所示：

""" Import rpy2 and the relevant package"""
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
DescTools = importr('DescTools')
""" You will probably need this if you want to work with numpy arrays"""
import rpy2.robjects.numpy2ri
rpy2.robjects.numpy2ri.activate()


def compute_skew(data, confidence_level=0.99):
    """ Compute the skew and confidence interval using rpy2, DescTools
        @param data
        @return dict with keys: skew, skew_ci_lower, skew_ci_upper"""
    d = {}
    d["skew"], d["skew_ci_lower"], d["skew_ci_upper"] = DescTools.Skew(data, conf_level=confidence_level)
    return d

""" Call the function on your data (assuming that is saved in a variable named sample)"""
print(compute_skew(sample))