这一篇讲的是区间估计…..因为这不是一个关于统计学的系列,所以对文中出现的公式不会给予任何证明…..就是这样。
就从一个最简单的正态分布的方差已知时,求均值的置信区间开始吧。
书上的公式告诉我们这个区间是 $\overline{x}\pm(\sigma/\sqrt{n})z_{1-\sigma/2}$,其中Zp表示的是正态分布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)求出Z1-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