《數學模組化：基于R》一一

本節書摘來自華章計算機《數學模組化：基于r》一書中的第1章，第1.2節,作者：薛　毅 更多章節内容可以通路雲栖社群“華章計算機”公衆号檢視。

參數的區間估計和假設檢驗可以看成一個問題的兩個方面.區間估計是用兩個統計量所構成的區間來估計一個未知的參數,并同時指明此區間可以覆寫住這個參數的可靠程度（置信度）.假設檢驗是統計推斷中的一個重要内容,它是利用樣本資料對某個事先做出的統計假設按照某種設計好的方法進行檢驗,判斷此假設是否正确.

最常用的區間估計和假設檢驗是針對正态分布而言的,是以,本節着重介紹正态總體下參數的區間估計與假設檢驗的内容.

1.2.1　單個總體的區間估計與假設檢驗

假設總體滿足x~n(μ,σ2),x1,x2,…,xn為來自總體x的一個樣本,樣本均值x=1n∑ni=1xi滿足x~nμ,σ2n,是以z=x-μσ/n~n(0,1)（1.19）即随機變量z為标準正态分布.樣本方差s2=1n-1∑ni=1(xi-x)2滿足(n-1)s2σ2~χ2(n-1)（1.20）如果x與s2互相獨立,則由式(1.19)、式(1.20)和t分布的定義,得到t=x-μσ/n(n-1)s2σ2/(n-1)=x-μs/n~t(n-1)（1.21）即随機變量t具有自由度為n-1的t分布.

1.均值μ的區間估計

對于正态總體來說,有兩個參數,一個是均值參數μ,另一個是方差參數σ2.在作均值μ的區間估計時,需要假設σ2已知或未知.為了與實際情況相吻合,這裡隻讨論σ2未知的情況.

在σ2未知的情況下,由式(1.21),得到px-μs/n≤tα/2(n-1)=1-α（1.22）其中tα(n-1)表示自由度為n-1的t分布的上α分位點.由式(1.22)得到均值μ的置信水準為1-α的雙側置信區間x-sntα/2(n-1),　x+sntα/2(n-1)（1.23）類似地,由式(1.21),有px-μs/n≤tα(n-1)=1-α,　p-tα(n-1)≤x-μs/n=1-α（1.24）得到μ的置信水準為1-α的單側置信區間x-sntα(n-1),+∞,　-∞,x+sntα(n-1)（1.25）2.均值μ的假設檢驗

在統計假設檢驗中,需要有一個作為檢驗的對象的假設,常稱為原假設或零假設,記為h0.與之相應,為使問題表述得更明确,還常提出一個與之對應的假設,稱為對立假設或備擇假設,記為h1.

與雙側區間估計相對應的檢驗是雙側檢驗,也稱為雙邊檢驗,這裡檢驗h0:μ=μ0,　h1:μ≠μ0如果h0為真,當σ2未知時,由式(1.21),得到t=x-μ0s/n~t(n-1)（1.26）是以當t>tα/2(n-1)時,拒絕原假設.是以,稱(-∞,-tα/2(n-1))和(tα/2(n-1),+∞)為拒絕域.

對應單側區間估計的檢驗是單側檢驗,也稱為單邊檢驗,即檢驗h0:μ≤μ0,　h1:μ>μ0

(h0:μ≥μ0　h1:μ<μ0)此時,當t>tα(n-1)(t<-tα(n-1))時,拒絕原假設.對應的拒絕域為(tα(n-1),+∞)　((-∞,-tα(n-1)))上述檢驗方法是通過t統計量構造的,是以,稱該方法為t檢驗.

1.2.2　兩個總體的區間估計與假設檢驗

設x1,x2,…,xn1是來自總體x～n(μ1,σ21)的樣本,y1,y2,…,yn2是來自總體y~n(μ2,σ22)的樣本,并且兩樣本獨立.這裡隻讨論σ2i未知的情況,令x和y分别為兩組樣本的均值,s21和s22分别為兩組樣本的方差.

如果兩個總體的方差相同,即σ21=σ22.類似于單個總體的推導,有t=x-y-(μ1-μ2)sw1n1+1n2~t(n1+n2-2)（1.27）其中sw=(n1-1)s21+(n2-1)s22n1+n2-2（1.28）當兩個總體的方差不同(σ21≠σ22)時,有t=x-y-(μ1-μ2)s21n1+s22n2~t(ν)（1.29）近似成立,其中ν=s21n1+s22n22(s21)2n21(n1-1)+(s22)2n22(n2-1)(1.30)1.均值差μ1-μ2的區間估計

