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