第四讲 统计分析及数据处理

第四讲 统计分析及数据处理

一、描述性统计分析

在统计工具箱中，一般将向量视为随机样本，如果分析对象是矩阵，则将矩阵的每一列视为随机变量的样本。

mean(x) 当x 为向量，返回向量的均值，当x 为矩阵，返回矩阵每列元素均值构成的行向量。 特例

min 、max 、sort 、median 、std 、var 、sum 、prod 、cumsum 、cumprod 等函数的用法与mean 类似。

cov(x,y) 求样本x 与y 的协方差。（1

1()()1n

x y i i i S x x y y

n ==---∑）

cov(X) 这里X 是矩阵，将矩阵各列看成一个样本，求样本的协方差矩阵。

corrcoef 用法与cov 类似，求相关系数或相关系数矩阵。 （xy xy x y

S R s s =

）

[y,I]=sort(x) 当x 为向量，y 返回x 的升序排列，I 返回y 各元素原来的位置，即y=x(I)，当x 为矩阵，分别对各列排序。如果要按降序排列，[y,I]=sort(x,’descend ’)，如果要将矩阵的各行排序，[y,I]=sort(x,2)。

[y,I]=sortrows(x,k) 将矩阵x 的各行按第k 列的升序排列.

y=prctile(x,p) 当x 为向量，y 返回x 的p%下分位数，当x 为矩阵，分别求各列的下分位数。

相关知识：样本百分位数（percentiles ） 将样本()12,,,n x x x 按升

序排列为()(1)(2)(),,,n x x x ，则样本的p%下分位数表示这样一个值p x ：比p x 小的样本占总样本数的p%，当p=50时的下分位数称为中位数。 trimmean(x,p) 剔除上下各(p/2)%数据以后的均值。

例4.1：data=[11 57 29;13 54 78;10 66 54;9 46 37;16 75 24;15 70 56;9 40 31];

观察下列指令的运行结果： mean(data) std(data)

var(data) %考虑var(data)与std(data)的关系 cov(data) corrcoef(data) [y,I]=sort(data) [y,I]=sortrows(data ，2) prctile(data,[25,50]) median(data) trimmean(data,20) 二、概率函数

y=pdf(dist,x,p1,p2,…) 返回以p1,p2,…为参数的服从dist 分布的概率密度函数（Probability density functions ）在x 处的值。

常用分布及相关参数为：

dist 说明 参数 ‘bino ’ or ‘Binomial ’ 二项分布B(n,p) n,p

‘chi2’or ‘Chisquare’χ2分布χ2(n) n

‘exp’or ‘Exponential’指数分布μ

‘f’or ‘F’F分布F(n1,n2) n1,n2

‘norm’or ‘Normal’正态分布N(μ,σ) μ,σ

‘poiss’or ‘Poisson’泊松分布P(λ)，λ

‘t’or ‘T’t-分布T(n) n

‘unif’or ‘Uniform’连续型均匀分布U(a,b) a,b

‘unid’or ‘Discrete Uniform’离散型均匀分布U(n)

y=cdf(dist,x,p1,p2,…) 返回dist分布在x处的分布函数值

(Cumulative distribution functions)。

y=icdf(dist,x,p1,p2,…) 求逆分布函数值（即下分位数）

y=normpdf(x,mu,sigma) 正态分布N(mu,sigma)的概率密度值

p=normcdf(x,mu,sigma) 正态分布N(mu,sigma)的分布函数值x=norminv(p,mu,sigma) 正态分布N(mu,sigma)的p下分位数。例4.2：观察下列指令，并解释其概率意义。

>> a=normpdf(90,80,10)-normpdf(70,80,10)

>> b=normcdf(90,80,10)-normcdf(70,80,10)

>> p=(1-b)/2;c=norminv(p,80,10)

>> d=norminv(1-p,80,10)

>> pdf('norm',90,80,10)

>> p=cdf('norm',90,80,10)

>> icdf('norm',p,80,10) >> icdf('t',0.95,20) >> pdf('bino',5,10,0.5) 三、统计推断 1、参数估计

[mu,sigma,muci,sigmaci]=normfit(x,alpha) 正态总体N(μ, σ2)的参数估计，x 为样本数据，alpha 为显著性水平（缺省值为0.05），输出mu 和sigma 为总体均值和标准差的点估计值，muci 和sigmaci 为相应的区间估计。