利用式(1.27),類似于單個總體置信區間的推導,當兩個總體的方差相同時,置信水準為1-α,均值差μ1-μ2的雙側置信區間為x-y-tα/2(n1+n2-2)sw1n1+1n2，　x-y+tα/2(n1+n2-2)sw1n1+1n2

(1.31)單側置信區間為x-y-tα(n1+n2-2)sw1n1+1n2,+∞(1.32)

-∞,x-y+tα(n1+n2-2)sw1n1+1n2(1.33)當兩個總體的方差不同時,利用式(1.29),可以得到置信水準為1-α,均值差μ1-μ2的雙側置信區間為x-y-tα/2(ν)s21n1+s22n2,　y-x+tα/2(ν)s21n1+s22n2(1.34)單側置信區間為x-y-tα(ν)s21n1+s22n2,+∞,　-∞,y-x+tα(ν)s21n1+s22n2(1.35)其中ν由式(1.30)得到.

2.均值差μ1-μ2的假設檢驗

均值差μ1-μ2的雙側檢驗為h0:μ1-μ2=0,　h1:μ1-μ2≠0當h0為真時,由式(1.27)(兩總體方差相同的情況),得到t=x-ysw1n1+1n2～t(n1+n2-2)（1.36）是以,當t>tα/2(n1+n2-2)時,拒絕原假設.對應的拒絕域為(-∞,-tα/2(n1+n2-2))　和　(tα/2(n1+n2-2),+∞)由式(1.29)(兩總體方差不同的情況),得到t=x-ys21n1+s22n2～t(ν)（1.37）近似成立.是以,當t>tα/2(ν)時,拒絕原假設.對應的拒絕域為(-∞,-tα/2(ν))　和　(tα/2(ν),+∞)對于均值差μ1-μ2的單側檢驗,h0:μ1-μ2≤0,　h1:μ1-μ2>0

(h0:μ1-μ2≥0,　h1:μ1-μ2<0)在方差相同的情況下,當t>tα(n1+n2-2)(t<-tα(n1+n2-2))時,拒絕原假設.對應的拒絕域為(tα(n1+n2-2),+∞)　（(-∞,-tα(n1+n2-2)）)在方差不同的情況下,當t>tα(ν)(t<-tα(ν))時,拒絕原假設.對應的拒絕域為(tα(ν),+∞)　（(-∞,-tα(ν))）此方法仍稱為t檢驗法.

3.成對資料均值差的區間估計與假設檢驗

如果資料是成對出現的,即(xi,yi)(i=1,2,…,n),則可以進行成對資料均值差的區間估計和假設檢驗,其估計和檢驗方法就是令zi=xi-yi(i=1,2,…,n),然後對z作單個總體均值的區間估計與假設檢驗.

1.2.3　區間估計與假設檢驗的計算

p值

在前面的論述中,需要用拒絕域來否定h0,即在計算完t統計量後,需要計算相應的臨界值.但在計算軟體中,通常會使用更友善的方法,也就是通過計算p值來作判斷.所謂p值,就是犯第一類錯誤的機率,即p值=p{否定h0h0為真}當p值<α時,拒絕原假設;否則,接受原假設.事實上,p值<α等價于t統計量屬于拒絕域,是以,兩種方法是等價的.

r中的t檢驗函數

在r中,用t.test()函數完成t檢驗的工作,并給出相應的置信區間,其使用格式為t.test(x, y = null,

　　alternative = c("two.sided", "less", "greater"),

　　mu = 0, paired = false, var.equal = false,

　　conf.level = 0.95, ...)參數x，y均為樣本構成的數值向量,對于單個總體的樣本,y為null(預設值).

alternative為備擇假設選項,取"two.sided"(預設值)表示雙側檢驗，取"less"表示備擇假設為“<”的單側檢驗，取"greater"表示備擇假設為“>”的單側檢驗.

mu為數值,表示原假設μ0,預設值為0.

paired為邏輯變量,用于說明是否完成成對資料的t檢驗.當資料是成對資料時,此參數取true;否則取false(預設值).

var.equal為邏輯變量,用于說明兩個樣本總體的方差是否相同.當兩總體的方差相同時,此參數取true,否則取false(預設值).

conf.level為0～1之間的數值(預設值為0.95),表示置信水準,它将用于計算均值μ的置信區間.

...為附加參數.

另一種使用格式是公式形式,其使用格式為t.test(formula, data, subset, na.action, ...)用于兩個總體樣本的檢驗,參數formula為形如value ~ group的公式,其中value為資料,group為資料的分組情況,通常是因子向量.

data為矩陣或資料框.subset為可選向量,表示使用的樣本子集.

