第5章 描述性统计.doc

第5章 描述性统计采集到大量的样本数据以后，常常需要用一些统计量来描述数据的集中程度和离散程度，并通过这些指标来对数据的总体特征进行归纳。描述样本数据集中趋势的统计量有算术平均值、中位数、众数、几何均值、调和均值和截尾均值等。描述样本数据离中趋势的统计量包括极差、平均值、平均绝对差、方差和标准差等，此外还有峰度、偏差、分位数和相关系数等统计量，也能描述样本数据的某些特征。5.1 描述集中趋势的统计量5.1.1 几何均值样本数据的几何均值m可以根据下式求得Matlab中利用geomean函数计算样本的几何均值，其语法格式如下：l m = geomean(X) 函数计算样本的几何均值。若X是矢量，则geomean(X)返回值为数据X中元素的几何均值。若X为矩阵，则geomean(X)返回值为一个行矢量，包含每列数据的几何均值。若X为N维数组，geomean函数沿X的第一个非单一元素维度进行计算。l m = geomean(X,dim) 计算X的第dim维的几何均值。【例 5-1】 样本均值大于或等于样本的几何均值。x = exprnd(1,10,6);geometric = geomean(x)geometric = 0.7466 0.6061 0.6038 0.2569 0.7539 0.3478average = mean(x)average = 1.3509 1583 0.9741 0.5319 1.0088 0.81225.1.2 调和均值样本数据的调和平均值m定义为Matlab中利用harmmean函数计算样本数据的调和平均值，其语法格式如下：l m = harrmmean 计算样本的调和平均值。若X为矢量，则harmmean(X)函数返回值为X中元素的调和平均值。若X为矩阵，则harmmean(X)函数返回值为包含每列元素调和平均值的行矢量。若X为N维数组，harmmean函数沿X第一个非单一元素维度进行计算。l m = harmmean (X,dim) 计算X的第dim维的几何均值。【例 5-2】 样本均值大于或等于样本的调和平均值。x = exprnd(1,10,6);harmonic = harmmean(x)harmonic = 0.3382 0.3200 0.3710 0.0540 0.4936 0.0907average = mean(x)average = 1.3509 1583 0.9741 0.5319 1.0088 0.81225.1.3 算术平均值样本数据的算术平均值可用下式定义Matlab中利用mean函数计算矢量和矩阵中元素的均值，其语法格式如下：l m = mean(X)，若X为矢量，mean(X)返回值为X中元素的均值；若X为矩阵，mean(X)返回值为包含X中每列元素均值的行矢量；l m = mean(X,dim) 计算X的第dim维元素的均值。【例 5-3】 下面的命令行生成5个包含100个服从正态分布的随机数的样本，然后计算每个样本的算术平均值。x = normmd(0,1,100,5);xbar = mean(x)xbar =0.727 0.0264 0.0351 00424 0.0752 5.1.4 中值median函数计算矢量和矩阵中元素的中值，其函数的调用格式为：m = median(X) 计算样本数据的中值。中值是样本数据中心趋势的稳健估计，因为异常值的影响较小。对于矢量，median(X)为矢量X中元素的中值。对于矩阵，median(X)为包含每一列中元素中值的行矢量。计算中值需要首先进行排序，因此计算大型矩阵的中值矢量时比较费时。【例 5-4】xodd = 1:5;modd = median(xodd)modd =3 meven = median(xeven)meven = 2.5000下列演示中值对于异常值的稳健性。xoutlier = x,10000;moutlier = median(xoutlier)moutlier = 35.1.5 截尾均值对样本数据进行排序以后，去掉两端的部分极值，然后对剩下的数据求算术平均值，得到截尾均值。Matlab利用trimmean函数计算截尾均值，其语法格式如下：l m = trimmean(X,percent) 若X为矢量，则m为X中元素的截尾均值，即剔除测量值中最大和最小的k个数据以后，计算样本X的均值，k=n*(percent/100/2)，n为X的维数；若X为矩阵，则m为行矢量，其元素为X中各列元素的截尾均值；若X为多维数组，则m沿X中的第一个非单一元素维度进行计算。Percent为0和100之间的数。l m = trimmean(X,percent,dim) 沿X的第dim维计算截尾均值。l m = trimmean(X,percent,flag)处理当k不是整数时如何截尾，flag取值及含义如下：round缺省值，对k就近取整。floor对k向下取整weight若k=i+f，其中i是整数部分，f是小数部分，以 (1-f) 为第(i+1)个和第（n-i）个元素权值，求两个元素之间所有元素的加权均值。l m = trimmean(X,percent,flag,dim)计算X的第dim维元素的截尾平均。截尾均值为样本位置参数的稳健性估计。若数据中有异常值，截尾均值为数据中心的一个更具代表性的估计。若所有数据取自服从同一分布的总体，则使用样本均值比使用截尾均值更有效。【例 5-5】下面用蒙特卡洛法模拟正态数据的10%截尾均值相对于样本均值的有效性。值小于1。说明正态条件下截尾均值不如算术平均值有效。rng(default);x = normmd(0,1,100,100);m = mean(x);trim = trimmean(x,10);sm = std(m);strim = std(trim);efficiency = (sm/strim).2输出为：efficiency = 0.96635.2 描述离中趋势的统计量描述离中趋势的统计量包括四分位差、均值绝对差、极差、方差和标准差等。5.2.1 四分位差四分位差指的是样本数据从小到大排序后75%与25%位置处的值之差。IQR是数据极差的稳健性估计。因为上下25%的数据变化对其没有影响。Matlab中用iqr函数计算样本的四分位差(IQR)，其语法格式如下：l y = iqr(X) 计算X的四分位差。若X为矢量，则y为X的75%与25%位置处的值之差；若X为矩阵，则y是行向量，其元素为对应X的各列的四分位差；若X为多维数组，iqr函数沿X的第一个非单一元素维进行计算。l m = iqr(X,dim) 计算X的第dim维元素的四分位差。若数据中没有异常值，则IQR用于衡量数据的极差比标准差更具代表性。当数据取自正态分布总体时，标准差比IQR有效。常用IQR*0.7413来代替标准差。【例 5-6】下面用蒙特卡洛法模拟正态数据的IQR相对于样本标准差的有效性。结果仅为0.33，说明正太条件下IQR不如标准差有效。x = normmd(0,1,100,100);s = std(x)s_IQR = 0.7413*iqr(x);efficiency = (norm(s-1)./norm(s_IQR - 1).2输出为：efficiency = 0.32975.2.2 均值绝对差利用mad函数可以计算数据样本的均值或中值绝对差(MAD).l y = mad(X) 计算X中数据的均值绝对差。若X为矢量，则y为mean(abs(X-mean(X)；若X为矩阵，则y为包含X中每列数据均值绝对差的行矢量；如果X为多维数组，则mad函数计算第一个非单一元素维的均值绝对差。l mad(X,0) 与mad(X)相同，使用均值。l mad(X,1) 基于中值计算y，即y=median(abs(X median(X)。l mad(X,flag,dim) 沿X的第dim维计算MAD。该函数将NaN视为缺失值并删除。对于正态分布数据，可以用下式估计标准差：sigma = 1.253*mad(X,0); sigma = 1.4826*mad(X,1); 【例 5-7】对正态分布数据添加离群点后，比较不同估计的鲁棒性。结果越接近1，说明正态条件下，对应度量的更有鲁棒性。x = normrnd(0,1,1,50);xo = x 10; % 添加离群点r1 = std(xo)/std(x)r1 = 1.7385r2 = mad(xo,0)/mad(x,0)r2 = 1.2306r3 = mad(xo,1)/mad(x,1)r3 = 1.06025.2.3 极差极差指的是样本中最大值与最小值之间的差值。用range函数计算样本的极差，其语法格式如下：l y = range(X) 返回极差。若X为矢量，range(X)为X中元素的极差；若X为矩阵，range(X)为行矢量，包含X中对应各列中元素的极差；若X为多维数组，range函数沿X的第一个非单一元素维进行计算。l m = range(X,dim) 计算X的第dim维元素的极差。用极差估计样本数据的范围具有计算简便的优点；缺点是异常值对它的影响较大，因此它是一个不可靠的估计值。【例 5-8】大样本标准正态分布随机数的极差近似为6。下面首先生成5个包含1000个服从正态分布的随机数的样本，然后进行求极差的运算。rv = normmd(0,1,1000,5);near6 = range(rv)输出为：near6 = 6.1451 6.4986 6.2909 5.8894 7.00025.2.4 方差用var函数计算样本的方差。其调用格式和描述如下。l y = var(X) 计算X中数据的方差。对矢量而言，var(X)为X中元素的方差。对于矩阵而言，var(X)是包含X中每一列元素方差的行矢量，通过除以n-1来达到标称化，其中n为样本大小。对于正态分布数据，这使var(x)成为的最小方差无偏估计量。l y = var(X,1) 通过处于n来标称化并生成样本数据的二级矩。l y = var(X,w) 使用权重矢量w计算方法。w中元素的个数必须等于矩阵X的行数，对于矢量X，w和X必须在长度上匹配。W的每个元素必须为正。l y = var(X,w,dim) 计算X的第dim维元素的方差。w为0时，使用默认的N-1进行标称化；w=1时，使用N进行标称化。注意：令SS为X矢量中元素与其均值之间的离差平安和，则var(X) = SS/(n-1)为的最小方差无偏估计量，var(X,1)=SS/n为的最大似然估计量。【例 5-9】x = -1 1;w = 1 3;v1 = var(x)v1 = 2v2 = var(x,1)v2 = 1v3 = var(x,w)v3 = 0.75005.2.5 标准差有两种样本数据的标准差计算公式，分别为：(1)(2)式中，样本均值为用std函数计算样本的标准差。l s = std(X) ，利用公式(1)计算X中数据样本的标准差。若X为矢量，std(X)为X中元素的标准差；若X为矩阵，std(X)为包含X中各列元素标准差的行矢量；若X是N维数组，则std(X)为X的第一个非单一元素维的标准差。l s = std(X,flag)，当flag=0时，等同于std(X)；当flag=1时，std(X,1)返回返回(2)式计算的标准差，用n对X进行标称化，结果Y为样本关于其均值二阶矩的平方根。l s = std(X,flag,dim) 计算X的第dim维元素的标准差。flag为0时，用n-1进行标称化；falg为1时，用n进行标称化。【例 5-10】下面首先生成6列服从标准正态分布的随机数，每列有100个数。每一列中，标准差y的期望值均为1。x = normmd(0,1,100,6);y = std(x)输出为：y = 0.9536 1.0628 1.0860 0.9927 0.9605 1.0254y = std(x,1)5.3 分组数据描述利用grpstats函数计算分组概括统计量。l means=grpstats(X,group)根据group参数，按照group分组的X中对应数据的各列均值。如果是单一的分组变量，means中对应每一个分组变量值给出一个行向量。Grpstats函数按照字母序排序(当分组变量是字符串数组时)或按数字升序排序(当分组变量时数值型时)。如果分组变量group是cell数组，则对应于每一个分组变量的观测唯一组合值，means中存在一行与之对应。Grpstats按照第一个分组变量对组排序，然后按照第二个分组变量对每一个具有相同第一分组变量值的数据进行排序，以此类推。如果X是矩阵，则means是具有相同列数的矩阵。Means的每一列给出对应X相应列的分组均值。l grpstats（x,group,alpha）按照分组变量值对X中的数据计算分组均值的置信区间，绘制每个均值100*（1-alpha）%置信区间的图形，横轴是分组变量值。 如果X是矩阵，则grpstats绘制X各列的均值及置信区间。如果group是cell数组，则grpstats关于grpstats的其他调用格式，请参考matlab帮助。【例5-11】Load(hospital)%载入数据，该数据具有100个观测，7个变量%创建只包含Sex,Age,Weight, andSmoker的数据数组。Sex是标称数组，取值为Male和Famale。Age和Weight是数值型变量，Smoker是逻辑型变量。ds = hospital(:,Sex,Age,Weight,Smoker);%按照Sex分组，计算数值型和逻辑型变量的均值Statsrray = grpstats (ds, Sex)输出为:statarray = SexGroupCountmean_Age mean_Weight mean_SmokerFemale Female5337.717 130.47 0.24528 MaleMale4738.915180.530.446815.4 包含缺失数据的样本描述MATLAB的统计工具箱中有一组名称以nan为前缀的函数，用于描述包含缺失数据的样本。表5-1示给了包含缺失数据的样本统计量的计算函数及其调用格式，具体调用格式的解释可以参考matlab帮助。表5-1 包含缺失数据的样本统计量函数名称功能调用格式nancov协方差,忽略缺失数据。Y = nancov(X)Y = nancov(X1,X2)Y = nancov(.,1)Y = nancov(.,pairwise)nanmax求样本数据的最大值,忽略缺失数据。y = nanmax(X)Y = nanmax(X1,X2)y = nanmax(X,dim)y,indices = nanmax(.)nanmin求样本数据的最小值,忽略缺失数据。y = nanmin(X)Y = nanmin(X1,X2)y = nanmin(X,dim)y,indices = nanmin(.)nanmean求样本数据的均值,忽略缺失数据。y = nanmean(X)y = nanmean(X,dim)nanmedian求样本数据的中位数,忽略缺失数据。y = nanmedian(X)y = nanmedian(X,dim)nanstd求样本数据的标准差,忽略缺失数据。y = nanstd(X)y = nanstd(X,1)y = nanstd(X,flag,dim)nansum求样本数据的和,忽略缺失数据。y = nansum(X)y = nansum(X,dim)nanvar求样本数据的方差,忽略缺失数据。y = nanvar(X)y = nanvar(X,1)y = nanvar(X,w)y = nanvar(X,w,dim)下面以nanmean函数为例介绍这一类函数的用法。nanmax函数的运算包含缺失数据的样本数据的均值，其调用格式如下：l m=nanmean(X)，移除NaN值，返回有效数据的均值。NaN表示缺失值，若X为矢量，nanmean(X)表示X的剩余元素的均值；若X为矩阵，nanmean(X)为X中各列有效元素的均值构成的行向量。l m=nanmean(X，dim)，返回X的第dim维数据的均值。注意：(1)NaN表示缺失值，缺失值与零不同，它表示对应的位置上没有观测值，不能简单的用零代替。(2)如果X的某一维全是NaN，则当移除NaN后，该列是空的，因此，剩余元素的和是0。由于该列均值是0做除数得到的，因此其输出结果是NaN，但是并不表示NaN值。【例5-12】求带有缺失值数据的列均值。X = magic(3);X(1 6:9) = repmat(NaN,1,5)X = NaN 1 NaN 3 5 NaN 4 NaN NaNy = nanmean(X)y =3.5000 0000 NaN5.5 百分位数和图形描述Prctile函数计算样本的百分位数，该函数的调用格式如下。l Y=prctile(X，p) 计算大X中p的值，p的值必须介于0至100之间，对于矢量而言，prctile(X，p)为X中元素的p百分位数。若p=50，则Y为X的中值，对于矩阵X和标量p，prctile(X，p)为包含每一列的p百分位数的行向量，若p为矢量。则Y的第i行为X的p(i)。对于多维数组，prctile函数沿X的第1个成对维计算。l Y=prctile(X，p,dim) 沿X的第dim维进行计算，Y的第dim维长度为length（p）。【例5-13】x=(1:5)*(1:5)x=1 2 3 4 52 4 6 8 103 6 9 12 154 8 12 16 20 5 10 15 10 25y= prctile(x，25 50 75)y=1.7500 3.5000 5.2500 7.0000 8.75003.0000 6.0000 9.0000 12.0000 15.00004.2500 .5000 12.7500 17.0000 21.25005.6 自助统计量用bootstrp函数计算数据重复取样的自助统计量。其调用格式和描述如下。l bootstat=bootstrap（nboot，bootfun，d1，d2）从输入数据集d1，d2等中提取nboot个自助数据样本并传递给bootfun函数进行分析。bootfun是一个函数句柄。nboot必须为正整数，并且每个输入数据集必须包含相同的行数n，每个自助样本包含n行，它们随机取自对应的输入数据集d1，d2等。输出bootstat的每一行包括将bootfun函数应用于一个自助样本时生成的结果。如果bootfun函数返回多个输出参数，只在bootstat中保存第1个。如果bootfun函数的第1个输出为矩阵，则该矩阵重塑为行矢量，以便保存到bootstat中。l bootstat，bootsam=bootstrap（）返回一个nxn的自助编号导入矩阵bootsam。bootsam中的每一列包含从原始数据集中提取出来组成对应自助样本的值的编号。例如，如果d1，d2等每个都包含了16个值，nboot=4，则bootsam是一个16x4的矩阵。第一列包含从d1，d2等数据集中提取出来形成前4个自助样本的16个值的编号，第2列包含4个自助样本的16个值的编号，以此类推。【例5-14】计算15个学生的LSAT分数和法学院GPA之间的关系。通过对这15个数据点进行重复采样，创建了1000个不同的数据集，然后计算每个数据集中这两个变量之间的相关关系。load lawdatabootstat,bootsam=bootstrp(1000,corrcoef,lsat,gpa);bootstat(1:5,:)ans =1 0.98743 0.98743 11 0.49182 0.49182 11 0.54587 0.54587 11 0.84577 0.84577 11 0.89587 0.89587 1bootsam(:,1:5)ans =13 3 11 8 1214 7 1 7 42 14 5 10 814 12 1 11 1110 15 2 12 142 10 13 5 155 1 11 11 99 13 5 10 315 15 15 3 315 11 1 2 43 12 7 8 1315 12 6 15 415 6 12 6 138 10 12 9 413 3 3 4 14hist(bootstat(:,2)生成直方图如图3-1所示。图3-1 直方图该直方图显示了整个自助样本的相关系数的变化。样本最小值为正，表示LSAT和GPA之间是相关的。5.7中心距K阶中心距可以用下式定义式中，为的期望。moment函数计算所有阶次的中心距。l m=moment(X,order) 返回由正整数order指定阶次的X中心距。对于矢量，moment(X,order) 函数返回X的元素的指定阶次的中心距。对于矩阵，moment(X,order) 返回每一列的指定阶次的中心距。对于多维数组，moment函数沿X的第1个成对维进行计算。l m=moment(X,order,dim) 沿X的第dim为进行计算。注意：一阶中心距为0，二阶中心距为用除数n(而非n-1)得到的方差，其中n为矢量X的长度或是矩阵X的行数。【例5-15】X=randn(6,5)X =0.53767 -0.43359 0.7254 1.409 0.488891.8339 0.34262 -0.063055 1.4172 1.0347-2.2588 3.5784 0.71474 0.6715 0.726890.86217 2.7694 -0.20497 -1.2075 -0.303440.31877 -1.3499 -0.12414 0.71724 0.29387-1.3077 3.0349 1.4897 1.6302 -0.78728m=moment(X,3)m =-1.1143 -0.99731 0.12339 -1.1023 -0.104555.8 相关函数用corrcoef函数计算样本数据的相关系数矩阵。该函数的调用格式如下。l R=corrcoef（X）返回输入矩阵X算的的相关系数矩阵R，输入矩阵的行为观测量，列为变量。相关系数矩阵R中的第(i，j)个元素与协方差矩阵C(=cov(X)有关，即有l R=corrcoef(x,y)与corrcoef(x y)一样，x和y为列矢量。l R,P=corrcoef()还返回矩阵P。P的值用于检验没有相关性的假设。如果P(i,j)较小，小于0.05，则相关性R(I,j)显著。l R,P,RLO,RUP=corrcoef()还返回大小与R相同的矩阵RLO和RUP，它们分别包含每个系数95%置信区间的下届和上届。l =corrcoef(,paraml,vall,param2,val2,)指定其他参数和它们的值。合法的值包括alpha和rows(或all/complete/pairwise)。【例5-16】生成4列随机数，使第4列与其他列有相关性。x=randn(30,4); %不相关的数据x(;,4)=sum(x,2) %引入相关性r,p=corrcoef(x) %计算样本相关性和p值i,j=find(p0.05); %查找相关性比较显著的数据i,j %显示它们(行或列)的编号r=1.0000 -0.3566 0.1929 0.3457-0.3566 1.0000 -0.1429 0.44610.1929 -0.1429 1.0000 0.51830.3457 0.4461 0.5183 1.0000p=1.0000 -0.0531 0.3072 0.0613-0.0531 1.0000 0.4511 0.01350.3072 0.4511 1.0000 0.00330.0613 0.0135 0.0033 1.0000ans=4 243 1 4 2 45.9 协方差矩阵用cov函数计算协方差矩阵。其语法格式为：l C=cov(X)l C=cov(X,Y)对于单一矢量，cov(X)返回一个包括方差的标量。对于行为观测量，列为变量的矩阵而言，cov(X)为协方差矩阵。计算方差的函数var(X)等价于diag(cov(X)。计算标准差的函数std(X)等价于sqrt(diag(cov(X)。cov(X,Y)等价于cov(X,Y)，其中，X，Y为长度相等的列矢量。cov函数的算法为：n,p=size(X);X=X-ones(n,1)*mean(X);Y=X*X/(n-1);5.10峰度和偏差5.10.1峰度样本的峰度由下式定义式中，为的期望。峰度用于度量样本数据偏离某分布的情况，正态分布的峰度为3。当样本数据的曲线峰值比正态分布的高时，峰度大于3；反之，比正态分布低时，峰度小于3。用kurtosis函数计算样本的峰度。l k=kurtosis(X) 返回X的样本峰度。对于矢量而言，kurtosis(X)函数为矢量X中元素的峰度。对于矩阵而言，kurtosis(X)函数为X的每一列返回一个样本峰度。对于多维数组，kurtosis函数沿X的第1个成对维进行计算。l k=kurtosis(X,flag) 指定是否校正系统偏差。将flag指定为0时，进行校正；指定1时，不校正。l kurtosis(X,flag,dim) 沿X的第dim维进行计算。注意：也有将峰度定义为计算值减3的，所以正态分布的峰度为0。【例5-17】X=randn(5,4)X =0.53767 -1.3077 -1.3499 -0.204971.8339 -0.43359 3.0349 -0.12414-2.2588 0.34262 0.7254 1.48970.86217 3.5784 -0.063055 1.4090.31877 2.7694 0.71474 1.4172k=kurtosis(X)k =2.7067 1.4069 2.3783 1.17