非参数统计:基于R语言案例分析
上QQ阅读APP看本书,新人免费读10天
设备和账号都新为新人

1.2 估计

统计量的一个基本目的是估计总体的未知性质。这些估计出来的未知性质通常是用数字表示的,并且包括可列举的一些项目,例如未知比率、均值、概率等。事实上,估计是基于样本的(如果有概率描述,则是随机样本),并且估计是关于随机变量分布未知性质的有根据推测,这里随机变量表示对总体研究感兴趣的量。例如,我们可以用产品中样本的不合格率来估计总体的不合格率。用来做估计的统计量自然叫做估计量(estimator)。本节我们将要讨论一些估计量,例如样本均值(sample mean)、样本方差(sample variance)和样本分位数(sample quantiles)。我们首先引入一个与众不同的估计量——经验分布函数(empirical distribution function)。

1.2.1 经验分布函数

一个随机变量的真实分布函数一般是未知的,有时我们只能够推测分布函数的形式,或将推测作为真实分布函数的一个近似。根据样本的观测值作经验分布函数图,以此来作为整个未知分布函数Fx)的估计,这是推测分布函数的一种好方法。下面将用具体的例子来介绍这种估计方法,由此我们给出定义:

【定义1.4】设X1X2,…,Xn是一组随机样本,经验分布函数Sx)(简称为EDF)是x的函数,它在x点的取值为小于或等于xXi在样本总数中所占的比例,其中-∞<x<∞。

【例1.3】在一项体能研究中,一高中随机抽取了5名男生,记录他们跑完1英里的时间。时间(转化成分钟后)分别为6.23,5.58,7.06,6.42,5.20。由于经验分布函数Sx)是小于或等于xXi在样本总数中所占的比例,根据这组特定样本有如下经验分布函数:

我们也可以很方便地画出该经验分布函数的图象,并且从例1.3中可以看出,经验分布函数总是阶梯函数,每阶的高度是1/n,并且只在样本取值处有变化。我们从左到右来考虑经验分布函数的值,注意到Sx)在样本最小值前均取值为零,在每个样本取值处会增加一阶的跃度,每个跃度是1/n。在样本最大之处Sx)取最大值1,并且在剩下所有比样本最大值大的x处都取1。Sx)很像非降的取值从0到1的分布函数。但Sx)只是由经验(来自样本)确定的,并由此而得名。

例1.3只描述了Sx)的一组观测值,其他的样本值将产生另外不同的经验分布函数Sx),对应的图象当然也不同。这表明了Sx)的随机性,从这个意义上讲,它是一个随机变量。但是,由于它是一个函数,且观测值是整个图象而不是单个值,所以称Sx)为随机函数(random function)更加合适。因为它能够相当好地估计随机变量的分布函数,所以它通常用作一个估计量。为了区分经验(或样本)分布函数,我们称随机变量的分布函数为总体分布函数。

从某种意义上讲,经验分布函数的观测值可以认为是总体分布函数的取值,准确地讲,基于样本观测值x1x2,…,xnSx)的一个观测值和取x1x2,…,xn中每个值的概率都是1/n的随机变量的分布是一样的。这种随机变量的分布函数是一个阶梯函数,且在每个数值x1x2,…,xn处的跃度为1/n,因此,我们很容易得到随机变量的均值、方差和分位数。

1.2.2 估计量

为了区分真实“总体”的均值、方差和分位数,由样本计算得到的均值、方差和分位数分别称为样本均值、样本方差和样本分位数。正如经验分布函数可以作为总体分布函数的估计量,样本均值、方差、分位数也可以分别作为总体均值、方差、分位数的估计量。

【定义1.5】设X1X2,…,Xn是一组随机样本,样本p分位数Qp满足以下两个条件:

(1)小于QpXi的比例小于等于p

(2)大于QpXi的比例小于等于1-p

正如总体分位数从总体分布函数得到的方式一样,每个样本分位数都可以由经验分布函数得到。样本p分位数是Sx)=p处的x值,如果有不止一个x值满足Sx)=p,我们取最大值与最小值的均值作为该样本p分位数,与总体分位数的处理一样。样本p分位数Qp取决于随机变量的取值,因此它是一个统计量。注意,为简便起见,我们只针对随机样本定义样本分位数。

一种直接由样本而不通过Sx)的图象而得到样本p分位数的方法是,用p乘以样本容量n,四舍五入得到一个相邻的较大整数,那么样本分位数是以(p·n)和(p·n+1)为秩的两个次序统计量的均值。