na.action為函數,表示樣本中出現缺失值(na)的處理方法,預設值為函數getoption("na.action").

例1.5　在某超市随機地選取了某品牌标重為10kg的大米10袋,稱得重量如下：

10.1,10,9.8,10.5,9.7,10.1,9.9,10.2,10.3,9.9

假設所稱出的重量服從正态分布,試分析該規格的大米重量是否與标重相符?

解　可以認為這10袋大米是來自總體的樣本,且總體服從正态分布,總體的均值為μ,方差未知,是以選擇t檢驗.根據題意,需要檢驗h0:μ=10,　h1:μ≠10程式(程式名:exam0105.r)和計算結果如下：> x <- c(10.1, 10, 9.8, 10.5, 9.7, 10.1, 9.9, 10.2,10.3, 9.9)

t.test(x, mu = 10)

　　　　one sample t-test

data: x

t = 0.6547, df = 9, p-value = 0.5291

alternative hypothesis: true mean is not equal to 10

95 percent confidence interval:

　9.877225 10.222775

sample estimates:

mean of x

　　10.05在程式中,x是10袋大米重量構成的向量,mu =10表示μ0=10.其他參數選用預設值,即作雙側檢驗,置信水準為0.95.

程式的計算結果表示:單個總體樣本的檢驗,資料由變量x輸入,t統計量為0.6547,也就是式(1.26)中的t值.自由度為9,p值為0.5291,是以接受原假設,說明該規格的大米重量與标重是相符的.

置信水準為95％的置信區間是［9.877,10.223］,這說明此規格大米的平均重量有95％的可能在9.877kg到10.223kg之間.樣本均值是10.05,說明這10袋大米的平均重量超過10kg.

例1.6　某種元件的壽命x（機關:h）服從正态分布n(μ,σ2),其中μ和σ2均未知.現測得16隻元件的壽命如下：

159　280　101　212　224　379　179　264

222　362　168　250　149　260　485　170

問是否有理由認為元件的平均壽命大于225h?

分析　該問題是作雙側檢驗還是作單側檢驗?實際上,并不需要考慮元件的壽命是不是等于225h,而是希望證明其壽命大于225h,是以雙側檢驗并不合适,應該作單側檢驗.

對于單側檢驗,是檢驗h0:μ≥μ0=225,　h1:μ<μ0=225(1.38)還是檢驗h0:μ≤μ0=225,　h1:μ>μ0=225？(1.39)從表面上看,好像應該作檢驗(1.38),但事實上,作檢驗(1.39)更合理.因為當原假設與需要得到的結論相反時,利用拒絕原假設,所得到的結論才更有說服力.

與原假設相反,就是與備擇假設一緻,可以從這方面了解為什麼t.test()函數輸入參數是備擇假設而不是原假設.

解　輸入資料,調用函數t.test()(程式名:exam0106.r)得到> x <- c(159, 280, 101, 212, 224, 379, 179, 264,

222, 362, 168, 250, 149, 260, 485, 170)

t.test(x, alternative = "greater", mu = 225)

t = 0.6685, df = 15, p-value = 0.257

alternative hypothesis: true mean is greater than 225

　198.2321　　inf

　　241.5在函數的輸出中,有t統計量(t)、自由度(df)、p值(p-value)和均值μ的置信區間,以及μ的估計值.由于這裡作的是單側檢驗,是以給出的置信區間也是單側的.

由于p值(=0.257)>0.05,不能拒絕原假設,接受h0,隻能認為這批元件的平均壽命不大于225h.

例1.7　從一批燈泡中随機地取5隻作壽命試驗,測得壽命（機關:h）為

1050,1100,1120,1250,1280

設燈泡壽命服從正态分布.(1)求燈泡壽命平均值的置信水準為0.95的單側置信下限;(2)求這批燈泡平均壽命大于1000h的機率.

解　計算單側置信下限,原假設應取“小于”,也就是說,備擇假設應取“大于”.同時打算計算平均壽命大于1000h的機率,是以,取μ0=1000.

輸入資料,調用函數t.test(),得到> x <- c(1050, 1100, 1120, 1250, 1280)

t.test(x, mu = 1000, al = "g")

t = 3.5867, df = 4, p-value = 0.01151

alternative hypothesis: true mean is greater than 1000

　1064.9　　inf

　　1160在程式中, al是alternative的縮寫,"g"是"greater"的縮寫.

注意:在r中,很多地方都可以使用縮寫,這樣便于程式的錄入.

單側置信下限為1064.9,即約有95％的燈泡能使用1065h以上.按照p值的意義,是以這批燈泡的平均壽命大于1000h的機率為1-0.0115=0.9885.

