第3章参数估计

 掌握参数估计的原理。
 理解总体方差、总体比例的区间估计。
 掌握统计量的分布: χ2分布、t分布和F分布。
 能运用参数估计的区间估计进行R编程计算。
 学会进行ShapiroWilk检验。






根据样本统计量来推断所关心的总体参数是非常必要的。本章重点介绍参数估计，这是推断统计的重要内容之一。通过统计描述，可以对样本数据的情况进行详细了解，但是引进统计量的真正目的是对感兴趣的问题进行统计推断，而在实际中，人们感兴趣的问题往往是与未知参数有关的。考察样本所代表的总体情况如何，也就是说，根据样本提供的信息来推断总体的特征才是研究的重点所在，涉及以下两个问题： 
问题1： 不同情况下如何根据样本统计量来估计总体的参数？
问题2： 如何确定参数估计问题的样本量？ 
3.1统计量的分布
3.1.1总体与样本
样本是从总体中抽出的部分单位的集合，这个集合中所包含的单位个数称为样本容量，一般用n表示。所研究对象的全体称为总体。研究只有一个特征的总体时，总体的这个指标(变量)可以看作一维随机变量。例如，生产线上生产的零件的直径、某年某专业考生的成绩等。
只要反映总体特征的变量的取值满足以下3个条件，就是一个随机变量。
(1) 在同一条件下可无限次重复的取值。 
(2) 取值的可能结果有多个，且不确定。 
(3) 事前不知取值结果。
满足上述3条，其实就是一个随机实验。反映总体特征的随机变量的取值的全体，也称为总体。这个总体，其实就是样本空间。反映总体特征的随机变量的概率分布，称为总体分布。
随机变量是表征一个随机试验结果的变量，其数值由一次试验结果所决定，但是在试验之前是不确定的(取值随机而定)。随机变量的所有可能取值就是所有基本事件对应的值(“值”不一定是数值，可以是字符串)，通常用英文大写字母或希腊字母表示，随机变量又分为连续型和离散型随机变量。分位数是统计推断中常用的一类数字特征，也称为分位点，主要用于连续型随机变量。给定0＜α＜1和随机变量X的概率密度函数f(x)，显然有： ∫+∞-∞f(x)dx=1。若从某个值开始，f(x)与x轴之间的面积等于α，这个值就是分位数Fα。
常用的分位数有上侧分位数和双侧分位数。对于给定的α(0＜α＜1)，满足∫∞Fαf(x)dx=α的Fα称为该分布的α水平上侧分位数。∫∞Fαf(x)dx=α等价于P(X＞Fα)=α。
样本个数又称样本可能数目，它指从一个总体中可能抽取的多少个样本。