phat=mle(data,’distribution ’,dist) 根据数据data ，求分布dist 中参数的极大似然估计phat 。 2、假设检验

ttest 单样本的t-检验，主要格式：

[h,p]=ttest(x,m,name,value) 检验零假设0:H m μ=，其中x 是来自方差未知的正态总体。输出h 的值为0或1，若h=1，则拒绝零假设，否则若h=0，则不能拒绝零假设；同样，p 值的大小也可以对假设进行推断，若p，tail=-1表示1:H m μ<。 ttest2 两个独立的正态总体的均值是否相等的假设检验。

[h,p]=ttest2(x,y,name,value) 其中x 是来自正态总体N(μ1, σ

12

)

的样本，y 是来自N(μ2, σ

2

2

)的样本，方差均未知，零假设012:H μμ=。

成对参数选项与ttest 相同之外，还有属性’vartype ’，可选值为’equal ’(表示方差相等)和’unequal ’（表示方差不等）。 3、分布的检验

h=jbtest(x) 检验样本数据x 是否来自均值和方差未知的正态总体，返回值h=1表示拒绝来自正态总体的假设，h=0则做出相反的结论。

h=kstest(x) 检验样本数据x 是否来自标准正态分布的总体。 h=kstest2(x,y) 检验两个数据样本x 和y 是否具有相同的连续分布，若h=1,则拒绝原假设，即x 和y 是来自不同的分布。

例4.3：>> x=randn(200,1); %生成来自标准正态分布的200个随机数

>> [h,p]=kstest(x) %检验是否服从标准正态分布 h = 0 %接受原假设，即服从标准正态分布 p = 0.7723 % p>0.05 ,同样表示接受原假设 >> [mu,sigma,muci,sigmaci]=normfit(x) %参数估计 mu = 0.0252 %均值的点估计值 sigma = 1.0883 %标准差的点估计值 muci = -0.1265

0.1769 %均值的区间估计

sigmaci = 0.9911

1.2068 %标准差的区间估计 >> h=ttest(x) %检验均值是否为0

h = 0 %接受原假设

例4.4：>> x=binornd(10,0.8,200,1); %生成服从二项分布B(10,0.8)

的200个随机数

>> mle(x,'distribution','bino','ntrials',10) %对B(10,p)中参数

p的极大似然估计

ans = 0.8065

例4.5：测得一批钢件的20个样品的屈服点（单位：T/mm2）为

4.98

5.11 5.20 5.11 5.00 5.61 4.88 5.27 5.38 5.20

5.46 5.27 5.23 4.96 5.35 5.15 5.35 4.77 5.33 5.54

并假设屈服点服从正态分布。已知总体均值为5.20，试判断该样本

均值与总体均值之间是否有显著差异。即检验原假设

0: 5.20

H x

>> x=[4.98 5.11 5.20 5.11 5.00 5.61 4.88 5.27 5.38 5.20 ...

5.46 5.27 5.23 4.96 5.35 5.15 5.35 4.77 5.33 5.54]; >> [h,p,ci]=ttest(x,5.20) %总体均值的t-检验

h = 0 %接受原假设，即样本均值等于5.20

p = 0.8796 % p>0.05 ，同样表示接受原假设

ci = 5.1052 5.3098 %样本均值的置信区间

例4.6：对两种不同的水稻品种A，B分别统计了8个地区的单位面积产量（单位：kg）。

品种A：86 87 56 93 84 93 75 79

品种B：80 79 58 91 77 82 76 66

试检验两个水稻品种的单位面积产量之间是否有显著差异。

即检验原假设

0:

A B

Hμμ

=

>> x=[86 87 56 93 84 93 75 79];

>> y=[80 79 58 91 77 82 76 66 ];

>> [h,p,ci]=ttest2(x,y) %两个正态总体均值相等的t-检验

h = 0 %接受原假设，即两个品种没有显著差异

p = 0.3393 %p>0.05，接受原假设

ci = -6.4236 17.4236 %均值差的置信区间

例4.7：下面列出了84个伊特拉斯坎（Etruscan）人男子的头颅的最大宽度（mm），试检验这些数据是否服从正态分布。

141 148 132 138 154 142 150 146 155 158

150 140 147 148 144 150 149 145 149 158

143 141 144 144 126 140 144 142 141 140

145 135 147 146 141 136 140 146 142 137

148 154 137 139 143 140 131 143 141 149

148 135 148 152 143 144 141 143 147 146

