本節書摘來自華章計算機《數學模組化：基于r》一書中的第1章，第1.4節,作者：薛 毅 更多章節内容可以通路雲栖社群“華章計算機”公衆号檢視。
本節要介紹的是另一類檢驗,其目标不是針對具體的參數,而是針對分布的類型.例如,通常假定總體分布具有正态性,而“總體分布為正态”這一斷言本身在一定場合下就是可疑的,有待于檢驗.
假設根據某理論、學說甚至假定,某随機變量應當有分布f,現對x進行n次觀察,得到一個樣本x1,x2,…,xn,要據此檢驗h0:x具有理論分布f這裡雖然沒有明确指出對立假設,但可以說,對立假設為h1:x不具有理論分布f本問題的真實含義是估計實測資料與該理論或學說符合得怎麼樣,而不在于當認為不符合時,x可能備擇的分布如何.故問題中不明确标出對立假設,反而使人感到提法更為貼近現實.
1.4.1 pearson拟合優度χ2檢驗
1.理論分布完全已知的情況
設x1,x2,…,xn為來自總體x的樣本,将數軸(-∞,∞)分成m個區間i1=(-∞,a1), i2=［a1,a2), …, im=［am-1,∞)記這些區間的理論機率分别為p1, p2, …, pm, pi=p{x∈ii}, i=1,2,…,m記ni為x1,x2,…,xn中落在區間ii内的個數,則在原假設成立的條件下,ni的期望值為npi,ni與npi的差距(i=1,2,…,m)可視為理論與觀察之間偏離的衡量,構造統計量k=∑mi=1(ni-npi)2npi（1.62）稱k為pearson χ2統計量.pearson證明了,在原假設成立的條件下,當n→∞時,k依分布收斂于自由度為m-1的χ2分布.
在r中,使用chisq.test()函數計算pearson拟合優度χ2檢驗,其使用格式為chisq.test(x, y = null, correct = true,
p = rep(1/length(x), length(x)), rescale.p = false,
simulate.p.value = false, b = 2000)參數x為數值向量或矩陣(用于列聯表檢驗),或者x和y同時為因子.y為數值向量,當x為矩陣時,y無效.如果x為因子,y必須為同樣長度的因子.
correct為邏輯變量,表示是否對2×2列聯表中的統計量作連續型修正,預設值為true.
p為向量,表示落在每個小區間内的理論機率,預設值為均勻分布.
rescale.p為邏輯變量,表示是否将輸入的p标準化,預設值為false.當p的和不是1時,需要取rescale.p = true,此時重新計算p,使其和為1.
simulate.p.value為邏輯變量,表示是否使用monte carlo方法仿真計算p值,預設值為false.
b為正整數,表示monte carlo仿真的次數,預設值為2000.
例1.21 某消費者協會為了确定市場上消費者對5種品牌啤酒的喜好情況,随機抽取了1000名啤酒愛好者作為樣本進行如下試驗:5種啤酒按分别标注為a、b、c、d和e，讓啤酒愛好者品嘗，每人分别品嘗過這5種啤酒後，選出一種作為他最喜歡的品牌.表1.8是根據樣本資料整理得到的各種品牌啤酒愛好者的頻數分布.試根據這些資料判斷消費者對這5種品牌啤酒的愛好有無明顯差異.
解 如果消費者對5種品牌啤酒喜好無顯著差異,那麼,就可以認為喜好這5種品牌啤酒的人呈均勻分布,即5種品牌啤酒愛好者人數各占20％.是以原假設為h0:喜好5種啤酒的人數服從均勻分布輸入資料,調用chisq.test()函數計算.> x <- c(210, 312, 170, 85, 223)
chisq.test(x) chi-squared test for given probabilities
data: x
x-squared = 136.49, df = 4, p-value < 2.2e-16x-squared為χ2統計量,df為自由度,p-value為p值(=2.2×10-16)0.05,拒絕原假設,故認為消費者對5種品牌啤酒的喜好是有顯著差異的.
例1.22 大麥的雜交後代關于芒性的比例應是無芒∶長芒∶短芒=9∶3∶4.實際觀測值為335∶125∶160.試檢驗觀測值是否符合理論假設.
解 根據題意h0:p1=916, p2=316, p3=416調用chisq.test()函數,其指令如下:> chisq.test(c(335, 125, 160), p = c(9, 3, 4)/16)
data: c(335, 125, 160)
x-squared = 1.362, df = 2, p-value = 0.5061p值(=0.5061)>0.05,接受原假設,即大麥的芒性比例符合9∶3∶4的比例.
例1.23 為研究電話總機在某段時間内接到的呼叫次數是否服從poisson分布,現收集了42個資料,如表1.9所示.通過對資料的分析,能否确認在某段時間内接到的呼叫次數服從poisson分布(α=0.1)？
解 編寫相應的計算程式(程式名:exam0123.r)如下：x <- 0:6; y <- c(7, 10, 12, 8, 3, 2, 0)
f <- ppois(x, mean(rep(x, y))); m <- length(y)
p <- f［1］; p［m］ <- 1 - f［m-1］
for (i in 2:(m-1)) p［i］ <- f［i］ - f［i-1］
chisq.test(y, p=p)但計算結果會出現警告. chi-squared test for given probabilities
data: y
x-squared = 1.5057, df = 6, p-value = 0.9591
warning message:
chi-squared近似算法有可能不準in: chisq.test(y, p = p)為什麼會出現這種情況呢?這是因為pearson拟合優度χ2檢驗要求:在分組後,每個組的理論頻數和實際頻數都要大于等于5,而後三組中出現的頻數分别為3,2,0,均小于5.解決問題的方法是将後三組合成一組,此時的頻數為5,滿足要求.下面給出相應的r程式.z <- c(7, 10, 12, 8, 5) #% 重新分組
m <- length(z); p <- p［1:m-1］; p［m］ <- 1 - f［m-1］
chisq.test(z, p=p) #% 作檢驗計算得到p值(= 0.9696)0.1,是以,能确認在某段時間内接到的呼叫次數服從poisson分布. chi-squared test for given probabilities
data: z
x-squared = 0.5389, df = 4, p-value = 0.96962.理論分布依賴于若幹個未知參數的情況
如果分布族f依賴于r個參數θ1,θ2,…,θr,要根據樣本x1,x2,…,xn去檢驗假設h0:x的分布屬于分布族{f(x;θ1,θ2,…,θr)}解決這個問題的步驟是,先通過樣本作出(θ1,θ2,…,θr)的極大似然估計(θ1,θ2,…,θr)再檢驗假設h0:x具有分布f(x;θ1,θ2,…,θr)然後再按理論分布已知的情況進行處理,所不同的是由式(1.62)得到的統計量k服從自由度為m-1-r的χ2分布,即自由度減少了r.
從這種角度來講,例1.23的計算是有問題的,因為在計算過程中,是用樣本均值代替總體的數學期望,是以,需要減少一個自由度.
為了便于此類問題的計算,這裡編寫一個小程式(程式名:chi2gof.r)來完成這項工作.chi2gof <- function(chi2test, nparams){
df <- chi2testparameter - nparams;
pval <- 1 - pchisq(chi2teststatistic, df)
names(pval) <- "p-val"
c(chi2teststatistic, pval, df)
}在程式中, chi2test為chisq.test()函數生成的對象,nparams為未知參數的個數.函數的傳回值有三個,分别是χ2統計量、p值和自由度.
用chi2gof()函數再計算例1.23,程式和計算結果如下：> source("chi2gof.r")
chi2gof(chisq.test(z, p = p), nparams = 1)
x-squared p-val df
0.5388945 0.9102670 3.0000000此時,自由度為3而不是4, 雖然p值有所減少但結論不變.
1.4.2 kolmogorov-smirnov檢驗
kolmogorov(科爾莫戈羅夫)拟合優度檢驗是kolmogorov于1933年提出的,它是針對單個總體的檢驗,檢驗樣本是否來自指定分布f0.h0:f(x)=f0(x), h1:f(x)≠f0(x) (雙側檢驗)（1.63）
h0:f(x)≥f0(x), h1:f(x)h0:f(x)≤f0(x), h1:f(x)>f0(x) (單側檢驗)（1.65）smirnov(斯米爾諾夫)拟合優度檢驗是smirnov于1939年提出的,它是針對兩個總體的檢驗,檢驗兩個總體的分布是否相同.h0:f(x)=g(x), h1:f(x)≠g(x) (雙側檢驗)（1.66）
h0:f(x)≥g(x), h1:f(x)h0:f(x)≤g(x), h1:f(x)>g(x) (單側檢驗)（1.68）可以将smirnov檢驗看成kolmogorov檢驗的兩樣本情形,是以,将兩個檢驗統稱為kolmogorov-smirnov檢驗,簡稱k-s檢驗.
在r中,用ks.test()函數作k-s檢驗,其使用格式為ks.test(x, y, ...,
alternative = c("two.sided", "less", "greater"),
exact = null)參數x為數值向量.y或者是數值向量(smirnov檢驗),或者是描述分布函數的字元串(kolmogorov檢驗)....描述特定分布參數的字元串.
alternative為備擇假設選項,取"two.sided"(預設值)表示雙側檢驗;取"less"表示備擇假設為“<”的單側檢驗;取"greater"表示備擇假設為“>”的單側檢驗.
exact為邏輯變量(預設值為null),表示是否精确計算p值.
例1.24 對一台裝置進行壽命檢驗,記錄10次無故障工作時間,并按從小到大的次序排列如下(機關：小時)：
420 500 920 1380 1510 1650 1760 2100 2300 2350
試用kolmogorov-smirnov檢驗方法檢驗此裝置的無故障工作時間是否服從指數分布,且平均無故障工作時間為1500h.
解 指數分布的參數為λ,均值為1/λ,是以,平均時間為1500h,則參數為λ=1/1500.輸入資料,調用ks.test()函數計算(程式名:exam0124.r)，p值(=0.2654)>0.05,無法拒絕原假設,是以,認為此裝置的無故障工作時間服從指數分布,且平均無故障工作時間為1500h.> x <- c( 420, 500, 920, 1380, 1510, 1650, 17<code>`</code>javascript
60, 2100, 2300, 2350)
ks.test(x, "pexp", 1/1500)
one-sample kolmogorov-smirnov test
d = 0.3015, p-value = 0.2654
can("exam0125.data", nlines = 3)
y <- scan("exam0125.data", skip = 3)
ks.test(x, y)
two-sample kolmogorov-smirnov test
data: x and y
d = 0.23, p-value = 0.5286
alternative hypothesis: two-sidedp值(=0.5286)>0.05,接受原假設,即認為兩個總體的分布函數是相同的.
kolmogorov-smirnov檢驗與pearson χ2檢驗相比,不需将樣本分組,少了一個任意性,這是其優點.其缺點是隻能用在理論分布為一維連續分布且分布完全已知的情形,适用面比pearson χ2檢驗小.研究也顯示:在kolmogorov-smirnov檢驗可用的場合下,其功效一般來說略優于pearson χ2檢驗.
1.4.3 正态性檢驗
正态性檢驗是非參數檢驗中最常用的檢驗,除了用前面介紹的兩種方法作檢驗外,還用許多檢驗方法,這裡介紹shapiro-wilk(夏皮羅威爾克)正态性檢驗,它是利用來自總體x的樣本x1,x2,…,xn,檢驗h0:總體x具有正态分布由于在檢驗中用w統計量作正态性檢驗,是以這種檢驗方法也稱為正态w檢驗方法.
在r中,shapiro.test()函數完成shapiro-wilk正态性檢驗,其使用格式為shapiro.test(x)參數x是由樣本構成的向量,并且向量的長度為3～5000.
例1.26 用shapiro-wilk正态性檢驗方法檢驗例1.1中考試成績是否服從正态分布。
解 讀取資料,調用shapiro.test()函數計算(程式名:exam0126.r).x <- scan("exam0101.data"); shapiro.test(x)
shapiro-wilk normality test
w = 0.8633, p-value = 0.0009852p值(=0.0009852)0.05,拒絕原假設,即考試成績不服從正态分布,這與前面的分析(見例1.3和例1.4)是相符的.