R語言系列—-區間估計

這一篇講的是區間估計…..因為這不是一個關於統計學的系列，所以對文中出現的公式不會給予任何證明…..就是這樣。

就從一個最簡單的正態分布的方差已知時，求均值的置信區間開始吧。

書上的公式告訴我們這個區間是 $\overline{x}\pm(\sigma/\sqrt{n})z_{1-\sigma/2}$,其中Z_p表示的是正態分布N(0,1)下側的p分位數。

我們用R來實現求得這一結果的過程。下面設x里存儲了給出的樣本，sigma表示已知的方差，n表示樣本的個數， alpha則是(1-置信水平)

       mean<-mean(x)

       ans<-c(mean-sigma*qnorm(1-alpha / 2)/sqrt(n) , mean+sigma*qnorm(1-alpha / 2)/sqrt(n))

這樣,ans就存儲了要求的置信區間。

來解釋一下吧，先用mean(x)求出樣本的平均值，然后用qnorm(1-alpha / 2)求出Z_1-a/2，（還記得么？前綴q是分位數函數，）剩下的就是套公式的加減法了。

這里的qnorm(1-alpha / 2)其實省略了很多參數，完整一些的寫法是

qnorm(1-alpha/2,mean=0,sd=1,lower.tail=TRUE)

第一個參數就不用解釋了，第二,三個參數mean=0,sd=1，表示這是一個標准正態分布(不同於前面，這里增加了mean=和sd=，這種做法的好處是可以改變參數的順序，但是結果是一樣的)，最后一個參數lower.tail這個參數的意思就比較有意思了，官方解釋如下：

 if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x].

明白了么?等於真的話，得出的就是X<=x的分位數，為假的話就是從X>x的方法尋找這個值。一般我們用默認的真就可以了。

接下來我們把它整理成一個函數，方便使用

z.test<-function(x,n,sigma,alpha){

mean<-mean(x)

ans<-c(

       mean-sigma*qnorm(1-alpha/2,mean=0,sd=1,lower.tail=TRUE)/sqrt(n),

       mean+sigma*qnorm(1-alpha/2,mean=0,sd=1,lower.tail=TRUE)/sqrt(n))

ans

}

這樣我們就可以直接使用z.test()完成對u的置信區間的計算。

比如，有10個樣本，分別是175,176,173,175,174,173,173,176,173,179。標准差為1.5，求均值95%的置信區間：

x<-c(175,176,173,175,174,173,173,176,173,179)

z.test(x,10,1.5,0.05)

則返回置信區間：

[1] 173.7703 175.6297