150 132 142 142 143 153 149 146 149 138

142 149 142 137 134 144 146 147 140 142

140 137 152 145

在命令窗口输入：

>> x=[141 148 132 138 154 142 150 146 155 158 ...

150 140 147 148 144 150 149 145 149 158 ...

143 141 144 144 126 140 144 142 141 140 ...

145 135 147 146 141 136 140 146 142 137 ... 148 154 137 139 143 140 131 143 141 149 ... 148 135 148 152 143 144 141 143 147 146 ... 150 132 142 142 143 153 149 146 149 138 ... 142 149 142 137 134 144 146 147 140 142 ... 140 137 152 145];

>> [h,p,jbtest,cv]=jbtest(x) %正态总体分布的假设检验 h = 0 %接受原假设，即服从正态分布 p = 0.5000 %p>0.05，与上述结论一致 jbtest = 0.7399 %检验统计量

cv = 5.3417 %拒绝零假设的临界值，即拒绝域为

（5.3417，+∞）

四、插值和曲线拟合

在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样本点，要求得到变量之间的函数关系或得到样本点之外的数据。解决此类问题的方法一般有插值、拟合和回归分析等。

设有一组实验数据0011(,),(,),(,)n n x y x y x y ，当原始数据精度较高，要求确定一个简单函数()y x ?=（一般为多项式或分段多项式）通过各数据点，即(),0,,i i y x i n ?== ，称为插值问题。

另一类是拟合问题，当我们已经有了函数关系的类型，而其中参数未知或原始数据有误差时，我们确定的初等函数()y x ?=并不要求

经过数据点，而是要求在某种距离度量下总体误差达到最小，即

(),0,,i i i y x i n ?ε=+= ，且20

n

i i ε=∑达到最小值。

对同一组实验数据，可以作出各种类型的拟合曲线，但拟合效果有好有坏，需要进行有效性的统计检验，这类问题称为回归分析。 1、插值(interpolation)

常用的插值方法有分段线性插值、分段立方插值、样条插值等。 ⑴一元插值

yi=interp1(x,y,xi,method)

对给定数据点(x,y)，按method 指定的方法求出插值函数在点（或数组）xi 处的函数值yi 。其中method 是字符串表达式，可以是以下形式：

'nearest' ——最邻近点插值

'linear' ——分段线性插值（也是缺省形式） 'spline' ——分段三次样条插值 'cubic' 分段立方插值

例4.8：在一天24小时内，从零点开始每间隔2小时测得环境温度数据分别为（℃）：

12，9，9，10，18，24，28，27，25，20，18，15，13 用不同的插值方法估计中午1点（即13点）的温度，并绘出温度变化曲线。 >> x=0:2:24;

>> y=[12 9 9 10 18 24 28 27 25 20 18 15 13];

>> y_linear=interp1(x,y,13),y_nearest=interp1(x,y,13,'nearest')

>> y_cubic=interp1(x,y,13,'cubic'),y_spline=interp1(x,y,13,'spline')

>> xx=0:0.1:24;

>> y1=interp1(x,y,xx); y2=interp1(x,y,xx,'nearest');

>> y3=interp1(x,y,xx,'cubic');y4=interp1(x,y,xx,'spline');

>> subplot(2,2,1),plot(x,y,'or',xx,y1)

>> subplot(2,2,2),plot(x,y,'or',xx,y2)

>> subplot(2,2,3),plot(x,y,'or',xx,y3)

>> subplot(2,2,4),plot(x,y,'or',xx,y4)

⑵二元插值

zi=interp2(X,Y,Z,xi,yi,method)/

已知数据点（X，Y，Z），求插值函数在（xi,yi）处的函数值zi,插值方法method同interp1。这里要求X,Y,Z是同维矩阵，且X，Y是网格矩阵，或者X是与Z列数相同的行向量，Y是与Z行数相同的列向量。

例4.9：测得平板表面5 3网格点处的温度分别为

试作出平板表面的温度分布图

>> x=1:5;y=1:3;z=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];

>> xx=1:0.1:5;yy=1:0.1:3;yy=yy';

>> zz=interp2(x,y,z,xx,yy,'cubic');

>> mesh(xx,yy,zz)

3、不规则点的插值

若数据是不规则的，即数据不能构成矩阵形式，从而不能用interp2函数进行插值。

zi=griddata(x,y,z,xi,yi,method)