3.1.2统计量的分布
正态分布是特别重要的连续型分布。一般说来，若某一数量指标受到很多相互独立的随机因素的影响，而每个因素所起的作用都很微小，则这个数量指标近似服从正态分布。许多自然现象和社会现象可用正态分布来描述。传统的统计理论主要建立在正态分布基础上。由正态分布出发，导出了一系列重要的抽样分布，如t分布、χ2分布、F分布等。
若随机变量X的概率密度为f(x)=12πσexp-12x-μσ2(-∞＜x＜∞)，则称X服从正态分布，记作： X~N(μ,σ2)。特别地，标准正态分布N(0,1)的概率密度函数为f(x)=12πexp-x22(-∞＜x＜∞)。 
由标准正态分布的随机样本所引出的重要的统计量的分布有χ2分布、t分布和F分布。
1.  χ2(n)分布
设随机变量X服从N(0,1)分布，X1，X2，…，Xn为X的样本，则χ2=X21+X22+…+X2n服从自由度为n的χ2分布。记为χ2~χ2(n)。
若X~χ2(n)，则E(X)=n，Var(X)=2n。
统计学上的自由度是指当以样本的统计量来估计总体的参数时，样本中独立或能自由变化的资料的个数，称为该统计量的自由度。χ2(n)分布的形状随自由度n的变化而变化，当n较小时，分布不对称； 当n增大时，分布逐渐趋向于正态分布。
例如，在估计总体的平均数时，样本中的个数全部加起来，其中任何一个数都和其他资料相独立，从其中抽出任何一个数都不影响其他资料(这也是随机抽样所要求的)。 因此一组资料中每一个资料都是独立的，所以自由度就是估计总体参数时独立资料的数目，而平均数是根据独立资料来估计的，因此自由度为n。
2. t分布
若X~N(0,1)，Y~χ2(n)，且X与Y相互独立，则t服从自由度为n的分布，t=XY/n~t(n)。
若X~t(n)，则E(X)=0，Var(X)=nn－2。t(n)分布的形状关于x=0对称，当自由度n很大时，t分布趋向于标准正态分布。
3. F分布
若X~χ2(n)，Y~χ2(m)且相互独立，则F=X/nY/m~F(n,m)。若X~F(n,m)，则E(X)=mm－2。F分布的曲线是右偏型的。并且随着n、m取值的变小，曲线的偏斜程度越来越大。
1) 由一般正态分布的随机样本所构成的若干重要统计量的分布 
如果X1，X2，…，Xn是正态总体N(μ,σ2)的一个随机样本，则样本均值函数与样本方差函数，满足如下性质。 
(1) 服从Nμ,σ2n分布。 
(2) z=-μσ/n服从N(0,1)分布。 
(3) (n-1)S2σ2服从χ2(n-1)分布。 
(4) t=-μS/n服从t(n-1)分布。 
(5) F=S21/σ21S22/σ22服从F(n1-1,n2-1)分布。