【例1.4】从市区妇女俱乐部中的已婚妇女中随机抽取6名妇女,记录每位妇女的孩子个数,分别为0,2,1,2,3,4,经验分布函数如图1-1所示。样本中位数Q0.5是2,样本分位数Q0.25Q0.75分别是1和3。按照我们的约定,1/3样本分位数Q1/3是1和2的平均值,即1.5。这些数可以用来估计未知的总体分位数。

图1-1 经验分布函数

注意此时样本均值和方差计算起来比较简单。对于这种简便的计算方法,我们给出下面的定义。

【定义1.6】设X1X2,…,Xn是一组随机样本,则样本均值定义如下:

S 2定义如下:

同时,它还等价于下式

【例1.5】例1.4中随机样本0,2,1,2,3,4的样本均值是

因此,未知均值的估计是2,未知方差的估计是1

除了经验分布函数外,上面所介绍的估计量还提供了未知总体参数的一种叫做点估计(point estimate)的方法,也就是前面的例子中我们得到的未知均值的估计方法。

我们经常更喜欢,同时也更谨慎地说:“我们以95%的置信水平认为未知均值落在1.3与2.7之间。”这种估计称为区间估计(interval estimate)。区间估计量由两个统计量组成,它们是区间的两个端点。置信系数(confidence coefficient)是区间估计量包含未知总体参数的概率。前面的叙述中置信系数是0.95。区间和置信系数一起称为未知量的置信区间。

点估计是比较容易的,因为点估计只需要考虑一个数,任意一个数。但是,有些点估计量比其余的估计量要好。为了比较哪种估计量更好,其点估计的比较标准几乎在任何一本概率统计导论的书中都可以找到。

一个好的估计量的标准之一是无偏性(unbiased)。在下面的讨论中,我们通常用希腊字母θμσρ表示参数,用表示θ的统计量。

【定义1.7】如果,则称统计量是总体参数θ的无偏估计(unbiased esti mator)。

下面的定理表明是总体均值的无偏估计。

【定理1.1】X1X2,…,Xn是一组来自均值为μ,方差是σ2总体的独立随机变量,那么

证明 由于X1X2,…,Xn是独立的,所以有下式成立

这表明是总体均值的无偏估计,同时有

经过简单的代数运算,可以得到下式

即完成了定理的证明。

1.2.3 标准误差

一个估计量的标准差通常称为标准误差(standard error),所以它不会和总体标准差

的概念混淆,因为这是一个完全不同的概念,如定理1.1所示,的标准误差是

1.2.4 无偏估计量s2

我们可以证明得到S2不是σ2的无偏估计,因此,习惯上我们使用无偏估计s2作为σ2的估计量。

但值得注意的是Ss都不是总体标准差σ的无偏估计。

1.2.5 渐近置信区间

由定理1.1及中心极限定理可知,如果X1X2,…,Xn是一组均值为μ,方差是σ2的独立随机变量,那么,当n趋于无穷时,

的分布函数趋于正态标准分布。在实际应用中,如果n足够大,

的概率近似是1,其中z1-α/2表示标准正态随机变量的(1/2)分位数。

上面的不等式经过代数整理,可写为

n很大时,上式给出了μ的近似置信区间。进一步讲,既然σ很少已知,那么通常情况下,当n足够大,对于来自非零方差有限总体的随机样本,根据中心极限定理我们可以用s来估计σ,从而得到μ的近似置信区间。对大多数的实际问题来说,当样本容量超过30就可以被认为是“足够大”了。

【例1.6】一窝猪的数量越多,那么农场主可获得的利润就越多。国家实验中心正在研究一种能够提高每窝猪产量的新技术,记录的55窝猪中平均每窝存活猪的数量是9.8,且s=1.4。平均每窝存活猪的总体均值95%的近似置信区间下限为

近似置信区间上限为

从而,可以以95%的置信水平认为真实的每窝猪数量的均值在9.43与10.17之间。

1.2.6 自助法

定理1.1给出了的均值和方差,S2的均值也不难得到。但是,S2的方差却不是那么容易计算。很多统计量作为总体参数的估计量是很难像定理1.1那样从理论上来推导的,甚至是不可能的,所以我们用其他的方法来估计它们的均值和方差。其中一种方法称为自助法(bootstrap)。

自助法是从原始容量为n的随机样本的观测值中有放回地抽取n个值,也就是说,一些原始随机样本的观测值在“自助样本”中可以出现一次,也可以出现多次或者根本不出现。自助样本的个数总是和原始随机样本中观测值的个数相等。