例1.8　在平爐上進行一項試驗以确定改變操作方法的建議是否會增加鋼的得率,試驗是在同一個平爐上進行的.每煉一爐鋼時除操作方法外,其他條件都盡可能做到相同.先用标準方法煉一爐,然後用新方法煉一爐,以後交替進行,各煉了10爐,其得率如表1.1所示.設這兩個樣本互相獨立,并且分别來自正态總體n(μ1,σ21)和n(μ2,σ22),其中μ1,μ2和σ21,σ22未知.問新的操作能否提高得率?（取α=0.05.）

表1.1　兩種方法的得率标準方法新方法标準方法新方法178.179.1678.479.1272.481.0776.079.1376.277.3875.577.3474.379.1976.780.2577.480.01077.382.1

解　作雙側檢驗h0:μ1=μ2,　h1:μ1≠μ2分别讨論σ21=σ22和σ21≠σ22的情況.計算過程如下：## %輸入資料

x <- c(78.1, 72.4, 76.2, 74.3, 77.4, 78.4, 76.0,75.5, 76.7, 77.3) y <- c(79.1, 81.0, 77.3, 79.1, 80.0, 79.1, 79.1,77.3, 80.2, 82.1)

t.test(x, y, var.equal = true)

two sample t-test

data: x and y

t = -4.2957, df = 18, p-value = 0.0004352

alternative hypothesis:

　　true difference in means is not equal to 0

-4.765026 -1.634974

mean of x mean of y

　　76.23　　79.43

t.test(x, y)

　　　　welch two sample t-test

t = -4.2957, df = 17.319, p-value = 0.000471

-4.769446 -1.630554

　　76.23　　79.43無論是方差相同模型還是方差不同模型,計算出的p值均小于0.05,拒絕原假設,即說明μ1≠μ2.進一步,從均值差μ1-μ2的區間估計來看,兩種方法得到的置信區間均小于0,說明μ1-μ2<0,即μ1<μ2,是以新的操作能夠提高得率.

那麼,兩種模型(方差相同與方差不同)哪一個更好呢?從p值和置信區間來看,方差相同模型的效果更好,因為它的p值較小,置信區間的長度也較短.

此問題也可以作單側檢驗,會得到同樣的檢驗結果,關于單側檢驗的計算留給讀者完成.

例1.9　為了調查應用克矽平治療矽肺的效果,今抽查應用克矽平治療矽肺的患者10名,記錄下治療前後血紅蛋白的含量資料,如表1.2所示.試求治療前後變化的區間估計(α=0.05).

表1.2　治療前後血紅蛋白的含量資料病人編号治療前治療後病人編号治療前治療後111.314.0610.012.0215.013.8711.014.7315.014.0812.011.4413.513.5913.013.8512.813.51012.312.0

解　此類資料可以看成是成對資料,因為每對資料是每個病人在治療前和治療後對比的資料.計算過程如下：> x<-c(11.3, 15.0, 15.0, 13.5, 12.8, 10.0, 11.0,12.0, 13.0, 12.3)

y<-c(14.0, 13.8, 14.0, 13.5, 13.5, 12.0, 14.7,11.4, 13.8, 12.0) t.test(x, y, paired = true)

　　　　paired t-test

t = -1.3066, df = 9, p-value = 0.2237

　　95 percent confidence interval:

　　-1.8572881　0.4972881

　　sample estimates:

　　mean of the differences

　　　　　　　　　　　　-0.68p值為0.2237>0.05,無法拒絕原假設,即μ1=μ2,也就是說,治療前和治療後沒有顯著性差異.另外,得到的置信區間包含0,這也說明治療前後沒有顯著性差異.

讀者也可以選擇方差相同模型和方差不同模型進行計算,并分析哪種模型效果最好及其原因。

例1.10　為了分析霧天與晴天對果園使用殺蟲劑的效果,現收集了12天噴灑殺蟲劑時果園的空氣樣本,如表1.3所示.利用這些資料檢驗空氣中的氧硫比值在霧天與晴天之間是否存在顯著差異.(α=0.05.)

表1.3　噴灑殺蟲劑時果園的空氣資料(機關:ng/m3)環境硫氧環境硫氧1霧38.210.37晴46.427.42霧28.66.98霧135.944.83霧30.26.29霧102.927.84霧23.712.410晴28.96.55晴74.145.811霧46.911.26霧88.29.912晴44.316.6