其中，S21是容量为n1的X的样本方差，S22是容量为n2的Y的样本方差。
2) 大样本均值函数的分布： 中心极限定理
设随机变量X服从任何均值为μ、标准差为σ的分布，是随机样本X1，X2，…，Xn的均值函数。
中心极限定理： 当n充分大时，近似地服从均值为μ、标准差为σn的正态分布。
中心极限定理要求n充分大，那么多大才叫充分大呢？这与总体的分布形状有关，总体偏离总体越远，则要求n越大。然而在实际应用中，总体的分布是未知的。本书约定n≥30表示大样本，n＜30表示小样本。这只是一种经验说法。顺便指出，大样本、小样本之间并不是以样本量大小来区分的。中心极限定理的另一种表述是，当n充分大时，-μσ/n近似地服从标准正态分布。
3.2参数估计的基本原理
3.2.1估计量与估计值
参数估计就是根据样本构造适当的统计量来估计总体的参数。在参数估计中，用来估计总体参数的统计量称为估计量，如样本均值、样本比例、样本方差等。样本均值就是总体均值的一个估计量。根据一个具体的样本计算出来的估计量的数值称为估计值。 
在实际中，用什么样的估计量来估计参数呢？实际上没有硬性限制。任何统计量，只要人们觉得合适就可以当成估计量。最常用的估计量是样本均值、样本标准差和Bernoulli实验的成功比例(即样本百分比)； 用它们来分别估计总体均值、总体标准差和Bernoulli实验的成功概率(或总体百分比)，而且分别是相应总体参数的“好”的估计量。
3.2.2点估计与区间估计
参数估计的方法有点估计和区间估计两种。
1. 点估计
点估计就是选择一个适当的样本统计量值作为总体参数的估计值，即用样本统计量的某一个取值θ^直接作为总体参数θ的估计值。
用于估计总体参数θ的估计量有很多，那么具体选用样本的哪种估计量作为总体参数的估计呢？就需要制定评价估计量的标准，统计学家给出了无偏性、有效性、一致性共3个原则。
(1) 无偏性： 设θ^是未知参数θ的估计，有E(θ^)=θ，则称θ^为θ的无偏估计，否则称有偏估计。因为E()=E(X)，总体样本均值是总体均值E(X)的无偏估计。因E(S2)=D(X)，总体样本方差S2是总体方差D(X)的无偏估计。
(2) 有效性： 一个方差较小的无偏估计量称为一个更有效的估计量。如果θ^1和θ^2都是θ的无偏估计，若D(θ^1)＜D(θ^2)，则称θ^1较θ^2有效。
(3) 一致性(相容性)： 被认为是对估计的一个最基本的要求，是指随着样本容量的增大，估计量越来越接近被估计的总体参数。如果一个估计量在样本量不断增大时，它都不能把被估参数估计到任意指定的精度，那么这个估计值是很值得怀疑的。通常，不满足一致性要求的估计一般不予考虑。
【点估计的实例： 点估计在战争中的一个经典实例】
在第二次世界大战期间，盟军非常想知道德军总共制造了多少辆坦克。德国人在制造坦克时是墨守成规的，他们把坦克从1开始进行了连续编号。在战争过程中，盟军缴获了一些敌军坦克，并记录了它们的生产编号。那么怎样利用这些号码来估计坦克总数呢？在这个问题中，总体参数是未知的坦克总数N，而缴获坦克的编号则是样本。
假设盟军手下的一位统计人员负责解决这个问题。针对该调查数据，得知制造出来的坦克总数肯定大于或等于记录的最大编号。为了找到它比最大编号大多少，先找到被缴获坦克编号的平均值，并认为这个值是全部编号的中点。因此样本均值乘以2就是总数的一个估计； 当然要特别假设缴获的坦克代表了所有坦克的一个随机样本。这种估计N的公式的缺点是，不能保证均值的2倍一定大于记录中的最大编号。N的另一个点估计公式是，用观测到的最大编号乘以因子1+1/n，其中n是被缴获的坦克个数。假如缴获了10 辆坦克，其中最大编号是50，那么坦克总数的一个估计是(1+1/10)50=55。此处认为坦克的实际数略大于最大编号。
从战后发现的德军记录来看，盟军的估计值非常接近所生产的坦克的真实值。记录仍然表明统计估计比通常通过其他情报方式做出估计要大大接近于真实数目。
参数的点估计给出了一个具体的数值，方便计算和使用，但是无法给出估计值接近总体参数程度的信息。由于样本是随机的，抽出一个具体的样本得到的估计值很可能不同于总体真值。一个点估计量的可靠性是由它的抽样标准误差来衡量的，这表明一个具体的点估计值无法给出估计的可靠性的度量。 
实际问题中，度量一个点估计的精度的最直观的方法就是给出未知参数的一个区间，这便产生了区间估计的概念。
2. 区间估计
当描述一个人的体重时，一般可能不会说这个人是76.35kg。人们会说这个人是七八十千克，或者是在70~80kg。这个范围就是区间估计的例子。
区间估计是在点估计的基础上，给出总体参数估计的一个区间范围，该区间通常是由样本统计量加减估计误差得到的。在区间估计中，由样本统计量所构造的总体参数的估计区间称为置信区间，其中区间的最小值称为置信下限，最大值称为置信上限。置信度是表示区间估计的可靠程度或把握程度，即所估计的区间包含总体参数真实值的可能性大小，一般以1-α表示，含义是，在同样的方法得到的所有置信区间中，有(1-α)%的区间包含总体参数。其
中，α表示显著性水平，即参数不落在区间内的概率。如图31所示为区间估计示意图。


图31区间估计示意图


3.3总体的区间估计
3.3.1用R进行总体均值的区间估计
由于在R软件中没有求方差已知条件下均值置信区间的内置函数，于是需要编写函数。
z.test()函数的定义如代码清单31所示。 

代码清单31

z.test<-function(x,n,sigma,alpha,u0=0,alternative="two.sided"){

options(digits=4)

result<-list()

mean<-mean(x)

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

p<-pnorm(z,lower.tail=FALSE)

result$mean<-mean

result$z<-z

result$p.value<-p

if(alternative=="two.sided")

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

else if (alternative=="greater")

result$p.value<-pnorn(z)

result$conf.int<-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))

result

}