这里,x,y,z为同维向量，表示已知数据点的坐标，xi,yi是行向量和列向量，返回值zi为在meshgrid(xi,yi)网格矩阵处的函数值。method 可选择’linear’,’nearest’,’cubic’。

例4.10：假如上例中的数据残缺不全

>> x=[3 4 5 1 3 4 1 2 5];y=[1 1 1 2 2 2 3 3 3];

>>z=[80 82 84 79 61 65 84 84 86];

>> xx=1:0.1:5;yy=1:0.1:3;

>> zz=griddata(x,y,z,xx,yy','cubic');

>> mesh(xx,yy',zz)

2、拟合(Fit) ⑴多项式拟合

p=polyfit(x,y,n) 用n 次多项式拟合向量数据(x,y)。 例4.11：拟合下列数据

>> x=[0.1 0.2 0.15 0 -0.2 0.3];y=[0.95 0.84 0.86 1.06 1.50 0.72]; >> p=polyfit(x,y,2);

>> xx=-0.2:0.01:0.3;yy=polyval(p,xx); >> plot(x,y,'or',xx,yy) ⑵曲线拟合

当经验函数不是多项式，而是其它类型的函数时，可以用lsqcurvefit 函数对拟合函数中的未知参数进行估计。

[c,Q]=lsqcurvefit(fun,c0,xdata,ydata)

fun 是经验拟合函数，含有未知参数，即具有形式fun(c,x)，c0是未知参数的预估计值，(xdata,ydata)是已知实验数据。输出c 是参数c 的估计值，Q 是残差平方和。 例4.12：已知数据表

用适当的曲线进行数据拟合。

先画散点图，根据散点图确定拟合曲线为对数函数ln b t y a += >> t=1:16;

>> y=[4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 ... 10.55 10.58 10.6]; >> plot(t,y,'or')

>> f=@(c,t)c(1)+c(2)*log(t) %建立拟合函数 >> c=lsqcurvefit(f,[1,1],t,y) %求未知参数 >> tt=1:0.1:16;yy=f(c,tt); >> hold on >> plot(tt,yy)

思考：上述问题是否可以转化为多项式拟合问题？

⑶拟合工具箱

Matlab 中的拟合工具箱是一个更方便、更直观进行曲线拟合的图形界面，用cftool 指令打开拟合工具箱。

拟合效果主要看2个参数：SSE （误差平方和）和R-Square ，SSE 越接近0，R-Square 越接近1，拟合效果越好。 五、回归分析 1、多元线性回归

问题：设有因变量y 和p 个自变量12,,p x x x L ，它们具有某种线性关系 1122p p y x x x βββε=+++L 其中12,,p βββ 为待定系数，ε为随机误差。现有容量为n 观测数

据,,1,2,,1,2,i ij y x i n j p ==L

L ，怎样确定待定系数12,,p βββ ，

并进行有效性检验？

将样本代入关系式，得Y X βε=+，其中

1111211221

22221

2,,p p p n n np p y x x x y x x x Y X y x x x ββββ??????

?????????

???===????????????????????????

L L M L

L M L

Matlab 求解：

[b,bint,r,rint,stats]=regress(Y ,X,alpha)

其中alpha 为显著性水平，省略时为默认值0.05；输出变量中，b 为β的参数估计值，bint 为b 的置信度为1-alpha 的置信区间，r 为残差向量Y X β-，也即ε，rint 为ε的置信区间，stats 是包含3个元素的检验统计量，分别是R-square ：相关系数R 的平方，F-统计量和p 值。 回归效果：R-Square 越接近1，p 值越接近0（一般要求p<0.05）。 例4.13：某种水泥在凝固时放出的热量（单位:卡/克）Y 与水泥中下列4种化学成分所占的百分比有关：

x1:233CaO Al O ?; x2: 23CaO SiO ?; x3:23234.CaO Al O Fe O ?; x4:22CaO SiO ? 现测得13组数据，见表，要求建立热量与水泥化学成分之间的经验回归关系式。

输入数据，可以先建立2个全零矩阵x=zeros(13,4);y=zeros(13,1);然后将表中的数据直接复制、粘贴到相应位置。 >> [b,bint,r,rint,stats]=regress(y,x) 最后得到的回归方程为：

12342.193 1.15330.75850.4863y x x x x =+++

如果回归方程是形式：01122p p y x x x ββββε=++++ ，相当于增加一个变量0x ，001122p p y x x x x ββββε=++++ ，而01x ≡。如上例

>> x=[ones(13,1),x];

>> [b,bint,r,rint,stats]=regress(y,x) 得回归关系式：

123462.4054 1.55110.51020.10190.1441y x x x x =+++-

2、非线性回归

[beta,R,J,CovB,MSE] = nlinfit(X,Y ,modelfun,beta0)

其中X 为因素数据矩阵，每一列一个变量，Y 为响应数据向量，

modelfun 表示模型的函数句柄或匿名函数，此函数的函数声明形式为：y=modelfun(beta,x)，其中beta 为待估参数；

beata0为参数迭代初值。输出部分：beta 为参数估计值，R 为残差，J 返回用于估计预测误差的Jacobi 矩阵，CovB 是参数beta 的协方差矩阵，MSE 是误差项ε的方差σ2估计值。

例4.14：（化学反应速度与反应物含量）在研究化学反应过程中，建立了一个反应速率和反应物含量的数学模型，形式为：

3

125

213243

1x x y x x x βββββ-

=

+++

其中12345,,,,βββββ是未知参数，123,,x x x 是三种反应物的含量，Y 是反应速率。今测得一组数据如下表，试确定参数12345,,,,βββββ。已给其参考值为（1,0.05,0.02,0.1,2）。

先编写经验函数 modelfun.m

function y=modelfun(beta,x)

b1=beta(1);b2=beta(2);b3=beta(3);b4=beta(4); b5=beta(5);

x1=x(:,1);x2=x(:,2);x3=x(:,3);

y=(b1*x2-x3/b5)./(1+b2*x1+b3*x2+b4*x3); 将表中123,,x x x 的三列数据导入到变量x ，反应速率Y 的数据导入到变量y 。

输入>> beta0=[1,0.05,0.02,0.1,2];

方法一：用非线性拟合函数lsqcurvefit 进行参数估计： >> [beta,Q]=lsqcurvefit(@modelfun,beta0,x,y) 方法二：用非线性回归函数nlinfit 。

>> [beta,R,J,CovB,MSE]=nlinfit(x,y,@modelfun,beta0) 注意这两个函数的调用格式的差异，而且Q=norm(R,2)^2。

上机练习

1、以下是100次刀具故障记录，即故障出现时该刀具完成的零件数。分析这批数据是否服从正态分布，并求其均值和均方差。注意，由于记录失误，其中可能有些数据是错误的，要对此进行适当处理。 459 362 624 542 509 584 433 748 815 505 612 452

434 982 640 782 742 565 706 593 680 926 653 164 487 734 608 428 1153 593 844 527 552 513 781 474 388 824 538 862 659 775 859 755 649 697 515 628 954 771 609 2 960 885 610 292 837 473 677 358 638 699 634 555 570 84 416 606 1062 484 120 447 654 564 339 280 246 687 539 790 581 621 724 531 512 577 496 468 499 544 645 764 558 378 765 666 763 217 715 310 851

2、下表给出了1930年各国人均年消耗的烟支数以及1950年男子死于肺癌的死亡率。（注：研究男子的肺癌死亡率是因为在1930年左右几乎极少的妇女吸烟，记录1950年的肺癌死亡率是因为考虑到吸烟的效应要有一段时间才能显现）

⑴画出该数据散点图；

⑵该散点图是否表明在吸烟多的人中间肺癌死亡率较高？

⑶计算两列数据的相关系数。

3、某保健食品商声称学生服用该保健食品一个月后能提高他们的数学能力和成绩，为了查明此保健食品是否真的有作用，设计了一次实验，随机地选取500名学生，并将他们随机地均分为两个组，甲组服用保健食品，乙组服用模样和品味与保健食品一样的葡萄糖丸，两组同学都以为自己在服用保健食品。一个月后进行一次数学考试，结果甲组的平均分是73分，标准差是18分，乙组的平均分是71分，标准差是17分，问其间的差异是由于机会变异引起的还是保健食品真的起了作用？

4、在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。

5、已知飞机下轮廓线上数据如下，求x每改变0.1时的y值。

X 0 3 5 7 9 11 12 13 14 15

Y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 6、在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要

避免进入。

7、山区地貌：在某山区测得一些地点的高程如下表：(平面区域1200<=x<=4000,1200<=y<=3600)，试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。

8、某商品的需求量与消费者的平均收入、上平的价格的统计数据如下，建立回归模型并进行检验，预测平均收入为1000，价格为6时的商品需求。