每组自助样本都可以计算出我们所关心的估计量。通过计算机模拟,我们可以由原始随机样本的观测值得到成百上千组自助样本,每组自助样本都会产生一个值。这成百上千个估计值的样本均值、样本标准差(sS)可以用来估计的总体均值和总体标准差(标准误差)。事实上,在自助法中,这些估计值的经验分布函数可以作为的真实的总体分布函数的一个估计。显然,自助法每一步都依赖于原始随机样本值,不同样本值的集合会产生不同估计值的集合。

对于简单估计一个估计量的均值和标准差,自助法的重复次数很少超过100或200,25次左右已经足够。但是,要得到置信区间则需要做大量的重复试验。得到θ近似置信区间的一种方法是利用自助样本估计量α/2和1/2的样本分位数。Efrom和Tib-shirane(1986)建议至少要做250次自助重复试验,他们同时也给出了另外一种更加精确获得置信区间的方法,这种方法需要更多的自助重复试验,至少要做1 000次。

【例1.7】在例1.6中,对X用中心极限定理,得到了每窝猪数量总体均值的置信水平为95%的近似置信区间。现在我们用自助法,求出每窝猪数量总体标准差σ的置信水平为95%的近似置信区间,这个参数对检验是非常有用的,因为每窝猪数量相差不大的情况(σ较小)要比一些窝数量很少而其他窝数量很多的情况(σ较大)好得多。

原始样本55个观测值是从1到55的编号。

观测号 1 2 3 4…55

每窝数量 9 9 8 6…11(55个) s=1.4

现在从1到55号进行有放回的抽样,得到55 个数,得到第一次自助样本,由它计算出估计量s

自助样本#1:

观测号 4 17 4 28…16

每窝数量 6 9 6 10…9(55个)=1.6

这个过程重复250次(自助法的这个过程可以重复所需要的尽可能多的次数,但是建议求置信区间至少需要250次)。

自助样本#2:

观测号 4 17 4 28……16

每窝数量 6 9 6 10……9(55个)=1.8

将这个过程继续重复下去,直到

自助样本#250:

观测号 6 1 55 14…17

每窝数量 10 9 11 11…9(55个)=1.1

s*的样本0.025分位数可以得到95%置信区间的近似置信下限,因为0.025(250)=6.25,四舍五入为6,所以置信区间的下限是第6个次序统计量s*(6)。由样本0.975分位数,即第244个次序统计量(0.975×250=243.75,四舍五入到244)得到,s*(244)即为置信上限。在这个问题中,我们需要将s*的250个值从小到大排序,即0.7,0.8,0.8,0.9,0.9,0.9,1.0,…,2.0,2.0,2.2,2.3,2.3,2.4,2.7,从而,95%置信区间是从1.0到2.0。通过计算250个s*的标准差S而得到s的标准误差的估计,正如由250个s*的均值可以算出s均值的估计值一样,我们还可以算出其他所关心的统计量。

几乎所有的计算机中的统计软件包,甚至许多便宜的手动计算器都可以计算出我们前面所讨论的点估计。但是,自助法可不太容易得到。R、S-Plus、SYSTAT、Resampling Stats和Stata等统计软件中都可以找到自助法的计算程序。

1.2.7 一般参数估计

可以用两个问题来概括估计一个未知参数θ的大致过程:

(1)使用哪个统计量?建议仿照我们前面所举的例子中计算以及分位数估计量的过程,看参数θ在总体分布函数Fx)上是如何定义的,相应地就可以从经验分布函数来定义参数θ的估计

(2)该统计量是否优良?这通常用一个估计量的标准差(称为标准误差)来衡量。定理1.1给出了的标准误差是。其他统计量的标准误差并不像均值这样容易计算,但可以通过自助法来估计,自助法还可以给出的整个分布函数的估计以及θ的近似置信区间。

1.2.8 生存函数

经验分布函数Sx)将随机样本X1X2,…,Xn与总体分布函数Fx)有机地结合起来,因为它用小于等于x观测值的频率来估计PXx)=Fx)。在寿命测试、医疗后续工作,以及其他领域中,和分布函数同样有用的是生命函数(survival function),即Px)=1-Fx),这里所关心的变量为事物的寿命(从发生到结束持续的时间),可以是人、动物或无生命产品的寿命,也可以简单地从某个起点开始,直到某事件发生,如痊愈、到达、离开等所持续的时间。

Px)的一个自然估计是它的经验生存函数,

