正态性检验
x2yline
2017年10月18日
正态性检验最好不使用直接计算法检验偏度和峰度
** SPSS默认的两个方法为Shapiro – Wilk（W 检验）和Kolmogorov – Smirnov（D 检验）**
参考文章：
https://www.yaolibio.com/2017/01/03/normality-test/
http://blog.sina.com.cn/s/blog_5efe47110100d28m.html
http://blog.sina.com.cn/s/blog_403aa80a01019ly5.html
但是其中所说的正态检验水准一般设置为0.10
当样本含量n ≤2000时，结果以Shapiro – Wilk（W 检验）为准
# faithful是R里的一个内置数据集
shapiro.test(faithful$eruptions)
##
## Shapiro-Wilk normality test
##
## data: faithful$eruptions
## W = 0.84592, p-value = 9.036e-16
当样本含量n >2000 时，结果以Kolmogorov – Smirnov（D 检验）为准
ks.test(faithful$eruptions, "pnorm")
## Warning in ks.test(faithful$eruptions, "pnorm"): ties should not be present
## for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: faithful$eruptions
## D = 0.94857, p-value < 2.2e-16
## alternative hypothesis: two-sided
Warning message:
In ks.test(faithful$eruptions, "pnorm") :
ties should not be present for the Kolmogorov-Smirnov test
出现警告，Kolmogorov - Smirnov检验里不应该有相等的值，可以忽略或加上参数exact=FALSE，默认使用近似算法计算近似p值。
直接计算法检验（纯练习用，可忽略）
1. 用R语言计算faithful数据框中的eruptions的3阶中心距(center moment)
e1071::moment(faithful$eruptions,
order=3, center=TRUE)
## [1] -0.6149059
2. 用R语言计算偏度，偏度计算公式为
μ2 和 μ3 为二阶和三阶中心距，用moment计算中心距时，其除数为n而不是n-1，这里需要注意一下
R语言验证:
u3 <- e1071::moment(faithful$eruptions, order=3, center=TRUE)
u2 <- e1071::moment(faithful$eruptions, order=2, center=TRUE)
u3 <- mean((faithful$eruptions - mean(faithful$eruptions))^3)
u3/(u2)^(3/2)
## [1] -0.415841
u3/(u2*nrow(faithful)/(nrow(faithful)-1))^1.5
## [1] -0.4135498
e1071::skewness(faithful$eruptions)
## [1] -0.4135498
agricolae::skewness(faithful$eruptions)
## [1] -0.4181505
moments::skewness(faithful$eruptions)
## [1] -0.415841
propagate::skewness(faithful$eruptions)
## [1] -0.415841
skew <- agricolae::skewness(faithful$eruptions)
propagate、moments、e1071包和agricolae的包计算结果的差异是由于e1071算的三阶中心距是有偏的估计，而propagate和moments的三阶和二阶距估计都是有偏的，所以agricolae的结果更为准确
偏度小于0说明均值小于中位数，呈左偏分布（负偏分布）
PS: 教材的第39页公式似乎更加精确
3. 用R语言计算峰度
u4 <- e1071::moment(faithful$eruptions, order=4, center=TRUE)
u2 <- e1071::moment(faithful$eruptions, order=2, center=TRUE)
u4/u2^2 -3
## [1] -1.5006
u4/var(faithful$eruptions)^2 - 3
## [1] -1.511605
e1071::kurtosis(faithful$eruptions)
## [1] -1.511605
agricolae::kurtosis(faithful$eruptions)
## [1] -1.506167
moments::kurtosis(faithful$eruptions)
## [1] 1.4994
moments::kurtosis(faithful$eruptions)-3
## [1] -1.5006
propagate::kurtosis(faithful$eruptions)
## [1] -1.5006
kurt <- agricolae::kurtosis(faithful$eruptions)
propagate、moments、e1071包和agricolae的包计算结果的差异是由于e1071算的四中心距是有偏的估计，而propagate和moments的都是有偏估计，所以agricolae的结果更为准确
峰度等于0说明标准化后的分布与正态分布的峰同样尖锐，大于0说明更加尖锐，尾部更粗；小于0说明更加平坦，尾部更细。
PS: 教材的第39页公式似乎更加精确
4. 偏度和峰度的假设检验
H0：skew = 0, krut = 0，即总体服从正态分布
H1：skew ！= 0或（和）krut！= 0，即总体不服从正态分布
计算出峰度和偏度以后，还能计算出峰度和偏度的标准误
- 峰度和偏度标准误
ses <- sqrt(6 / length(faithful$eruptions))
sek <- sqrt(24/length(faithful$eruptions))
- 检验统计量
totests <- skew/ses
totestk <- kurt/sek
- 计算p值
pvals <- pt(totests, (length(faithful$eruptions) - 1))
pvalk <- pt(totestk, (length(faithful$eruptions) - 1))
pvals
## [1] 0.002614556
pvalk
## [1] 3.68164e-07
所以按0.01水平，拒绝H0，接受H1，这些样本不服从正态分布
其他正态性检验R包nortest和normtest
R自带的stats包里有Kolmogorov-Smirnov检验和Shapiro-Wilk检验，同时nortest包里的Lilliefor检验、Anderson-Darling检验···和tseries包中的Jarque-Bera检验也能进行检验
# Anderson-Darling test for normality
nortest::ad.test(faithful$eruptions)
##
## Anderson-Darling normality test
##
## data: faithful$eruptions
## A = 17.305, p-value < 2.2e-16
# Lilliefors (Kolmogorov-Smirnov) test for normality
# Although the test statistic obtained from lillie.test(x) is the same as that obtained from
# ks.test(x, "pnorm", mean(x), sd(x)), it is not correct to use the p-value from the latter for
# the composite hypothesis of normality (mean and variance unknown), since the distribution of the
# test statistic is different when the parameters are estimated.
nortest::lillie.test(faithful$eruptions)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: faithful$eruptions
## D = 0.18135, p-value < 2.2e-16
# jarque.bera.test: This test is a joint statistic using skewness and kurtosis coefficients.
tseries::jarque.bera.test(faithful$eruptions)
##
## Jarque Bera Test
##
## data: faithful$eruptions
## X-squared = 33.36, df = 2, p-value = 5.702e-08
normtest::skewness.norm.test(faithful$eruptions)
##
## Skewness test for normality
##
## data: faithful$eruptions
## T = -0.41584, p-value = 0.006
normtest::kurtosis.norm.test(faithful$eruptions)
##
## Kurtosis test for normality
##
## data: faithful$eruptions
## T = 1.4994, p-value = 0.001
https://cran.r-project.org/web/packages/nortest/nortest.pdf
https://cran.r-project.org/web/packages/normtest/normtest.pdf
https://stats.stackexchange.com/questions/3136/how-to-perform-a-test-using-r-to-see-if-data-follows-normal-distribution
https://cran.r-project.org/web/packages/tseries/tseries.pdf
参考：
http://www.r-tutor.com/elementary-statistics/numerical-measures/variance
http://www.r-tutor.com/elementary-statistics/numerical-measures/moment
https://stat.ethz.ch/pipermail/r-help/1999-July/004529.html
http://www.r-tutor.com/elementary-statistics/numerical-measures/skewness
http://www.r-tutor.com/elementary-statistics/numerical-measures/kurtosis
https://artax.karlin.mff.cuni.cz/r-help/library/agricolae/html/skewness.html
https://artax.karlin.mff.cuni.cz/r-help/library/agricolae/html/kurtosis.html
https://www.researchgate.net/post/Could_anyone_tell_me_how_to_calculate_skewness_and_kurtosis_of_a_numeric_variable_in_R