利用上述程序就可以求出总体均值的置信区间。上述t.test()函数的定义同时完成了区间估计和假设检验，这是为了与R中的t检验函数相对应。如果从上述程序中抽取出区间估计的部分，则得到代码清单32中求置信区间的程序。

代码清单32

conf.int<-function(x,n,sigma,alpha){

options(digits=4)

mean<-mean(x)

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))

}

下面举例说明如何用R软件求解置信水平为1-α的置信区间。
例31
一个人10次称自己的体重(单位： 斤(1斤＝500g))，得到的数值分别为175、176、173、175、174、173、173、176、173、179。现在需要估计一下他的体重，假设此人的体重服从正态分布，标准差为1.5，要求体重的置信水平为95%的置信区间。
解： 根据上述函数z.test()编写R程序如下。

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

> result<-z.test(x,10,1.5,0.05)

> result$conf.int

输出结果如下：

［1］ 173.8 175.6

因此，这里得到体重的置信水平为95%的置信区间为［173.8,175.6］。
注意： R语言中的赋值符号一般是由一个尖括号与一个负号组成的箭头形标志, 该符号可以是从左到右的方向，也可以相反。赋值也可以用函数assign( )实现, 还可以用“=”实现，但它们很少使用。
如果运行如下代码：

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

将同时获得假设检验的结果： 

$mean

［1］ 174.7

$z

［1］ 368.3

$p.value

［1］ 0

$conf.int

［1］ 173.8 175.6

方括号中的数字1表示从mean等变量的第一个元素开始显示。其实该命令的功能在这里与函数print( )相似，输出结果与print(mean) 相同。而上面的程序仅提供了区间估计的部分，这相当于执行了如下代码： 

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

> conf.int(x,10,1.5,0.05)

方差未知时，可直接利用R语言的t.test()函数来求置信区间。t.test()的调用格式如下： 

t.test(x,y=NULL,

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

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

conf.level = 0.95, …)

说明： 若仅出现数据x，则进行单样本t检验； 若出现数据x和y，则进行二样本的t检验； alternative=c("two.sided", "less", "greater")用于指定所求置信区间的类型； alternative="two.sided"是默认值，表示求置信区间； alternative="less"表示求置信上限； alternative="greater"表示求置信下限。mu表示均值，它仅在假设检验中起作用，默认值为零。
在上例中如果不知道方差，就需要使用函数t.test( )来求置信区间，下面看一下在R中是如何实现的。
R程序如下： 

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

> t.test(x)

运行结果如下： 

One Sample t-test

data:  x

t = 280, df = 9, p-value <2e-16

alternative hypothesis: true mean is not equal to 0

95 percent confidence interval:

173.3 176.1

sample estimates:

mean of x 

174.7 

可以看到置信水平为0.95的置信区间为(173.3076, 176.0924)。
由于只需要置信区间的结果，因此R程序如下：

> t.test(x)$conf.int

提取出置信区间的部分, 结果如下： 

［1］ 173.3 176.1

attr(,"conf.level")

［1］ 0.95

例32下面是某市某类12个微型企业的销售收入数据(单位： 元)把下面的数据保存到文件income.txt中。
283192.80232700.0051000.00181827.0016281.17449066.00
152450.00673669.50200275.20315000.00293515.10331624.00
现要根据这个样本对这类企业销售收入的总体均值做出区间估计，置信度取95%。
假定该总体服从正态分布，根据样本数据，利用统计软件得到置信区间。
求置信区间的R程序如下： 

One Sample t-test

data:  x

t = 5.2114, df = 11, p-value = 0.0002893

alternative hypothesis: true mean is not equal to 0

95 percent confidence interval:

153108.1 376992.0

sample estimates:

mean of x 

265050.1 

下面绘制自由度为n-1的t分布密度函数图。
R软件读取数据的语句如下：

readLines("income.txt")

#读取income.txt为数据集 

y=scan("income.txt")

t. test(y,con=0.95)$con


如图32所示为相应于例32的自由度为n-1的t分布密度函数图。


图32t分布密度函数图



图32中标明了t0.025和-t0.025及中间阴影部分(1-α=0.95)的面积。
图32是由下面的R语句得到的： 