它是样本X1X2,…,Xn中超过x值的频率。

1.2.9 Kaplan-Meier估计

Kaplan与Meier(1985)用x表示死亡时间(应当注意,由于试验中某些项的缺失,死亡时间在某些情况下是不可观测的,比如试验中研究对象中途离开、进入试验较晚,或者试验结束后才死亡等),提供了一种从缺失数据中获取信息的方法,即在缺失前死亡还没有发生。他们利用了下面的事实:如果死亡发生在x后,那么在x前的任意时刻后死亡仍然会发生。下面是他们的理由。

由条件概率的定义,我们可以得到,对于x0x1,有

假设在第一年年初有100个研究对象参加测试,第一年年底只有30个存活,我们用下式估计P(1),

这里X表示研究对象个体的寿命。

接着,假设第二年年初又有另外1 000个个体参加试验,第二年年底,1 000个中有250个存活,而最初100个个体中存活的30个只剩下10个。我们可以用最初的100个个体来估计P(2),

此时,我们可以用第二年新参加的1 000个个体的信息来更新估计P(1),因为两年中参加试验的个体共有1 100个,其中共有250+30=280个存活,改进后P(1)的估计为

进一步讲,由1 100个个体中有280个存活和等式(1-13),我们可以用改进后的P(1)来改进P(2),得到估计式

改进后的估计^ Px>1)为0.255。不幸的是,我们无法改进PX>2│X>1)的估计值,因为我们不知道在接下来一年的试验中1 000个观测个体有多少个存活。所以我们用下面的估计量

因为它仅用到了已知信息,即第一年年底有30个存活,第二年年底有10个存活,那么P(2)的改进估计就是

Kaplan与Meier把上面的方法推广应用到了一般情况,设u1u2<…<uk表示k个个体“寿命”,这里寿命是指从开始到死亡,或者到研究缺失所持续的时间,并令

pi=P(XuiXui-1)

用下式估计

在时刻ui缺失的个体,可认为在时刻ui以后仍然存活,因为知道他们在时刻ui还是活着的。而在时刻ui死亡的个体不可能在时刻ui后是存活的。在第一次死亡或缺失的计算时,的分母是参加试验个体的总数。

PX)的Kaplan-Meier估计为

其中,乘积是对所有寿命uixi进行的。注意,这个估计量是一个递降的阶梯函数,且只在观测的死亡时间取值发生变化,而且由,这种方法可定义更一般的经验分布函数Sx),用于针对缺失数据的研究。

R,Minitab,S-Plus和SYSTAT软件都可以制作生存曲线,特别是Kaplan-Meier估计。在估计生存曲线时,有时我们需要求出缺失数据的方差,在本书中没有介绍,但这里提到的软件包可以解决这些问题。

【例1.8】要测试10个汽车上风扇皮带的质量,我们把它们装到车上,并记录每辆车上皮带所能承受的里程数。测试结束后,5个皮带都断裂了。寿命(以千英里计)分别为77,47,81,56,80。另外5个没有断裂,寿命(以千英里计)分别是62,60,43,71,37。生存函数的Kaplan-Meier估计如表1-1所示。

表1-1

对所有的x>0,的图象如图1-2所示。

图1-2 生存函数Px)的Kaplan-Meier估计

如果没有缺失数据,只有死亡数据,那么Kaplan-Meier估计和1-Sx)是一样的,从开始,在每个死亡时刻以1/n为阶梯高度下降,直到。如果既有缺失数据又有死亡数据,那么从1.0开始,下降的阶梯高度不再一致。如果在最后知道死亡时刻之后仍有缺失数据,那么^Px)不会下降到0,并且对于那些在最后已知缺失数据后面的x,它没有定义。在这种情况下,通过本章前面介绍的常用方法,用Sx)来估计Fx)的某些相关参数,如均值、方差等是不太合适的,但是可以估计它的一些分位数。Kaplan和Meier(1958)也介绍了一些特殊的方法。

在这种意义下,点估计是一种非参数统计方法,因为不用了解任何关于未知分布函数的形式就可以做出点估计。本章中的例子足以说明这一点。

很难说清构造置信区间的方法是参数的还是非参数的。如果构造置信区间时不需要任何分布函数的形式,这显然是非参数方法,例1.7和1.8所示的近似方法就是非参数方法。另外,如果方法要求未知分布函数是正态分布,或是其他的特殊形式,那么这种方法就是参数的。我们将在接下来的章节中介绍其他几种构造置信区间的非参数方法。