第四讲数据分析方法

第四讲 数据分析方法

第一节、数据拟合

问题：给定一批数据点（输入变量与输出变量的数据），需确定满足特定要求的曲线或 曲面。如果输入变量和输出变量都只有一个，则属于一元函数的拟合和插值；而若输入变量 有多个，则为多元函数的拟合和插值（有点回归分析的意思）

解决方案： （1） 若要求所求曲线（面）通过所给所有数据点，就是插值问题；

（2） 若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就

是数据拟合，又称曲线拟合或曲面拟合。 注意：插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同， 二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时 容易确定，有时则并不明显。

例 1：下面数据是某次实验所得，希望得到 X和 f之间的关系？ x f

1 1.5

2 3.9

4 6.6

7 11.7

9 15.6

12 18.8

13 19.6

15 20.6

17 21.1

曲线拟合问题最常用的解法——最小二乘法的基本思路

第一步：确定拟合的函数类型 y f (x;a1,a2,\",am )，其中a1,a2,\",am为待定系数。 （函数类型的确定可以根据内在的规律确定，如果无现成的规则，则可以通过散点图，联系 曲线的形状进行分析）

第二步：确定a1,a2,\",am的最小二乘准则：要求n个已知点(xi, yi )与曲线 y f (x) 的距离di的平方和 (yi f x ( )) i 2 最小 。

i1 n

用 MATLAB作拟合

1．多项式拟合。作多项式 y a0x

a=polyfit(x,y,m)—其中 x,y为给出的数据，m为多项式的次数。

a1x 多项式在 x处的值 y可用以下命令计算： m1 \"am拟合，可利用

y=polyval（a,x）

2．用 MATLAB作非线性最小二乘拟合

Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和 lsqnonlin。两个命令 都要先建立 M-文件 fun.m，在其中定义函数 f(x)。 （1）x = lsqcurvefit (‘fun’,x0,x‚ata,y‚ata);

（2）x =lsqcurvefit (‘fun’,x0,x‚ata,y‚ata,options);

（3）x = lsqcurvefit (‘fun’,x0,x‚ata,y‚ata,options,’‛ra‚’); （4）“x, options” = lsqcurvefit (‘fun’,x0,x‚ata,y‚ata,„);

（5）“x, options,funval” = lsqcurvefit (‘fun’,x0,x‚ata,y‚ata,„);

m

1

（6）“x, options,funval, ‟acob” = lsqcurvefit (‘fun’,x0,x‚ata,y‚ata,„); （7）x=lsqnonlin（‘fun’，x0）；

（8）x= lsqnonlin （‘fun’，x0，options）；

（9）x= lsqnonlin （‘fun’，x0，options，‘‛ra‚’）； （10）[x，options]= lsqnonlin （‘fun’，x0，„）；

（11）[x，options，funval]= lsqnonlin （‘fun’，x0，„）；

在使用 lsqcurvefit与 lsqnonlin命令时，共同的问题是要先知道函数的类型，而拟合其实是决 定函数中的待定系数。

第二节 插值

插值的基本问题：给出n个数对，(Pi, f (Pi)),i 1,2,\",n，求点 P处对应的函数值 f (P)。 一、一维插值

y *。常用的插值方法 x *处的函数值 已知 n 1个节点 (xi, yi ),i 0,1,\",n，求任意点

有拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite插值和三次样条插值。

分段线性插值：将各数据点用折线连接起来

多项式插值：求一个多项式通过所有数据点，可以假设出多项式的系数，最后通过 求解方程得到每个系数（拉格朗日插值，用n次多项式描述 n 1个点） 样条插值：分段多项式的光滑连接（三次样条插值） 牛顿插值：利用节点之间的各阶差商和差分构造多项式

Hermite插值：对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有 相同的一阶、二阶甚至更高阶的导数值

（1）MATLAB命令：y=interp1(x0,y0,x,'method')

method指定插值的方法，默认为线性插值。其值可为： 'nearest' 最近项插值 'linear' 'spline' 'cubic'

线性插值

立方样条插值 立方插值。

所有的插值方法要求 x0是单调的。

当 x0为等距时可以用快速插值法，使用快速插值法的格式为'*nearest'、'*linear'、'*spline'、 '*cubic'。

（2）三次样条插值在 Matlab中的实现

在 Matlab中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就 是采用非扭结（not-a-knot）条件。这个条件强迫第 1个和第 2个三次多项式的三阶导数相 等。对最后一个和倒数第 2个三次多项式也做同样地处理。

Matlab中三次样条插值也有现成的函数： y=interp1(x0,y0,x,'spline')； y=spline(x0,y0,x)； pp=csape(x0,y0,conds),

pp=csape(x0,y0,conds,valconds)，y=ppval(pp,x)。

其中 x0,y0是已知数据点，x是插值点，y是插值点的函数值。

对于三次样条插值，我们提倡使用函数 csape，csape的返回值是 pp形式，要求插值点

2

的函数值，必须调用函数 ppval。

pp=csape(x0,y0)：使用默认的边界条件，即 Lagrange边界条件。

pp=csape(x0,y0,conds,valconds)中的 conds指定插值的边界条件，其值可为：

'complete' 边界为一阶导数，一阶导数的值在 valconds参数中给出，若忽略 valconds 参数，则按缺省情况处理。

'not-a-knot' 非扭结条件

周期条件 'periodic'

'second' 边界为二阶导数，二阶导数的值在 valconds参数中给出，若忽略 valconds 参数，二阶导数的缺省值为[0, 0]。

'variational' 设置边界的二阶导数值为[0,0]。

对于一些特殊的边界条件，可以通过 conds的一个12矩阵来表示，conds元素的取值 为 0，1，2。

conds(i)=j的含义是给定端点i的 j阶导数，即 conds的第一个元素表示左边界的条件， 第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对 应的值由 valconds给出。

例 2 机床加工

待加工零件的外形根据工艺要求由一组数据(x, y)给出（在平面情况下），用程控铣床 加工时每一刀只能沿 x方向和 y方向走非常小的一步，这就需要从已知数据得到加工所要求 的步长很小的(x, y)坐标。

表中给出的 x, y数据位于机翼断面的下轮廓线上，假设需要得到 x坐标每改变 0.1时的

y坐标。试完成加工所需数据，画出曲线，并求出 x 0处的曲线斜率和13x 15范围

内 y的最小值。

x y

0 3 5 7 9 11 12 13 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2

14 15 1.0 1.6

要求用分段线性和三次样条计算。 解 编写以下程序：

x0=[0 3 5 7 9 11 12 13 14 15];

y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6]; x=0:0.1:15;

