本文參考：https://zhuanlan.zhihu.com/p/29399126（詳細參考該鏈接，本文只是記錄，以備后續查看）
6. 例子 1 對一台設備進行壽命檢驗，記錄10次無故障工作時間，並按照從小到大的次序排列如下：（單位小時）
420 500 920 1380 1510 1650 1760 2100 2300 2350
用Kolmogorov-Smirnov檢驗方法檢驗這些數據的分布是否為參數為1/1500的指數分布？
解：（1） 用ks.test()函數進行檢驗
x<-c(420,500,920,1380 ,1510, 1650,1760,2100,2300, 2350)
ks.test(x,"pexp",1/1500)
運行結果如下：
One-sample Kolmogorov-Smirnov test
data: x
D = 0.30148, p-value = 0.2654
alternative hypothesis: two-sided
因為p值大於0.05，無法拒絕原假設，因此接受這些數據服從為參數為1/1500的指數分布。
（2）用統計模擬方法計算p值.
- 作業：（1）用統計模擬方法計算例子1中的p值。
x<-c(420, 500 ,920 ,1380, 1510, 1650, 1760 ,2100 ,2300 ,2350)
n<-length(x)
d<-max(1:n/n-pexp(x,1/1500),pexp(x,1/1500)-0:(n-1)/n)
r<-10000
D<-0
for(i in 1:r){u<-runif(n);y<-sort(u);D[i]<-max(1:n/n-y,y-0:(n-1)/n)}
p<-length(D[D>d])/r
p
p值的運行結果為：
[1] 0.2656
即模擬的p值為0.2656.
－－－－－－－－－－－－－－－－－－
（2）假設從一個總體中抽取容量為14的樣本：164, 142, 110, 153, 103, 52, 174, 88, 178, 184, 58, 62, 132, 128，問該總體是否服從區間為（50，200）的均勻分布？
題目要求：首先用ks.test()函數進行檢驗，然后用統計模擬方法計算p值，說明檢驗結果。
解：首先用ks.test()函數進行檢驗，
x<-c(164, 142, 110, 153, 103, 52, 174, 88, 178, 184, 58, 62, 132, 128)
ks.test(x,"punif",50,200)
結果如下：
One-sample Kolmogorov-Smirnov test
data: x
D = 0.13429, p-value = 0.9332
alternative hypothesis: two-sided
由於p值= 0.9332>0.05，因此接受該總體服從區間為（50，200）的均勻分布。
－－－－－－－－－－－
統計模擬的方法
x<-c(164, 142, 110, 153, 103, 52, 174, 88, 178, 184, 58, 62, 132, 128)
n<-length(x)
y<-sort(x)
d<-max(1:n/n-punif(y,50,200),punif(y,50,200)-0:(n-1)/n)
r<-10000
D<-0
for(i in 1:r){
u<-runif(n)
y<-sort(u)
D[i]<-max(1:n/n-y,y-0:(n-1)/n)
}
p<-length(D[D>d])/r
p
p值的運行結果為：
[1] 0.9357
即模擬的p值為0.9357.
（3）假設從一個總體中抽取容量為20的樣本：-1.4716958 0.2431494 0.8265337 1.8980804 -4.0423784 3.4327539 5.9290206 1.8254857 -2.4829388 2.5908620 -2.7393827 7.1431654 2.8001338 -0.8907646 6.5693761 -4.2833803 -1.1415025 -8.5100896 2.6901414 10.2162229
該總體是否服從N(1,16)？
x<-c(-1.4716958, 0.2431494, 0.8265337, 1.8980804, -4.0423784, 3.4327539 , 5.9290206, 1.8254857, -2.4829388 , 2.5908620, -2.7393827, 7.1431654, 2.8001338 ,-0.8907646, 6.5693761, -4.2833803, -1.1415025, -8.5100896 , 2.6901414 ,10.2162229)
n<-length(x)
y<-sort(x)
d<-max(1:n/n-pnorm(y,1,4),pnorm(y,1,4)-0:(n-1)/n)
r<-10000
D<-0
for(i in 1:r){
u<-runif(n)
y<-sort(u)
D[i]<-max(1:n/n-y,y-0:(n-1)/n)
}
p<-length(D[D>d])/r
p
KS檢驗是如何工作的？
首先觀察下分析數據
對於以下兩組數據：
controlB={1.26, 0.34, 0.70, 1.75, 50.57, 1.55, 0.08, 0.42, 0.50, 3.20, 0.15, 0.49, 0.95, 0.24, 1.37, 0.17, 6.98, 0.10, 0.94, 0.38}
treatmentB= {2.37, 2.16, 14.82, 1.73, 41.04, 0.23, 1.32, 2.91, 39.41, 0.11, 27.44, 4.51, 0.51, 4.50, 0.18, 14.68, 4.66, 1.30, 2.06, 1.19}
對於controlB，這些數據的統計描述如下：
Mean = 3.61
Median = 0.60
High = 50.6 Low = 0.08
Standard Deviation = 11.2
可以發現這組數據並不符合正態分布， 否則大約有15%的數據會小於均值-標准差（3.61-11.2），而數據中顯然沒有小於0的數。
觀察數據的累計分段函數（Cumulative Fraction Function）
對controlB數據從小到大進行排序：
sorted controlB={0.08, 0.10, 0.15, 0.17, 0.24, 0.34, 0.38, 0.42, 0.49, 0.50, 0.70, 0.94, 0.95, 1.26, 1.37, 1.55, 1.75, 3.20, 6.98, 50.57}。10%的數據（2/20）小於0.15，85%（17/20）的數據小於3。所以，對任何數x來說，其累計分段就是所有比x小的數在數據集中所占的比例。下圖就是controlB數據集的累計分段圖
可以看到大多數數據都幾種在圖片左側（數據值比較小），這就是非正態分布的標志。為了更好的觀測數據在x軸上的分布，可以對x軸的坐標進行非等分的划分。在數據都為正的時候有一個很好的方法就是對x軸進行log轉換。下圖就是上圖做log轉換以后的圖：
將treatmentB的數據也做相同的圖（如下），可以發現treatmentB和controlB的數據分布范圍大致相同（0.1 - 50）。但是對於大部分x值，在controlB數據集中比x小的數據所占的比例比在treatmentB中要高，也就是說達到相同累計比例的值在treatment組中比control中要高。KS檢驗使用的是兩條累計分布曲線之間的最大垂直差作為D值（statistic D）作為描述兩組數據之間的差異。在此圖中這個D值出現在x=1附近，而D值為0.45（0.65-0.25）。
值得注意的是雖然累計分布曲線的性狀會隨着對數據做轉換處理而改變（如log轉換），但是D值的大小是不會變的。