x=c(seq(-4,4,length=1000))

r1=qt(0.025,11);

r2=r1

x2=c(r1,r1,x［x<r2&x>r1］,r2,r2)

y2=c(0,dt(c(r1,x［x<r2&x>r1］,r2),11),0)

y=dt(x,11)

plot(x,y,ylab = "f(x)");

abline(0,0)

polygon(x2,y2,col="yellow")

title("Density of t(11)")

text(locator(1),expression(t［0.025］==2.200985))

#在图表上的任务位置单击，文字就出来了

text(locator(1),expression(-t［0.025］))

text(locator(1),expression(1-alpha==0.95))

3.3.2总体方差的区间估计
在R中也没有直接求σ2的置信区间的函数，需要编写自己需要的函数，下面的函数chisq.var.test()可用来求σ2的置信区间，具体实现代码如代码清单33所示。

代码清单33

chisq.var.test <- function (x,var,alpha,alternative="two.sided"){

options(digits=4)

result<-list( )

n<-length(x)

v<-var(x)

result$var<-v

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

result$chi2<-chi2

p<-pchisq(chi2,n-1)

if(alternative == "less"|alternative=="greater"){

result$p.value<-p

} else if (alternative=="two.sided") {

if(p>.5)

p<-1-p

p<-2*p

result$p.value<-p

} else return("your input is wrong")

result$conf.int<-c(

(n-1)*v/qchisq(alpha/2, df=n-1, lower.tail=F),

(n-1)*v/qchisq(alpha/2, df=n-1, lower.tail=T))

result

}

将此函数用到例31，对应的R语句如下：

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

chisq.var.test(x,1.5,0.05)

运行显示σ2的0.95置信区间为［1.793,12.628］。
运行结果如下： 

$var

［1］ 3.789

$chi2

［1］ 22.73

$p.value

［1］ 0.01365

$conf.int

［1］  1.793 12.628

3.3.3总体比例的区间估计
在R中，可直利用函数prop.test()对p进行估计与检验，其调用格式如下： 

prop.test(x, n, p = NULL,

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

conf.level = 0.95, correct = TRUE)

说明： x为样本中具有某种特性的样本数量； n为样本容量； correct选项为是否做连续性校正。根据抽样理论，p的1-α的近似置信区间为

p±zα/2(1-f)p(1-p)n-1-12n

其中，f为抽样比。由于假设样本容量很大，因此修正后p的置信度为1-α的置信区间近似地为

p±zα/2p(1-p)n-12n

它与刚才用中心极限定理推出的结论相比，区间长了1n，这是由于用连续分布去近似离散分布(超几何分布)引起的。
例33从一份共有3042人的人名录中随机抽200人，发现38人的地址已变动，试以95%的置信度估计这份名录中需要修改地址的比例。
解： 在R中输入如下语句：

> prop.test(38,200,correct=TRUE)

得到如下的结果:

1-sample proportions test with continuity correction

data:  38 out of 200, null probability 0.5

X-squared = 75.645, df = 1, p-value < 2.2e-16

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

0.1394851 0.2527281

sample estimates:

p 

0.19 

所以，以95%的置信水平认为这份名录中需要修改地址的比例p落在(0.1395, 0.2527)中, 其点估计为0.19。
如果不进行校正，相应的R语句如下：

> prop.test(38,200,correct=FALSE)

输出的结果如下： 

1-sample proportions test without continuity correction

data:  38 out of 200, null probability 0.5

X-squared = 76.88, df = 1, p-value < 2.2e-16

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

0.1416717 0.2500124

sample estimates:

p 

0.19 

此时，p的95%置信区间为(0.1417，0.2500)，其长度比修正的缩短了。
3.3.4两个总体均值之差的区间估计
在R语言中可以编写函数求置信区间，函数two.sample.ci()的定义如下： 