y2=interp1(x0,y0,x);

y3=interp1(x0,y0,x,'spline'); pp1=csape(x0,y0); y4=ppval(pp1,x);

pp2=csape(x0,y0,'second'); y5=ppval(pp2,x);

[x',y1',y2',y3',y4',y5'] subplot(2,2,1)

plot(x0,y0,'+',x,y2)

title('Piecewise linear') subplot(2,2,2) plot(x0,y0,'+',x,y3)

3

title('Spline1') subplot(2,2,3)

plot(x0,y0,'+',x,y4) title('Spline2') dx=diff(x); dy=diff(y3); dy_dx=dy./dx; dy_dx0=dy_dx(1) ytemp=y3(131:151); ymin=min(ytemp);

index=find(y3==ymin); xmin=x(index); [xmin,ymin]

二、二维插值：

前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若节点 是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的高程（节 点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插 值）。

（1）

网格节点：已知m

n个节点 (xi , y j , zij ),i 1,2,\",m; j 1,2,\",n，构造

一个二元函数 z f (x, y)通过全部已知节点，即 zij f (xi , y j )，再利用

z f (x, y)插值。

（2）

散乱节点：已知 n 个节点 (xi, yi, zi ),i 1,2,\",n ，构造一个二元函数

z f (x, y)通过全部已知节点，即 zi f (xi, yi )，再利用 z f (x, y)插值。

二维插值方法：

（1）最邻近插值：二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值 即为所求。

