R語(yǔ)言與區(qū)間估計(jì)學(xué)習(xí)筆記

鑒于區(qū)間估計(jì)的理論與方法可以在任意一本統(tǒng)計(jì)學(xué)教程中找到，故這里只是單純的介紹R語(yǔ)言中區(qū)間估計(jì)的函數(shù)與一些自己編寫的區(qū)間估計(jì)函數(shù)。

一、單正態(tài)總體的參數(shù)估計(jì)

1、 方差已知時(shí)的均值估計(jì)

z.test<-function(x,n,sigma,a,u0,alt){

result<-list()

mean<-mean(x)

result$interval<-c(mean-sigma*qnorm(1-a/2,0,1)/sqrt(n),mean+sigma*qnorm(1-a/2,0,1)/sqrt(n))

z<-(mean-u0)/(sigma/sqrt(n))

p<-pnorm(z,lower.tail=F) #函數(shù)筆記：lower.tail是真的話，得出的就是X<=x的分位數(shù)，為假的話就是用P(X>x)的辦法尋找這個(gè)值。一般我們用默認(rèn)的真就可以了

result$z<-z

result$p.value<-p #通過(guò)P值判定參數(shù)估計(jì)效果

if(alt==2)

reslut$p.value<-2*pnorm(abs(z),lower.tail=F)

else

reslut$p.value<-pnorm(z)

reslut#函數(shù)筆記：如果函數(shù)的結(jié)果需要有多個(gè)返回值，可以創(chuàng)建一個(gè)list()，并返回該對(duì)象。也可以用return()函數(shù)，設(shè)定返回值。但是一個(gè)函數(shù)的返回的對(duì)象只有一個(gè)。

}

2、 方差未知時(shí)的均值估計(jì)

在小樣本中，我們通常使用R語(yǔ)言的內(nèi)置函數(shù)t.test()調(diào)用格式：

t.test(x, y = NULL,

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

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

conf.level = 0.95, ...)

對(duì)于大樣本，我們可以使用樣本方差代替總體方差，使用z.test()處理

3、 方差的區(qū)間估計(jì)

chisq.var.test<-function(x,n,a,alt=2,sigma0=1)

{

result<-list()

v<-var(x)

result$interval<-c((n-1)*v/qchisq(1-a/2,n-1,lower.tail=T),(n-1)*v/qchisq(a/2,n-1,lower.tail=T))

chi2<-(n-1)*v/sigma0

result$chi2<-chi2

p<-pchisq(chi2,n-1)

if(alt==2)

result$p.value<-2*min(pchisq(chi2,n-1),pchisq(chi2,n-1,lower.tail=F))

else

result$p.value<-pchisq(chi2,n-1,lower.tail=F)

result

}

這里雖然用fisher引理知道利用卡方分布來(lái)處理，但是我們不用chisq.test()來(lái)命名這個(gè)函數(shù)，因?yàn)镽的內(nèi)置函數(shù)中有chisq.test().如果我們這樣命名函數(shù)，會(huì)導(dǎo)致卡方檢驗(yàn)時(shí)使用有些許不便。

二、兩正態(tài)總體參數(shù)的區(qū)間估計(jì)

1、 兩方差都已知時(shí)兩均值差的置信區(qū)間

two.sample.sigmaknown<-function(x,y,conf.level=0.95,sigma1,sigma2,alt=c("twosides","less","greater"))

{

n1<-length(x)

n2<-length(y)

x_<-mean(x)-mean(y)

a<-1-conf.level

z1<-qnorm(1-a/2)*sqrt(sigma1/n1+sigma2/n2)

z2<-qnorm(1-a)*sqrt(sigma1/n1+sigma2/n2)

if(alt=="two sides")

x_ +c(-z1,z1)

else if(alt=="less")

x_ -z2

else

x_ +z2

}

注：對(duì)于大樣本，我們可以以樣本標(biāo)準(zhǔn)差代替總體方差來(lái)進(jìn)行區(qū)間估計(jì)

2、 兩方差都未知但相等時(shí)兩均值差的置信區(qū)間

直接使用t.test()函數(shù)即可

注：由于對(duì)于一般情形估計(jì)的方法特別多，可以使用neyman的樞軸量法亦可以使用fisher的信仰推斷(通常認(rèn)為后者較好)。故在此不予介紹

3、 兩方差比的置信區(qū)間 數(shù)據(jù)分析培訓(xùn)

仔細(xì)閱讀方差比的區(qū)間估計(jì)內(nèi)容，我們應(yīng)該注意到，兩樣本在做方差比估計(jì)時(shí)應(yīng)該需要做正態(tài)性檢驗(yàn)，在R中，你可以使用函數(shù)shapiro.test()來(lái)實(shí)現(xiàn)，該檢驗(yàn)對(duì)數(shù)據(jù)的正態(tài)性是給與保護(hù)的。

使用var.test()函數(shù)，調(diào)用格式如下：

var.test(x, y, ratio = 1,

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

conf.level = 0.95, ...)

最后，我想解釋一下置信水平的含義，我們所說(shuō)的置信水平是指用這樣的辦法對(duì)數(shù)據(jù)進(jìn)行100次估計(jì)，包含真值的次數(shù)為100*conf.level。這里的估計(jì)是指對(duì)不同數(shù)據(jù)用同樣方法進(jìn)行估計(jì)。我們可以編寫一個(gè)R函數(shù)來(lái)驗(yàn)證一下：

judge<-rep(0,1000)

for(i in 1:1000){

set.seed(5*i)

if(t.test(rnorm(100000,5,17))$conf.int[1]<=5&5<=t.test(rnorm(100000,5,17))$conf.int[2])

judge[i]<-0

else

judge[i]<-1

}

table(judge)

輸出結(jié)果：

#judge

# 0 1

#954 46

#從這里來(lái)看，估計(jì)達(dá)到95%的置信水平