3

请告诉我如何计算偏度和峰度以及它们各自的标准误差和与之相关的置信区间（即偏度的 SE 和峰度的 SE）我找到了两个包

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

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

4

5 回答 5

4

我只是在这里复制并粘贴 Howard Seltman 发布的代码：

# 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
于 2019-01-25T16:49:42.683 回答
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) 进行比较；例如，使用假设检验。它们不能用于构建真实峰度的置信区间。

于 2019-06-13T12:19:08.240 回答
1

@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))
于 2020-08-26T15:18:45.133 回答
0 
def skew_kurt(dataframe: pd.DataFrame) -> pd.DataFrame:

out = []
for col in dataframe:
    
    x = dataframe[col]
    sd = x.std()
    
    if sd == 0:
        out.append({name: np.nan for name in ['skew stat', 'skew se', 'kurt stat', 'kurt se']})
        continue
    
    w, m1 = len(x), x.mean()
    dif = x - m1
    m2, m3, m4 = tuple([(dif ** i).sum() for i in range(2, 5)])
    
    skew = w * m3 / (w - 1) / (w - 2) / sd ** 3
    skew_se = np.sqrt(6 * w * (w - 1) / ((w - 2) * (w + 1) * (w + 3)))
    
    kurt = (w * (w + 1) * m4 - 3 * m2 ** 2 * (w - 1)) / ((w - 1) * (w - 2) * (w - 3) * sd ** 4)
    kurt_se = np.sqrt(4 * (w ** 2 - 1) * skew_se ** 2 / ((w - 3) * (w + 5)))

    out.append({'skew stat': skew, 'skew se': skew_se, 'kurt stat': kurt, 'kurt se': kurt_se})
    
dataframe = pd.DataFrame(out, index = list(dataframe))
dataframe['skew<2'] = np.absolute(dataframe['skew stat']) < 2 * dataframe['skew se']
dataframe['kurt<2'] = np.absolute(dataframe['kurt stat']) < 2 * dataframe['kurt se']
    
return dataframe[['skew stat', 'skew se', 'skew<2', 'kurt stat', 'kurt se', 'kurt<2']]
于 2022-02-10T05:02:20.340 回答
0

尝试包装心理：

> a <- data.frame(cola=rep(c('A','B','C'),100),colb=sample(1:1000,300),colc=rnorm(300))
> describe(a)
      vars   n   mean     sd median trimmed    mad   min    max  range  skew kurtosis    se
cola*    1 300   2.00   0.82   2.00    2.00   1.48  1.00   3.00   2.00  0.00    -1.51  0.05
colb     2 300 511.76 285.59 506.50  514.21 362.50  1.00 999.00 998.00 -0.04    -1.17 16.49
colc     3 300   0.12   1.04   0.05    0.10   1.07 -2.54   2.91   5.45  0.12    -0.24  0.06
> describe(a)$skew
[1]  0.00000000 -0.04418551  0.11857609
于 2016-07-12T12:05:51.833 回答