（2）分片线性插值：将四个插值点（矩形的四个顶点）处的函数值依次简记为： f (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f4 分两片的函数表达式如下： 第一片(下三角形区域)： f (x, y) 第二片(上三角形区域)： f (x, y)

f1 ( f 2 f1)(x xi ) ( f3 f 2)(y y j ) f1 ( f 4 f1)(y y j ) ( f3 f 4)(x xi )

（3）双线性插值：双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形 式如下： f (x, y) (ax b)(cy d)。其中有四个待定系数，利用该函数在矩形的四个顶点 （插值节点）的函数值，得到四个代数方程，正好确定四个系数。

用 MATLAB作网格节点数据的插值：

4

z=interp2(x0,y0,z0,x,y,’metho‚’) —method 可以选择 ‘nearest’：最邻近插值；‘linear’： 双线性插值；‘cubic’：双三次插值；缺省时, 双线性插值。

如果是三次样条插值，可以使用命令 pp=csape({x0,y0},z0,conds,valconds)，z=fnval(pp,{x,y}) 其中 x0,y0分别为m维和n维向量，z0为mn维矩阵， z为矩阵，它的行数为 x的维数， 列数为 y的维数，表示得到的插值，具体使用方法同一维插值。

例3 在一丘陵地带测量高程， x和 y方向每隔100米测一个点，得高程如下表，试插值 一曲面，确定合适的模型，并由此找出最高点和该点的高程。

x y

100 200 300 400

100 636 698 680 662

200 697 712 674 626 300 624 630 598 552 400 500 478 450 478 420 412 400 334 310

解 编写程序如下：

x=100:100:500; y=100:100:400;

z=[636 697 624 478 450; 698 712 630 478 420; 680 674 598 412 400; 662 626 552 334 310]; pp=csape({x,y},z')

xi=100:10:500;yi=100:10:400 cz1=fnval(pp,{xi,yi})

cz2=interp2(x,y,z,xi,yi','spline') [i,j]=find(cz1==max(max(cz1))) x=xi(i),y=yi(j),zmax=cz1(i,j)

用 MATLAB作散点数据的插值计算

cz =griddata（x,y,z,cx,cy,‘metho‚’）—method可以选择‘nearest’：最邻近插值；‘linear’ ： 双线性插值；‘cubic’：双三次插值；'v4'：Matlab提供的插值方法；缺省时, 双线性插值 例 4 在某海域测得一些点（x,y）处的水深 z由下表给出，在矩形区域（75,200）×(-50,150) 内画出海底曲面的图形。

x 129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5 y 7.5 141.5 23 147 22.5 137.5 85.5 –6.5 -81 3 56.5 –66.5 84 -33.5 z 4 8 6 8 6 8 9 9 8 8 9 9 8 4

解 编写程序如下：

x=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];

y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5]; z=[4,8,6,8,6,8,8,9,9,8,8,9,4,9]; xi=75:1:200;

5

yi=-50:1:150;

zi=griddata(x,y,z,xi,yi','cubic') subplot(1,2,1) plot(x,y,'*') subplot(1,2,2) mesh(xi,yi,zi)

第三节、回归分析方法

“回归”问题最早来源于生物界，英国生物学家兼统计学家高尔顿（Galton,1822-1911） 发现同一种族中儿子的平均高度介于其父亲的高度与种族平均高度之间。儿子的身高有返归 于种族平均身高的趋势，即回归于种族的平均身高。

回归分析是指对具有相关关系的现象，根据其关系形态，选择一个合适的数学模型，用 来近似地表示变量间的平均变化关系的一种统计方法。

回归分析的分类：按照回归模型中变量个数分（一元回归，多元回归）；按照回归曲线 的形态分（线性回归，非线性回归）；按照是否要求总体分布类型已知分（参数回归，非参 数回归）

一元线性回归特点

两个变量中，一个是自变量，一个是因变量

回归方程不是抽象的数学模型，而随机方程，可以进行实证 因果关系不明显时，应同时作两个回归方程 回归系数具有较强的经济含义

作为回归模型的因变量是随机变量，而自变量是确定性变量，即可控变量 1．一元线性回归模型 一般地，称由 y 0

1x 0

1

确定的模型为一元线性回归模型，记为

y E

x

2

0, D

其中固定的未知参数 0， 1称为回归系数，自变量 x也称为回归变量， Y 为 y 对 x的回归直线方程。

一元线性回归分析的主要任务

1．用试验值（样本值）对 0， 1和 作点估计； 2．对回归系数 0， 1作假设检验；

3．在 x x0处对 y作预测，对 y作区间估计。

回归系数的最小二乘估计

有n组观测值，(xi, yi ),i 1,2,\",n。

0 1

x称

yi

设 E

0 i

x1 i,i 1,2,...,n 0, Di 2 且12，...,

n相互

6

记Q Q(

01

,)

n i1

n

0和 1的估

i1

yi 0 1xi2 。最小二乘法就是选择

ˆ Q(ˆ0,ˆ1) min Q(0,1)。计算得到 计 ˆ0 ， 1使得

2

i ,0

1

ˆ0 y ˆ1x ˆ

1 1

其中 x xi, y n yi ， x 2

n i1 i1

n

n

xy xy x x

n

2 2

1

1 1 n

xi 2, xy xi yi 。（经验）回归方程为： n i1 n i1

yˆ ˆ0 ˆ1x y ˆ1(x x) 。

记 Q Q(ˆ0,ˆ )

e

1

n i1

yi

ˆ0

ˆ 1xi

n 2

i1

(yi y ˆ ) 2 ，称 Q e为残差平方和或剩

i

ˆ e Qe (n 2) 。

余平方和。 2的无偏估计为

回归方程的显著性检验 对回归方程 Y 检验。

假设 H 0 :

1

0

1

2

x的显著性检验，归结为对假设 H 0 : 1 0;H1 : 1 0进行

0被拒绝，则回归显著，认为 y与 x存在线性关系，所求的线性回归方

程有意义；否则回归不显著，y与 x的关系不能用一元线性回归模型来描述，所得的回归方 程也无意义。

U

~ F(1,n 2)，其中 U F检验法：当 H 0成立时， F

Qe /(n 2)

归平方和）。若 F F1(1,n 2)，拒绝 H 0，否则就接受 H 0。

n i1

yˆi y2（回

回归系数的置信区间

0和 1置信水平为1 的置信区间分别为

ˆ0 t (n 2)ˆ e

12

2

1 x ,ˆ0 t (n 2)ˆ e

1n Lxx

2

1 x 2

和

n Lxx

ˆ t (n 2)ˆ e / Lxx ,ˆ1 t (n 2)ˆ e / L ；

12 12 1 xx

7

2

的置信水平为1 的置信区间为

Qe Qe 2 (n 2) , (n 2 2)

1 2

2

。

预测

用 y 的回归值 yˆ0 ˆ0 ˆ1x

0

0作为 0的预测值。 0的置信水平为 的预测区间为

2

y y 1

yˆ0 (x0), yˆ0 (x0)，其中(x0) ˆ et (n 2) 11

12

x0 xn

Lxx

。特别，当 n很

大 且 x0 在 x 附 近 取 值 时 , y 的 置 信 水 平 为 1的 预 测 区 间 近 似 为

yˆ ˆ eu 12 , yˆ ˆ eu 12 。

2．可线性化的一元非线性回归（需要配曲线）

先对两个变量 x和 y作n次试验观察得(xi , yi ),i 1,2,...,n画出散点图，根据散点图确 定须配曲线的类型.然后由n对试验数据确定每一类曲线的未知参数a和b。采用的方法是通 过变量代换把非线性回归化成线性回归，即采用非线性回归线性化的方法。

通常选择的六类曲线如下 （1）双曲线 1 a

y

b x

（2）幂函数曲线 y axb，其中 x 0,a 0 （3）指数曲线 y aebx其中参数 a 0 （4）倒指数曲线 y aeb / x其中 a 0 （5）对数曲线 y a bln x, x 0

1

（6）S型曲线 y a bex

3．多元线性回归

一般称

Y X

E() 0,COV(,)

2

2

I n

为高斯—马尔柯夫线性模型(k 元线性回归模

型)，并简记为(Y, X,I n )。

8

其 中 Y

y1 ... ... y n

， X

1 x11 x12 ... x1k 1 x21 x22 ... x2k

... ... ... ... ... 1 xn1 xn2 ... xnk

0

1

，

1

...

k

，

2

...

n

，

y

0 11

x...

kk称为回归平面方程。

2

x

线性模型(Y, X,

I n )考虑的主要问题是：

（1）用试验值（样本值）对未知参数 和2 作点估计和假设检验，从而建立 y与 x1, x2,..., xk 之间的数量关系；

（2）在 x1 x01, x2 x02,..., xk x0k处对 y的值作预测与控制，即对 y作区间估计。

对 i和

2

作估计，用最小二乘法求 0,...,

n i1

k 的估计量：作离差平方和

Q

选择 0,...,

k 使

yi

0 1i1

x...

k ik 2

x

Q达到最小。解得估计值 ˆ X T X

1

X TY，得到的 ˆ 代入回归平面 i

方程得 y ˆ0 ˆ1x1 ...ˆk xk ，称为经验回归平面方程. ˆi称为经验回归系数

4．多元非线性回归

多项式回归

设变量 x,Y 的回归模型为

Y

0 1

x

i

2

x ...

p

x

2

其中 p 是已知的， (i 1,2,\", p) 是未知参数， 服从正态分布 N(0,x

k

2 p

) 。

Y

0 1

x

2

2 x ...

k称为回归多项式，上面的回归模型称为多项式回归。

令 xi xi，i=1，2，„，k多项式回归模型变为多元线性回归模型.

多元线性回归中的检验与预测

假设 H 0 :

0

1

...

k

0。 U / k

~ F(k,n k 1)。如果 F F1(k,n k 1)，

当 H 0成立时， F

Qe /(n k 1)

则拒绝 H 0，认为 y 与 x1,\", xk之间显著地有线性关系；否则就接受 H 0，认为 y 与 x1,\", xk 之间线性关系不显著。

9

求 出 回 归 方 程 yˆ ˆ0 ˆ1x1 ...ˆk xk ， 对 于 给 定 自 变 量 的 值

* yˆ ˆ0 ˆ1x1

*

*

,..., xk x 1

， 用

...

ˆk xk 来预测 y

*

0 11

x

*

...

k k

x

*

y

.称 *为 *的点预测。

yˆ

5．逐步回归分析

“最优”的回归方程就是包含所有对Y 有影响的变量, 而不包含对Y 影响不显著的变量 回归方程。选择“最优”的回归方程有以下几种方法：

（1）从所有可能的因子（变量）组合的回归方程中选择最优者； （2）从包含全部变量的回归方程中逐次剔除不显著因子； （3）从一个变量开始，把变量逐个引入方程； （4）“有进有出”的逐步回归分析

逐步回归分析法的思想：

（1）从一个自变量开始，视自变量Y 作用的显著程度，从大到小地依次逐个引入回归方程。 当引入的自变量由于后面变量的引入而变得不显著时，要将其剔除掉。

（2）引入一个自变量或从回归方程中剔除一个自变量，为逐步回归的一步。

（3）对于每一步都要进行Y 值检验，以确保每次引入新的显著性变量前回归方程中只包含 对Y 作用显著的变量。

（4）这个过程反复进行，直至既无不显著的变量从回归方程中剔除，又无显著变量可引入 回归方程时为止。 6．回归 MATLAB命令： （1）多元线性回归

确定回归系数的点估计值：b=regress( Y, X )

ˆ 0 ˆ

1

b

Y

... ˆ

Y 1 Y 2 X ... Y n

1 x11 x12 ... x1p

1 x21 x22 ... x2 p ... ... ... ... ... 1 xn1 xn2 ... xnp

p

求回归系数的点估计和区间估计、并检验回归模型：

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

Stats给出用于检验回归模型的统计量，有三个数值：相关系数r2、F值、与F对应的概率p。 判别规则：相关系数r2越接近 1，说明回归方程越显著； F F1(k,n k 1)时拒绝H ，

0

F越大，说明回归方程越显著；与F对应的概率p时拒绝H0，回归模型成立。

（2）多元二项式回归

rstool（x，y，’mo‚el’, alpha）

‘mo‚el’：由下列 4个模型中选择 1个（用字符串输入，缺省时为线性模型） linear（线性）： y

0

11

x\"

mm；

x

10

purequadratic（纯二次）： y

n

0

11

x\"x

mm

x

jj2j j1

x

interaction（交叉）： y

0 11

x\"

11

mm

1jkm

jk

x jxk

jk

quadratic（完全二次）： y （3）非线性回归

确定回归系数的命令：

0

x\"

mm

1j,km

x

x jxk

[beta，r，J]=nlinfit（x，y，’mo‚el’, beta0）：r-返回残差。 nlintool（x，y，’mo‚el’, beta0，alpha）

[Y，DELTA]=nlpredci（’mo‚el’, x，beta，r，J）

求 nlinfit 或 nlintool所得的回归函数在 x处的预测值 Y及预测值的显著性为 1-alpha的置信 区间。

11