two.sample.ci<-function(x,y,conf.level=0.95, sigma1,sigma2 ){

options(digits=4)

m= length(x); n = length(y)

xbar=mean(x)-mean(y) 

alpha = 1 - conf.level

zstar= qnorm(1-alpha/2)* (sigma1/m+sigma2/n)^(1/2)

#函数qnorm()的返回值是给定概率p后的下分位点

xbar +c(-zstar, +zstar)

}

在R中各种概率函数都有统一的形式，即一套统一的前缀+分布函数名： 
d——表示密度函数； p——表示累积分布函数； q——表示分位数函数，能够返回特定分布的分位数； r——表示随机函数，生成特定分布的随机数。
一共提供了4类有关统计分布的函数(密度函数、累计分布函数、分位函数、随机数函数)。
(1) 正态分布的函数是norm，使用命令dnorm(0)即可获得正态分布的密度函数在0处的值为0.3989(默认为标准正态分布)。
(2) 同理，pnorm(0)是0.5，就是正态分布的累计密度函数在0处的值。
(3) 而qnorm(0.5)则得到的是0，即标准正态分布在0.5处的分位数是0(比较常用的qnorm(0.975)，就是估计中经常用到的1.96)。
(4) 最后一个rnorm(n)，则是按正态分布随机产生n个数据。
例34为比较两个小麦品种的产量，选择18块条件相似的试验田，采用相同的耕作方法做实验，结果播种甲品种的8块实验田的单位面积产量和播种乙品种的10块实验田的单位面积产量分别为甲品种，628、583、510、554、612、523、530、615； 乙品种，535、433、398、470、567、480、498、560、503、426。
假定每个品种的单位面积产量均服从正态分布，甲品种产量的方差为2140，乙品种产量的方差为3250，试求这两个品种平均面积产量差的置信区间(取α=0.05)。
用R语言进行两个总体均值的区间估计直接利用上面编写的函数two.sample.ci( )：

x<-c(628,583,510,554,612,523,530,615)

y<-c(535,433,398,470,567,480,498,560,503,426)

sigma1<-2140

sigma2<-3250

two.sample.ci(x,y,conf.level=0.95, sigma1,sigma2)

结果显示如下：

［1］34.67130.08

所以，这两个品种平均面积产量差的置信水平0.95的置信区间为(34.67,130.08)。
小样本估计的R语言如下：

x<-c(28.3,30.1,29.0,37.6,32.1,28.8,36.0,37.2,38.5,34.4,28.0,30.0)

y<-c(27.6,22.2,31.0,33.8,20.0,30.2,31.7,26.0,32.0,31.2,33.4,26.5)

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

#如果不知道两种品种产量的方差但已知两者相等，此时需在t.test( )中指定#选项var.equal=TRUE。

运行结果如下：

Two Sample t-test

data:  x and y

t = 2.2, df = 22, p-value = 0.04

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

0.14037.2597

sample estimates:

mean of x mean of y 

32.528.8 

例35
现有两种方法可以组装产品。假定第一种方法随机安排12个工人，第二种方法随机安排8个工人，即n1=12、n2=8，组装产品所需时间的具体数据如表31所示。


表31两种方法组装产品所需的时间

单位： min


方法1方法2方法1方法2

28.327.636.031.7
30.122.237.226.5
29.031.038.5—
37.633.834.4—
32.120.028.0—
28.830.230.0—
R程序如下：

x<-c(28.3,30.1,29.0,37.6,32.1,28.8,36.0,37.2,38.5,34.4,28.0,30.0)

y<-c(27.6,22.2,31.0,33.8,20.0,30.2,31.7,26.5)

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

运行结果如下：

Welch Two Sample t-test

data:  x and y

t = 2.3, df = 13, p-value = 0.04

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

0.19899.0511

sample estimates:

mean of x mean of y 

32.5027.88

例36甲、 乙两台机床分别加工某种轴承，轴承的直径分别服从正态分布N(μ1,σ21)和N(μ2,σ22)，从各自加工的轴承中分别抽取若干个轴承测其直径，结果如表32所示。试求两台机床加工的轴承直径的方差比σ21/σ22的0.95置信区间。