解　本例的資料适合用資料框的形式輸入,是以用公式格式計算.airdata <- data.frame(

　　e = factor(c("fog", "fog", "fog", "fog", "air", "fog",

　　　　　　　　　"air", "fog", "fog", "air", "fog", "air")),

　　s = c( 38.2, 28.6, 30.2, 23.7, 74.1, 88.2, 46.4, 135.9,

　　　　　102.9, 28.9, 46.9, 44.3),

　　o = c(10.3, 6.9, 6.2, 12.4, 45.8, 9.9, 27.4, 44.8,

　　　　　27.8, 6.5, 11.2, 16.6)

)

airdata$r <- airdata$o/airdata$s

t.test(r ~ e, data = airdata, var.equal = true)在程式中,data.frame表示資料為資料框結構,factor為因子函數,說明變量是因子.這裡采用方差相同模型,計算結果如下：　　　　　　two sample t-test

data: r by e

t = 2.0446, df = 10, p-value = 0.06812

-0.01600593　0.37255524

mean in group air mean in group fog

　　　　0.4520581　　0.2737834計算結果表明,在α=0.05的條件下,并不能說明霧天與晴天空氣中的氧硫比有顯著差異.

1.2.4　兩個正态總體方差比σ21/σ22的估計與檢驗

在兩個正态總體的計算中,要判斷兩總體的方差是否相同,通過方差比的估計和檢驗,可以完成這項工作.

1.方差比σ21/σ22的區間估計

設x1,x2,…,xn1是來自總體x～n(μ1,σ21)的樣本,y1,y2,…,yn2是來自總體y～n(μ2,σ22)的樣本.令x和y分别為來自總體x和y的樣本均值,s21和s22分别為來自總體x和y的樣本方差.

當兩樣本獨立時,有f=s21/σ21s22/σ22=s21/s22σ21/σ22～f(n1-1,n2-1)（1.40）是以pf1-α/2(n1-1,n2-1)≤s21/s22σ21/σ22≤fα/2(n1-1,n2-1)=1-α(1.41)則σ21/σ22的置信水準1-α的置信區間為s21/s22fα/2(n1-1,n2-1),　s21/s22f1-α/2(n1-1,n2-1)（1.42）2.方差比σ21/σ22的假設檢驗

兩個總體方差比的檢驗是檢驗h0:σ21/σ22=1,　h1:σ21/σ22≠1當h0為真時,由式(1.40),有f=s21s22～f(n1-1,n2-1)（1.43）當f≥fα/2(n1-1,n2-1)或者f≤f1-α/2(n1-1,n2-1)時,拒絕原假設.對應的拒絕域為［fα/2(n1-1,n2-1),+∞)和［0,f1-α/2(n1-1,n2-1)］.

類似于前面的讨論,方差比也可以作單側檢驗h0:σ21/σ22≤1,　h1:σ21/σ22>1

(h0:σ21/σ22≥1,　h1:σ21/σ22<1)上述檢驗是根據f統計量完成的檢驗,是以稱為f檢驗.由于是判斷方差是否相同,是以也稱為方差齊性檢驗.

3.計算

在r中,用var.test()函數作f檢驗,其使用格式為var.test(x, y, ratio = 1,

　　conf.level = 0.95, ...)參數x和y分别為兩個樣本構成的數值向量.ratio為兩個總體的方差比,預設值為1.

conf.level為0～1之間的數值(預設值為0.95),表示置信水準,它将用于計算方差比σ21/σ22的置信區間.

另一種使用格式是公式形式,其使用格式為var.test(formula, data, subset, na.action,...)參數formula為形如value ~ group的公式,其中value為資料,group為資料的分組情況,通常是因子向量.

data為矩陣或資料框.subset為可選向量,表示使用樣本的子集.

例1.11　試對例1.8中的資料進行假設檢驗：h0:σ21/σ22=1,　h1:σ21/σ22≠1解　調用var.test()函數:> source("exam0108.r"); var.test(x, y)

　　　　f test to compare two variances

f = 1.4945, num df = 9, denom df = 9, p-value = 0.559

　　true ratio of variances is not equal to 1

　0.3712079 6.0167710

ratio of variances

　　　　1.494481在函數的輸出中,有f統計量(f)、第1自由度或分子自由度(num df)、第2自由度或分母自由度(denom df)、p值(p-value)和方差比σ21/σ22的置信區間,以及f比.

由于p值(=0.559)0.05,無法拒絕原假設,認為兩總體的方差是相同的.從方差比的置信區間［0.37,6.02］來看,它包含1,也就是說,有可能σ21/σ22=1,是以認為兩總體的方差是相同的.是以,在例1.8中,使用兩總體方差相同模型的計算效果更好是有道理的.