logo标签列表

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

rstandard-errorkurtosis
3
请告诉我如何计算偏度和峰度，以及与之相关的标准误差和置信区间（即偏度的SE和峰度的SE）。我找到了两个包：
1）包：'measure' 只能计算偏度和峰度
2）包：'rela' 可以计算偏度和峰度，但默认使用bootstrap，没有命令在计算过程中关闭它。
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
2
请注意，使用的SE公式假定x服从正态分布。 - Remi Cuingnet
3

标准误差仅适用于正态分布，而不适用于其他分布。要了解原因，您可以运行以下代码（使用上面显示的spssSkewKurtosis函数），以估计通过取kurtosis估计值加减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）进行比较，例如使用假设检验。它们不能用于构建真实峰度的置信区间。

1

@HBat 正确：如果您的样本数据符合高斯分布，您可以使用wikipedia中的方程式计算标准误差。

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))
0
尝试使用 psych 包：
> 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
1
我认为OP想要偏度和峰度本身的标准误差。 - Ben Bolker
2
SES和SEK可以直接从样本大小n计算得出：其中SES=分数的平方根，6n乘以n-1在上面，n-2乘以n+1，乘以n+3在下面。 - Zahiro Mor
请查看以下链接：https://estatistics.eu/what-is-statistics-standard-error-of-skewness-standard-error-of-kurtosis/ - Zahiro Mor
2
你找到的是一个粗略估计值（sqrt(6/n) - 我建议使用：ses <- function(n) {sqrt((6*n*(n-1))/((n-2)*(n+1)*(n+3)))} - Zahiro Mor
@ZahiroMor，这个公式也适用于非正态样本吗？当数据来自正态分布时，可以使用此公式，但是当分布是偏斜的时，我不确定是否可以使用该公式来计算偏度估计的标准误差。 - Skumin
显示剩余2条评论
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']]
1
你的回答可以通过添加进一步的支持信息来改善，请[编辑]以添加更多细节，例如引用或文档，以便其他人可以确认您的答案是正确的。您可以在帮助中心中找到有关如何撰写良好答案的更多信息。 - Community
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接