表32机床加工的轴的直径数据


总体样本量直径/mm

X(机床甲)820.519.819.720.420.120.019.019.9

Y(机床乙)720.719.819.520.820.419.620.2

R程序如下：

> x<-c(20.5,19.8,19.7,20.4,20.1,20.0,19.0,19.9)

> y<-c(20.7,19.8,19.5,20.8,20.4,19.6,20.2)

> var.test(x,y)

运行结果如下：

F test to compare two variances

data:  x and y

F = 0.79, num df = 7, denom df = 6, p-value = 0.8

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

0.1393 4.0600

sample estimates:

ratio of variances 

0.7932 

可见，两台机床加工的轴承的直径的方差比σ21/σ22 的0.95置信区间为(0.1393,4.0600)。 结果中sample estimates给出的是方差比σ21/σ22 的矩估计值0.7932。
3.3.5两个总体比例之差的区间估计
例37据一项市场调查，在A地区被调查的1000人中有478人喜欢品牌K，在B地区被调查的750人中有246人喜欢品牌K，试估计两地区人们喜欢品牌K比例差的95%置信区间。
解： 可以利用R中的内置函数prop.test( )求两总体的比例差的置信区间，在R中运行如下代码。

> like<-c(478, 246)

> people<-c(1000, 750)

> prop.test(like, people)

得到的结果如下：

2-sample test for equality of proportions with continuity correction

data:  like out of people

X-squared = 39, df = 1, p-value = 4e-10

alternative hypothesis: two.sided

95 percent confidence interval:

0.1031 0.1969

sample estimates:

prop 1 prop 2 

0.478  0.328 

可以看出A地区喜欢品牌K的人更多，且A、B两地区喜欢品牌K的比例之差的95%的置信区间为(0.1031, 0.1969)。
注意： 例37的结果实际上是经过连续性修改后得到的。
根据prop.test()函数也可以自己编写没有修正的两比例之间的区间估计函数ratio.ci()：

ratio.ci<-function(x, y, n1, n2, conf.level=0.95){

xbar1=x/n1;

xbar2=y/n2;

xbar=xbar1-xbar2;

alpha=1-conf.level;

zstar=qnorm(1-alpha/2) *(xbar1*(1-xbar1)/n1+xbar2*(1-xbar2)/n2)^(1/2);

xbar +c(-zstar, +zstar)

}

将ratio.ci()函数用到例37中，运行如下代码：

ratio.ci(478,246,1000,750,conf.level=0.95)

得到结果如下：

［1］ 0.1043 0.1957

这时，两比例之差的95%的置信区间为(0.1043，0.1957)，长度修正的结果略小了些。
3.4估计总体均值时样本量的确定
在R中可以定义如下函数size.norm1()求样本容量。

size.norm1<-function(d,var,conf.level) {

alpha = 1 - conf.level

((qnorm(1-alpha/2)*var^(1/2))/d)^2

}

在R中可以通过循环确定样本容量，size.norm2()的定义如下：

size.norm2<-function(s,alpha,d,m){

t0<-qt(alpha/2,m,lower.tail=FALSE)

n0<-(t0*s/d)^2

t1<-qt(alpha/2,n0,lower.tail=FALSE)

n1<-(t1*s/d)^2

while(abs(n1-n0)>0.5){

n0<-(qt(alpha/2,n1,lower.tail=FALSE)*s/d)^2

n1<-(qt(alpha/2,n0,lower.tail=FALSE)*s/d)^2

}

n1

}

说明： m是事先给定的一个很大的数。
例38某公司生产了一批新产品，产品总体服从正态分布，现要估计这批产品的平均质量，最大允许误差为2，样本标准差S=10，试问α=0.01下要抽取多少样本。
解： 在R中运行如下程序。

> size.norm2(10,0.01,2,100)

运行结果如下：

［1］ 169.7


