R：用于正态性检验的 Kolmogorov-Smirnov 检验可视化经验 CDF 与理论 CDF 的比较

问题描述

我想使用 Kolmogorov-Smirnov 检验 (ks-test) 检查正态性并将其可视化以支持我的演示。

我使用了这个默认代码

sample1 <- rnorm(10000,10,5)
sample2 <- rnorm(10000,1,5)
group <- c(rep("sample1",length(sample1)),rep("sample2",length(sample2)))
dat <- data.frame(KSD = c(sample1,sample2),group = group)
# create ECDF of data
cdf1 <- ecdf(sample1) 
cdf2 <- ecdf(sample2) 
# find min and max statistics to draw line between points of greatest distance
minMax <- seq(min(sample1,max(sample1,length.out=length(sample1)) 
x0 <- minMax[which( abs(cdf1(minMax) - cdf2(minMax)) == max(abs(cdf1(minMax) - cdf2(minMax))) )] 
y0 <- cdf1(x0) 
y1 <- cdf2(x0) 

ggplot(dat,aes(x = KSD,group = group,color = group))+
  stat_ecdf(size=1) +
  theme_bw(base_size = 28) +
  theme(legend.position ="top") +
  xlab("Sample") +
  ylab("ECDF") +
  #geom_line(size=1) +
  geom_segment(aes(x = x0[1],y = y0[1],xend = x0[1],yend = y1[1]),linetype = "dashed",color = "red") +
  geom_point(aes(x = x0[1],y= y0[1]),color="red",size=8) +
  geom_point(aes(x = x0[1],y= y1[1]),size=8) +
  ggtitle("K-S Test: Sample 1 / Sample 2") +
  theme(legend.title=element_blank())

结果是这样的：

现在，我想用我自己的数据替换它

mydata <- c(33,28,23,27,22,25,32,26,14,37,16,30,34,18,31,20,24,12,40,17,29,36,46,33,21,35,13,19,39,11,15,42,38,50,44,41,27)

我可以将 mydata 用作经验 CDF，但如何将 sample2 更改为理论 CDF？

目标：

• 2 cdf 图：经验 CDF（mydata 替代 sample1）和理论 CDF（替代 sample2 分布
• 出现最大差异的垂直线 + geom_pint 用于加法

提前致谢。

解决方法

R 中的理论累积概率函数由分布函数的 p* 版本给出。如果你想测试正态性，函数应该是 pnorm。在下面，您将找到一种绘制此图的方法。黑线表示最大差异点。

library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.0.5

mydata <- c(33,28,23,10,27,22,25,32,26,14,37,16,30,34,18,31,20,24,12,40,17,29,36,46,33,21,35,13,19,39,11,15,42,38,50,44,41,27)

# Find x position of maximum difference.
# You could precalculate `mean(mydata)` and `sd(mydata)` if you wish,# we use it later in plotting too.
maxdiff <- mydata[which.max(
  abs(ecdf(mydata)(mydata) - pnorm(mydata,mean(mydata),sd = sd(mydata)))
)]

df <- data.frame(x = mydata)


ggplot(df,aes(x)) +
  stat_ecdf(aes(colour = "Empirical")) +
  stat_function(fun = pnorm,args = list(mean = mean(mydata),sd = sd(mydata)),aes(colour = "Theoretical")) +
  annotate(
    "segment",x = maxdiff,xend = maxdiff,y = ecdf(mydata)(maxdiff),yend = pnorm(maxdiff,sd(mydata)),size = 1
  )

^{由 reprex package (v1.0.0) 于 2021 年 4 月 25 日创建}