也就是说，在最大允许误差为2时，应抽取170个样本。
例39抽查用克矽平质量的矽肺病患者10人，得到他们治疗前后的血红蛋白差(单位： g)如下： 2.7，-1.2，-1.0，0，0.7，2.0，3.7，-0.6，0.8，-0.3。现要检验治疗前后血红蛋白差是否服从正态分布(α=0.05)。
解： 把观测值按非降序排列成-1.2，-1.0，-0.6，-0.3，0，0.7，0.8，2.0，2.7，3.7。
n=10，α=0.05，W1-0.05=0.842，所以拒绝域为{W≤0.842}。
将排列后的数据填表，其中ai(W)这一列的值是由GB/T 4882—2001(《数据的统计处理和解释正态性检验》)根据n的值查得的。为了计算统计量，列出如表33所示的计算表。


表33矽肺病患者的相关数据


ix(i) x(n+1-i)x(n+1-i)-x(i)ai(W)

1-1.2
3.7
4.9
0.5739
2
-1.0
2.7
3.7
0.3291
3
-0.6
2.0
2.6
0.2141
4
-0.3
0.8
1.1
0.1224
5
0
0.7
0.7
0.0399

其中，x(i)为观察值的前一半按照升序排列，x(n+1-i)为大的一半按照降序排列，计算出： 

∑10i=1x(i)=6.8, ∑10i=1x2(i)=29, ∑10i=1x(i)-2=24.376,
∑5i=1ai(W)x(n+1-i)-x(i)=4.74901，
W=4.74901224.376=0.9252

对于显著性水平α=0.05，查表知，n=10时，p0.05=0.842。由于0.9252＞0.842，所以不拒绝正态性的原假设。
3.5R语言中的ShapiroWilk检验
shapiro.test(x)函数只有一个参数，即数据集x。x可以是数值型向量，但是非丢失数据需要在3~5000范围内。
例310现有11个随机抽取的样本的体重数据为148、154、158、160、161、162、166、170、182、195、236。
编写R语言如下：

> k<-c(148 ,154, 158, 160, 161, 162, 166, 170, 182, 195, 236)

> shapiro.test(k)

Shapiro-Wilk normality test

data:k

W = 0.7888, p-value = 0.006704

结果中W统计量为0.7888，非常接近于0，但是其P值小于0.05，因此可以拒绝原假设，即该数据不符合正态分布。因此，在这个例子中也可以发现较大的W统计量的情况下，依然可能拒绝其符合正态分布。

>shapiro.test(rnorm(100, mean = 5, sd = 3))

Shapiro-Wilk normality test

data: rnorm(100, mean=5, sd=3)

W=0.9926, p-value=0.863

>shapiro.test(runif(100, min=2, max=4))

Shapiro-Wilk normality test

data: runif(100, min=2, max=4)

W=0.9561, p-value=0.00214

在这个例子中，第一个命令是检验一个随机生产的100个数据，该数据集符合均值为5、标准差为3的正态分布，其W统计量接近1，P值显著大于0.05，不能拒绝其符合正态分布。第二个命令检验一个随机产生的100个数据，但是该数据是符合最小值为2、最大值为4的均匀分布，其W统计量也是接近1，但是其P值小于0.05，所以有足够理由拒绝其符合正态分布的说法。
本章小结
本章都是关于参数估计的内容。首先，介绍了参数估计的原理，总体方差、总体比例的区间估计，统计量的分布(χ2分布、t分布和F分布)； 其次，介绍了如何运用参数估计的区间估计进行R编程计算； 最后，介绍了ShapiroWilk检验的R语言编程。
思考与练习
1. 为了解居民用于服装消费的支出情况，随机抽取90户居民组成一个简单的随机样本，计算得到样本均值为810元，样本标准差为85元，试建立该地区每户居民平均用于服装消费支出的95%的置信区间。
2. 某市为了解居民住房情况，抽查了n=2000户家庭，其中人均不足5m2困难户有x=214个，通过样本信息计算该市困难户比率P的置信区间(置信度为0.95)。
3. 从一份共有3062人的人名名录中随机抽取300人，发现42人的地址已变动，试以95%的置信度估计这份名录中需要修改地址的比例。