1 / 116

数据处理专题

数据处理专题. 数据处理是指用简明而严格的方法把获得的实验数据所代表的事物内在的规律提炼出来,得出结果的加工过程,包括数据记录、描绘曲线,从带有误差的数据中提取参数,验证和寻找经验规律,外推实验数据等等。本章介绍一些最基本的数据处理方法。. 数据处理的过程:. 1 、获得数据(标准化处理)。. 2 、将数据分类(聚类分析)。. 3 、提取主要影响因素( 主成分分析)。. 4 、数据分析(相关性分析,回归分析)。. 聚类分析.

demont
Télécharger la présentation

数据处理专题

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. 数据处理专题 数据处理是指用简明而严格的方法把获得的实验数据所代表的事物内在的规律提炼出来,得出结果的加工过程,包括数据记录、描绘曲线,从带有误差的数据中提取参数,验证和寻找经验规律,外推实验数据等等。本章介绍一些最基本的数据处理方法。

  2. 数据处理的过程: 1、获得数据(标准化处理)。 2、将数据分类(聚类分析)。 3、提取主要影响因素( 主成分分析)。 4、数据分析(相关性分析,回归分析)。

  3. 聚类分析 聚类也就是分类,在社会经济领域中存在大量的分类问题,比如三十个省市自治区独立核算工业企业经济效益进行分析,一般不是逐个省市自治区去分析,而较好的做法是选取具有代表性的指标如,百元固定资产实现利税,资金利税率、产值利税率、百元销售收入实现利润、全员劳动生产率等等,根据这些指标对省市自治区分类,然后根据分类结果对企业经济效益进行综合评价。 聚类分析方法 系统聚类法、有序样品聚类法、动态聚类法、模糊聚类法、图论聚类法、聚类预报法等。我们着重讲述系统聚类法。对样品分类成Q分类,对指标分类称R分类

  4. 聚类的三种尺度: 1、间隔尺度:变量是用连续量来表示,如长度、重量等 2、有序尺度:用一些等级来表示。如上中下三等。 3、名义尺度:既没有数量表示也没有次序表示。如红黄蓝三色等 我们通过距离来分类。方法有:最短距离法、最长距离法、中间距离法、重心法等。我们用最短距离法来讲述,其它方法读者自己翻阅相关的多元统计教材。

  5. 最短距离法步骤如下: 【1】定义样品之间的距离,计算样品两两距离,得一距离记为D(0) 开始每个样品自成一类,显然这时Dij =dij。其中D表示类G之间的距离,d表示样品之间的距离。 【2】找出D(0)的非对角线最小元素,设为Dpq,则将Gp和Gq合并为一新类,记为Gr 。 【3】给出计算新类与其他的类的距离公式: 距离公式有:欧氏距离,马氏距离,兰氏距离等。我们一般用马氏距离,应为它即排除了各指标之间相关性的干扰,而且还不受各指标量纲的影响。 两个样本间的距离定义: 其中,Xi 为样品的p个指标组成的向量。 协方差阵的逆矩阵

  6. 协方差阵定义如下: 样品到总体的距离定义: 总体均值向量

  7. Dkr=min{Dkp,Dkq}将D(0)中的第p、q行及p、q列用上面公式并成一个新行新列,新行新列对应Gr,所得到得矩阵记为D(1) 【4】对D(1)重复上述对D(0)的(2)(3)两步得D(2);如此下去,直到所有的元素并为一类。 注意:如果某一步中非对角线最小的元素不止一个,则对应这些最小元素的类可以同时合并。 为了大家便于掌握我们举例如下:

  8. 例:设抽取五个样品,每个样品只测一个指标,它们是1,2,3.5,7,9,试用最短距离法对这五个样品进行分类。例:设抽取五个样品,每个样品只测一个指标,它们是1,2,3.5,7,9,试用最短距离法对这五个样品进行分类。 解:我们距离选用我们所熟悉的绝对值距离。

  9. 最终我们分为两类比较合适,{x1,x2,x3}与{x4,x5}

  10. matlab做聚类分析 分步聚类:(1)找到数据集合中变量两两之间的相似性和非相似性,用pdist函数计算变量之间的距离;(2)用 linkage函数定义变量之间的连接;(3)用 cophenetic函数评价聚类信息;(4)用cluster函数创建聚类。 Step1 寻找变量之间的相似性 用pdist函数计算相似矩阵,有多种方法可以计算距离,进行计算之前最好先将数据用zscore函数进行标准化。 X=[1,2,3.5,7,9] X2=zscore(X); %标准化数据 Y2=pdist(X2); %计算距离 Step2定义变量之间的连接 Z2=linkage(Y2); Step3 评价聚类信息 C2=cophenet(Z2,Y2); //0.94698 Step4 创建聚类,并作出谱系图 T=cluster(Z2,2); H=dendrogram(Z2);%画出聚类图

  11. 为了更深入了解我国人口的文化程度状况,1990年全国人口普查数据对全国30个省直辖市、自治区进行聚类分析。分析选用了三个指标:【1】大学以上文化程度的人口占全部人口的比例(DXBZ);【2】初中以上文化程度的人口占全部人口的比例(CZBZ);【3】文盲半文盲的人口占全部人口的比例(WMBZ);分别用来反映较高、中等、较低文化程度人口的状况,原始数据如附件: clear clc X=load('data1.txt') Y2=pdist(X);%计算距离 Z2=linkage(Y2); C2=cophenet(Z2,Y2); T=cluster(Z2,4); H=dendrogram(Z2);%画出聚类图

  12. pdist函数 调用格式:Y=pdist(X,’metric’) 说明:用 ‘metric’指定的方法计算 X 数据矩阵中对象之间的距离。 X:一个m×n的矩阵,它是由m个对象组成的数据集,每个对象的大小为n。 metric’取值如下: ‘euclidean’:欧氏距离(默认);‘seuclidean’:标准化欧氏距离; ‘mahalanobis’:马氏距离;‘cityblock’:布洛克距离; ‘minkowski’:明可夫斯基距离;‘cosine’: ‘chebychev’:Chebychev距离。 linkage函数 调用格式:Z=linkage(Y,’method’) 说 明:用‘method’参数指定的算法计算系统聚类树。 Y:pdist函数返回的距离向量; method:可取值如下: ‘single’:最短距离法(默认); ‘complete’:最长距离法; ‘average’:未加权平均距离法; ‘weighted’: 加权平均法; ‘centroid’:质心距离法; ‘median’:加权质心距离法; ‘ward’:内平方距离法(最小方差算法)

  13. 练习题 根据信息基础设施的发展状况,对二十个国家的地区进行分类。

  14. 主成分分析 在实际问题中,研究多指标的问题是经常遇到的,然而在多数情况下,不同指标之间是有一定关系的。由于指标较多再加上指标之间有一定的相关性,势必增加了分析问题的复杂性。主成分分析就是设法将原来指标重新组合成一组新的互相无关的几个综合指标来代替原来指标,同时根据实际需要从中可取几个较少的综合指标尽可能多滴反映原来指标的信息。这种多个指标化为少数互不干扰的综合指标的统计方法叫做主成分分析法,如某人要做一件上衣要测量很多尺寸,如身长、袖长、胸围、腰围、肩宽、肩厚等十几项指标。但是某服装产生产一批新型服装绝不可能吧尺寸型号分的过多。而是从其中选取几个综合性的指标作为分类型号。1、反映胖瘦。2、反映特体。3反映长度。

  15. 计算步骤 设有n个样品,每个样品观测p个指标,将原始数据写成矩阵形式 特征值大的贡献大。 贡献率=特征值/所有特征值和 1、将原始数据标准化 2、建立变量的相关系数阵 3、求R的特征根及相应的单位特征向量a1,a2,.....ap 4、写出主成分 一般取累计贡献率达85—95%的特征值 所对应的第一、第二,…,第m(m≤p)个主成分。

  16. pcacov 功能:运用协方差矩阵进行主成分分析 格式:PC=pcacov(X) [PC,latent,explained]=pcacov(X) 说明:[PC,latent,explained]=pcacov(X)通过协方差矩阵X进行主成分分析,返回主成分(PC)、协方差矩阵X的特征值(latent)和每个特征向量表征在观测量总方差中所占的百分数(explained)。 例 中国大陆35个大城市某年的10项社会经济统计指标指标做主成分分析数据见下表。 相关系数矩阵: std = 1.0000 -0.3444 0.8425 0.3603 0.7390 0.6215 0.4039 0.4967 0.6761 0.4689 -0.3444 1.0000 -0.4750 0.3096-0.3539 0.1971 0.3571 0.2600 0.1570 0.3090 0.8425 -0.4750 1.0000 0.3358 0.5891 0.5056 0.3236 0.4456 0.5575 0.3742 0.3603 0.3096 0.3358 1.0000 0.1507 0.7664 0.9412 0.8480 0.7320 0.8614 0.7390 -0.3539 0.5891 0.1507 1.0000 0.4294 0.1971 0.3182 0.3893 0.2595 0.6215 0.1971 0.5056 0.7664 0.4294 1.0000 0.8316 0.8966 0.9302 0.9027 0.4039 0.3571 0.3236 0.9412 0.1971 0.8316 1.0000 0.9233 0.8376 0.9527 0.4967 0.2600 0.4456 0.8480 0.3182 0.8966 0.9233 1.0000 0.9201 0.9731 0.6761 0.1570 0.5575 0.7320 0.3893 0.9302 0.8376 0.9201 1.0000 0.9396 0.4689 0.3090 0.3742 0.8614 0.2595 0.9027 0.9527 0.9731 0.9396 1.0000

  17. 特征值(val) val = 0.0039 0 0 0 0 0 0 0 0 0 0 0.0240 0 0 0 0 0 0 0 0 0 0 0.0307 0 0 0 0 0 0 0 0 0 0 0.0991 0 0 0 0 0 0 0 0 0 0 0.1232 0 0 0 0 0 0 0 0 0 0 0.2566 0 0 0 0 0 0 0 0 0 0 0.3207 0 0 0 0 0 0 0 0 0 0 0.5300 0 0 0 0 0 0 0 0 0 0 2.3514 0 0 0 0 0 0 0 0 0 0 6.2602 特征根排序: 6.260222.351380.5300470.3206990.2566390.1232410.09909150.03070880.02403550.00393387

  18. 特征向量(vec): -0.1367 0.2282 -0.2628 0.1939 0.6371 -0.2163 0.3176 -0.1312 -0.4191 0.2758 -0.0329 -0.0217 0.0009 0.0446 -0.1447 -0.4437 0.4058 -0.5562 0.5487 0.0593 -0.0522 -0.0280 0.2040 -0.0492 -0.5472 -0.4225 0.3440 0.3188 -0.4438 0.2401 0.0067 -0.4176 -0.2856 -0.2389 0.1926 -0.4915 -0.4189 0.2726 0.2065 0.3403 0.0404 0.1408 0.0896 0.0380 -0.1969 -0.0437 -0.4888 -0.6789 -0.4405 0.1861 -0.0343 0.2360 0.0640 -0.8294 0.0377 0.2662 0.1356 -0.1290 0.0278 0.3782 0.2981 0.4739 0.5685 0.2358 0.1465 -0.1502 -0.2631 0.1245 0.2152 0.3644 0.1567 0.3464 -0.6485 0.2489 -0.4043 0.2058 -0.0704 0.0462 0.1214 0.3812 0.4879 -0.5707 0.1217 0.1761 0.0987 0.3550 0.3280 -0.0139 0.0071 0.3832 -0.7894 -0.1628 0.1925 0.2510 -0.0422 0.2694 0.0396 0.0456 0.1668 0.3799 于是的三个指标为: Y1=-0.1312*x1-0.5562*x2+0.3188*x3+....+0.0456*x10 Y2=-0.4191*x1+0.5487*x2+...........+0.1668*x10 Y3=0.2758*x1+ 0.0593*x2+............+0.3799*x10

  19. 通过观察我们发现Y1当中x2,x5的系数比较大,即影响Y1比较明显因此我们可将Y1看做反映非农业人口比与客运总量的综合指标。通过观察我们发现Y1当中x2,x5的系数比较大,即影响Y1比较明显因此我们可将Y1看做反映非农业人口比与客运总量的综合指标。 进一步还可做因子分析。 练习、我们给出了各地的企业的经济效益状况,通过相关的方法对各地的经济效益做分析。数据如下表:

  20. 相关性分析 在一元统计分析中,研究两入随机变量之间的线性相关关系、 可用相关系数(称为简单相关系数);研究一个随机变量与多个随 机变量之间的线性相关关系,可用复相关系数(称为全相关系 数)将它推广到研究多个随机变量与多个随机变量之间的相关关系的讨论中, 提出了典型相关分析。 实际问题中,两组变量之间具有相关关系的问题很多,例如 几种主要产品如猪肉、牛肉、鸡蛋的价格(作为第一组变量)和 相应这些产品的销售量(作为第二组变量)有相关关系;投资性 变量(如劳动各人数、货物周转量、生产建设投资等)与国民收 入变量(如工农业国民收入、运输业国民收入、建筑业国民收入 等)只有相关关系;患某种疾病的病人的各种症状程度(第一组 变量)和用物理化学方法检验的结果(第二组变量)具有相关关 系;运动员的体力测试指标(如反复横向跳、纵跳、背力、握力 等)与运动能力测试指标(如耐力跑、跳远、投球等)之间具有 相关关系等等。

  21. 典型相关分析就是研究两组变量之间相关关系的一种多元统计方法,设两组变量用x1,x2,…xn和y1,y2…yn表示,要研究两组变量的相关关系,一种方法是分别研究X和Y之间的相关关系,然后列出相关系数表进行分析,当两组变量较多时,这样做法不仅烦琐也不易抓住问题的实际;另一种方法采用典型相关分析就是研究两组变量之间相关关系的一种多元统计方法,设两组变量用x1,x2,…xn和y1,y2…yn表示,要研究两组变量的相关关系,一种方法是分别研究X和Y之间的相关关系,然后列出相关系数表进行分析,当两组变量较多时,这样做法不仅烦琐也不易抓住问题的实际;另一种方法采用 类似主成分分析的做法在每一组变量中都选择若干个有代表性的综合指标(变量的线性组合),通道研究两组的综合指标之间的关系来反映两组变量之间关系比如猪肉价格和牛肉价格用x1,X2表示,它们的销售售量用X,xl表示,研究它们之间的相又关系,从经济学观点就是希望构造一个X1、x2的线性函数入y1=a11X1十a12x2称为价格指数及x3、x4的线性函数y2=a21x3十a22X4称为销售指数,要求它们之间具有最大相关性,这就是一个典型相关分析问题。

  22. 1.插值拟合 2.线性回归 4.灰色分析 5.神经网络

  23. 在解决实际问题的生产(或工程)实践和科学实验过程中,通常需要通过研究某些变量之间的函数关系来帮助我们认识事物的内在规律和本质属性,而这些变量之间的未知函数关系又常常隐含在从试验、观测得到的一组数据之中。因此,能否根据一组试验观测数据找到变量之间相对准确的函数关系就成为解决实际问题的关键。在解决实际问题的生产(或工程)实践和科学实验过程中,通常需要通过研究某些变量之间的函数关系来帮助我们认识事物的内在规律和本质属性,而这些变量之间的未知函数关系又常常隐含在从试验、观测得到的一组数据之中。因此,能否根据一组试验观测数据找到变量之间相对准确的函数关系就成为解决实际问题的关键。 例如在工程实践和科学实验中,常常需要从一组试验观测数据(xi ,yi ) ,i = 0,1,....,n之中找到自变量x与因变量y 之间的函数关系,一般可用一个近似函数y = f (x)来表示。函数y = f (x)的产生办法因观测数据和要求不同而异,通常可采用数据拟合与函数插值两种办法来实现。

  24. 数据拟合主要是考虑到观测数据受随机观测误差的影响,进而寻求整体误差最小、能较好反映观测数据的近似函数y = f (x),此时并不要求所得到的近似函数y = f (x)满足yi= f (xi) , i = 0,1,…,n。 函数插值则要求近似函数y = f (x)在每一个观测点i x 处一定要满足y i= f (xi) , i = 0,1,…,n,在这种情况下,通常要求观测数据相对比较准确,即不考虑观测误差的影响。

  25. 在实际问题中,通过观测数据能否正确揭示某些变量之间的关系,进而正确认识事物的内在规律与本质属性,往往取决于两方面因素。其一是观测数据的准确性或准确程度,这是因为在获取观测数据的过程中一般存在随机测量误差,导致所讨论的变量成为随机变量。其二是对观测数据处理方法的选择,即到底是采用插值方法还是用拟合方法,插值方法之中、拟合方法之中又选用哪一种插值或拟合技巧来处理观测数据。插值问题忽略了观测误差的影响,而拟合问题则考虑了观测误差的影响。但由于观测数据客观上总是存在观测误差,而拟合函数大多数情况下是通过经验公式获得的,因此要正确揭示事物的内在规律,往往需要对大量的观测数据进行分析,尤为重要的是进行统计分析。统计分析的方法有许多,如方差分析、回归分析等。在实际问题中,通过观测数据能否正确揭示某些变量之间的关系,进而正确认识事物的内在规律与本质属性,往往取决于两方面因素。其一是观测数据的准确性或准确程度,这是因为在获取观测数据的过程中一般存在随机测量误差,导致所讨论的变量成为随机变量。其二是对观测数据处理方法的选择,即到底是采用插值方法还是用拟合方法,插值方法之中、拟合方法之中又选用哪一种插值或拟合技巧来处理观测数据。插值问题忽略了观测误差的影响,而拟合问题则考虑了观测误差的影响。但由于观测数据客观上总是存在观测误差,而拟合函数大多数情况下是通过经验公式获得的,因此要正确揭示事物的内在规律,往往需要对大量的观测数据进行分析,尤为重要的是进行统计分析。统计分析的方法有许多,如方差分析、回归分析等。

  26. 数据拟合虽然较有效地克服了随机观测误差的影响,但从数理统计的角度看,根据一个样本计算出来的拟合函数(系数),只是拟合问题的一个点估计,还不能完全说明其整体性质。因此,还应该对拟合函数作区间估计或假设检验,如果置信区间太大或包含零点,则由计算得到的拟合函数系数的估计值就毫无意义。这里所采用的统计分析方法就是所谓的回归分析。另外还可用方差分析的方法对模型的误差作定量分析。数据拟合虽然较有效地克服了随机观测误差的影响,但从数理统计的角度看,根据一个样本计算出来的拟合函数(系数),只是拟合问题的一个点估计,还不能完全说明其整体性质。因此,还应该对拟合函数作区间估计或假设检验,如果置信区间太大或包含零点,则由计算得到的拟合函数系数的估计值就毫无意义。这里所采用的统计分析方法就是所谓的回归分析。另外还可用方差分析的方法对模型的误差作定量分析。 对于插值方法,本章简单介绍最常用的插值法的基本结论及其Matlab实现问题。由于数据拟合问题必须作区间估计或假设检验,所以除了在本章介绍最基本的数据拟合方法——最小二乘法的基本结论及其Matlab实现问题外,我们在专门介绍了对数值拟合问题进行区间估计或假设检验的统计方法,

  27. 即介绍回归分析方法及其Matlab实现。 数据处理问题通常情况下只是某个复杂实际问题的一个方面或部分内容,因而这里所介绍的数据处理方法——函数插值和数据拟合的方法(包括回归分析)通常只能解决实际问题中的部分问题——计算问题。一般来说,对实际问题进行数学建模需要用到多方面知识,只有很少的情况下可以单独使用本章所介绍的内容,故我们只在本章最后一节以修改后的美国91年数学建模A题为例说明如何使用数值计算知识建立数学模型,从而解决实际问题的方法。

  28. 插值方法 1、拉格朗日插值法 分段线性插值的Matlab实现 用Matlab实现分段线性插值不需要编制函数程序,Matlab中有现成的一维插值函数interp1。 y=interp1(x0,y0,x,'method') method指定插值的方法,默认为线性插值。其值可为: 'nearest' 最近项插值 'linear' 线性插值 'spline' 立方样条插值 'cubic' 立方插值。 2、分段线性插值法

  29. 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是插值点的函数值。 3、三次样条插值法 对于三次样条插值,我们提倡使用函数csape,csape的返回值是pp形式,要求插值点的函数值,必须调用函数ppval。

  30. 例1 机床加工 待加工零件的外形根据工艺要求由一组数据(x, y)给出(在平面情况下),用程控铣床加工时每一刀只能沿x方向和y 方向走非常小的一步,这就需要从已知数据得到加工所要求的步长很小的(x, y)坐标。 表中给出的x, y数据位于机翼断面的下轮廓线上,假设需要得到x坐标每改变0.1时的y坐标。试完成加工所需数据,画出曲线,并求出x = 0处的曲线斜率和13 ≤x ≤ 15范围内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 要求用分段线性和三次样条two种插值方法计算。

  31. 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,’linear’); y3=interp1(x0,y0,x,'spline'); pp1=csape(x0,y0); y4=ppval(pp1,x); pp2=csape(x0,y0,'second'); y5=ppval(pp2,x); subplot(2,2,2) plot(x0,y0,'+',x,y2) title('Piecewise linear') subplot(2,2,3) plot(x0,y0,'+',x,y3) title('Spline1') subplot(2,2,4) 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] 计算结果略。 可以看出,分段线性插值的光滑性较差(特别是在x = 14附近弯曲处),建议选用三次样条插值的结果。

  32. 五 一维插值总结 插值函数一般是已知函数的线性组合或者称为加权平均。在已知数据点较少时,插值技术在工程实践和科学实验中有着广泛而又十分重要的应用。例如在信息技术中的图像重建、图像放大过程中为避免图像失真、扭曲而增加的插值补点,建筑工程的外观设计,化学工程试验数据与模型分析,天文观测数据、地理信息数据的处理,社会经济现象的统计分析等方面,插值技术的应用是不可或缺的。 插值技术(或方法)远不止这里所介绍的这些,但在解决实际问题时,对于一维插值问题而言,前面介绍的插值方法已经足够了。剩下的问题关键在于什么情况下使用、怎样使用和使用何种插值方法的选择上。 拉格朗日插值函数在整个插值区间上有统一的解析表达式,其形式关于节点对称,光滑性好。但缺点同样明显,这主要体现在高次插值收敛性差(龙格现象);增加节点时前期计算作废,导致计算量大;一个节点函数值的微小变化(观测误差存在)将导致整个区间上插值函数都发生改变,因而稳定性差等几个方面。因此拉格朗日插值法多用于理论分析,在采用拉格朗日插值方法进行插值计算时通常选取n < 7。分段线性插值函数(仅连续)与三次样条插值函数(二阶导数连续)虽然光滑性差,但他们都克服了拉格朗日插值函数的缺点,不仅收敛性、稳定性强,而且方法简单实用,计算量小。因而应用十分广泛。

  33. 数据拟合 在科学计算中经常要建立实验数据的数学模型。给定函数的实验数据,需要用比较简单和合适的函数来逼近(或拟合)实验数据。这种逼近的特点是: (a) 适度的精度是需要的; (b) 实验数据有小的误差; (c) 对于某些问题,可能有某些特殊的信息能够用来选择实验数据的数学模型。 逼近离散数据的基本方法就是曲线拟合,常采用最小二乘拟合 曲线拟合问题的数学描述是,已知一组(二维)数据(xi,yi ) ,i = 1,2,。。。,n(即平面上的n个点(xi, yi ) ,i = 1,2,。。,n), x i 互不相同。寻求一个函数(曲线) y = f (x),使f (x)在某种准则下与所有数据点最为接近,即曲线拟合得最好。 最小二乘拟合分为线性最小二乘拟合和非线性最小二乘拟合。

  34. 一 、线性最小二乘拟合 线性最小二乘法是解决曲线拟合问题最常用的方法,基本思路是,令 f(x)=a1*r1(x)+a2*r2(x)+…….+am*rm(x) 其中rk(x)是一组事先选定的线性无关的函数,(ak)是一组待定系数。寻求系数(ak)使得yi与f(xi)的距离d i(i = 1,2,,n)的平方和最小。这种准则称为最小二乘准则,其求系数的方法称为线性最小二乘拟合方法。 1、 系数(ak)的求法 若记 则J 为a1‥am的二次函数。由高等数学的极值理论,J 达到最小的充分必要条件是a1‥am满足dJ/dak =0(k = 1,。。。,m)。于是得到求使J 达到最小的a1‥am的方法是求解线性方程组(称为法方程组)

  35. 即求解线性方程组 若记 则以上方程组可表示为 RT RA = RTY 。 由于当{ r1(x ), ‥ , rm( x)} 线性无关时,R列满秩,RT R可逆,所以方程组(10.7)有唯一解 A = (RT R)-1RTY 。

  36. 用以上方法作线性最小二乘拟合的误差通常考虑以下两种形式:用以上方法作线性最小二乘拟合的误差通常考虑以下两种形式: 最小平方误差: 最大偏差: 2、 函数组{ r1(x ), ‥ , rm( x)}的选取 面对一组数据(x i , y i ), i = 1,2,。。。,n ,用线性最小二乘法作曲线拟合时,首要的、也是关键的一步是恰当地选取{ r1(x ), ‥ , rm( x)} 。

  37. 如果通过机理分析,能够知道y 与x之间应该有什么样的函数关系,则{ r1(x ), ‥ , rm( x)}容易确定。若无法知道y 与x之间的关系,通常可以将数据(x i , y i ), i = 1,2,。。,n,作图,直观地判断应该用什么样的曲线去作拟合。人们常用的曲线有: (i)直线 y = a1 x + a2; (ii)多项式 y = a1 x m + ‥ +amx+am+1, (一般m = 2,3,不宜太高) (iii)双曲线(一支)y=a1/x+a2拟合前作变量替换t = 1/x求解 a 1,a2较简单 (iv)指数曲线 y =a1ea2x 拟合前作变量代换z = ln y,t = 1/x,则指数曲线y =a1ea2x 转化为关于lna1,a2的线性函数z = ln a1+ a2t,这样做求 解 a1,a2 较简单。

  38. 在实际计算过程中,面对一组已知数据,到底用什么样的曲线拟合最好,可以在直观判断的基础上,选几种曲线分别拟合,然后比较,看哪条曲线的最小二乘指标J 最小。 二 、非线性最小二乘拟合 非线性最小二乘法是假设f (x)是待定系数{ak}的任意非线性函数,在最小二乘准则下求其系数{ak} 。 例如上述人们常用的双曲线和指数曲线就是非线性最小二乘拟合中最常用的非线性函数,只不过在上面使用中将它们转变成线性最小二乘拟合方法。 对于给定的实验数据,通常应根据实验数据的走向、趋势选择合适的数学模型,即拟合函数。例如当实验数据具有单调性和凸性时,可选择下述适当的数学模型y = f (x)来拟合实验数据。 f (x) = aebx,f (x) = aeb/x, f (x) = axb, f (x) = a + b / x

  39. 在有可能的情况下,一般将非线性拟合函数转化为线性拟合函数求解,这一方面是如此求解简单,另一方面也是因为一般情况下求解法方程组dJ/dak=0得到的(a1, ‥, am) 通常仅是J 的驻点,不一定是极值点。也可以直接解J 极小化问题。 三、 最小二乘拟合法的Matlab实现 命令为 A = R \Y 例3 用最小二乘法求一个形如y = a1 + bx2的经验公式,使它与下表所示的数据拟合. x 19 25 31 38 44 y 19.0 32.3 49.0 73.3 97.8 解 编写程序如下 x=[19 25 31 38 44]'; y=[19.0 32.3 49.0 73.3 97.8]'; r=[ones(5,1),x.^2]; ab=inv(r’*r)*r’*y x0=19:0.1:44; y0=ab(1)+ab(2)*x0.^2; plot(x,y,'o',x0,y0,'r')

  40. 2 线性最小二乘拟合(多项式拟合)方法 在线性最小二乘拟合中,用的较多的是多项式拟合。如果取 { r1( x), ‥, rm+1( x)} ={1, ‥ ,xm } ,即用m 次多项式拟合给定数据,则Matlab中有现成的函数 a=polyfit(x0,y0,m), 其中输入参数x0,y0为要拟合的数据,m为拟合多项式的次数,输出参数a为拟合多项式 y=amxm+…+a1x+a0系数a=[ am, …, a1, a0]。 多项式在x处的值y可用下面的函数计算 y=polyval(a,x)。 例4 某乡镇企业1990-1996年的生产利润如下表: 年份 1990 1991 1992 1993 1994 1995 1996 利润(万元) 70 122 144 152 174 196 202 试预测1997年和1998年的利润。

  41. 解 作已知数据的的散点图, x0=[1990 1991 1992 1993 1994 1995 1996]; y0=[70 122 144 152 174 196 202]; plot(x0,y0,'*') 发现该乡镇企业的年生产利润几乎直线上升。因此,我们可以用y = a1 x + a0作为拟合函 数来预测该乡镇企业未来的年利润。编写程序如下: x0=[1990 1991 1992 1993 1994 1995 1996]; y0=[70 122 144 152 174 196 202]; a=polyfit(x0,y0,1) y97=polyval(a,1997) y98=polyval(a,1998) 求得a1 = 20 ,a0= -4.0705×104,1997年的生产利润y97=233.4286,1998年的生产利润y98=253.9286。

  42. 3 非线性最小二乘拟合 Matlab的优化工具箱中提供了两个求非线性最小二乘拟合的函数:curvefit和leastsq。使用这两个命令时,都要先建立M文件fun.m,但它们定义f (x)的方式是不同的。 1 curvefit 设已知xdata=(xdata1,xdata2,…,xdatan ),ydata=(ydata1,ydata2,…,ydatan ), curvefit用以求含参量x(向量)的向量值函数F(x,xdata)=(f(x,data1), …,f(x,xdata n )) T中的参变量x(向量),使得 Sum(F(x,xdatai)-ydatai)2最小 输入格式为: (1)x=curvefit('fun',x0,xdata,ydata); (2)x=curvefit('fun',x0,xdata,ydata,options); (3)x=curvefit('fun',x0,xdata,ydata,options, 'grad'); (4)[x,options]=curvefit('fun',x0,xdata,ydata,…); (5)[x,options,funval]=curvefit('fun',x0,xdata,ydata,…); (6)[x,options,funval,Jacob]=curvefit('fun',x0,xdata,ydata,…). 输出目标函数值格式:f=fun(x,xdata). 其中x0为迭代初值,options为控制参数。

  43. 2 leastsq 设已xdata=(xdata1,xdata2,…,xdatan ),ydata=(ydata1,ydata2,…,ydatan ), leastsq 用以求含参量x(向量)的向量值函数 输入格式为: (1)x= leastsq ('fun',x0); (2)x= leastsq ('fun',x0,options); (3)x= leastsq ('fun',x0,options, 'grad'); (4)[x,options]= leastsq ('fun',x0,…); (5)[x,options,funval]= leastsq ('fun',x0,…); 例5 用下面一组数据拟合函数 c(t) = a + be-0.02kt中的参数a,b, k 。 t 100 200 300 400 500 600 700 800 900 1000 cj×103 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59

  44. 1 用命令curvefit。此时 F(x,tdata)=(a+b e-0.02kt1,…,a+be-0.02kt10)T,x=(a,b,k) (1) 编写M文件curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中x(1)=a;x(2)=b;x(3)=k; (2) 输入命令 tdata=100:100:1000 cdata=1e03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]; x0=[0.2,0.05,0.005]; x=curvefit(‘curvefun1’,x0,tdata,cdata) f=curvefun1(x,tdata) (3)运算结果为: x=0.0070 -0.0030 0.1012 f= Columns 1 through 7 0.0045 0.0050 0.0054 0.0057 0.0059 0.0061 0.0063 Columns 8 through 10 0.0064 0.0065 0.0066 即拟合得a=0.0070,b=-0.0030,k=0.0066

  45. 2 用命令leastsq。此时 f(x)=F(x,tdata,cdata)=(a+be-0.02kt1-c1,…,a+be-0.02kt10-c10)T,x=(a,b,k) (1) 编写M文件curvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]; f=cdata-x(1)-x(2)*exp(-0.02*x(3)*tdata) %其中x(1)=a;x(2)=b;x(3)=k; (2) 输入命令 x0=[0.2,0.05,0.005]; x=leastsq(‘curvefun2’,x0) f=curvefun2(x) (3)运算结果为: x=0.0070 -0.0030 0.1012 f=1.0e-005* Columns 1 through 7 0.0221 0.2081 -0.3933 -0.2872 0.2973 0.3561 0.0693 Columns 8 through 10 -0.2327 -0.0970 0.0296 可以看出,两个命令的计算结果是相同的

  46. 回归分析 回归分析是处理很难用一种精确方法表示出来的变量之间关系的一种数学方法,它是最常用的数理统计方法,能解决预测、控制、生产工艺优化等问题。它在工农业生产和科学研究各个领域中均有广泛的应用。回归分析一般分为线性回归分析和非线性回归分析。本节着重介绍线性回归分析的基本结论及其在Matlab中的相应命令。线性回归分析是两类回归分析中较简单的一类,也是应用较多的一类。 一 、一元线性回归分析 针对一组(二维)数据( xi,yi ),i = 1,2,。。,n(其中xi 互不相同),其最简单的数据拟合形式为寻求直线y = b1 + b2*x ,使b1 + b2*x在最小二乘准则下与所有数据点最为接近。但由于随机观测误差的存在,满足上述数据点的直线应该是 y = b1 + b2* x +e, (1.1) 其中x, y是准确的,b1 ,b2 是两个未知参数,e 是均值为零的随机观测误差,具有不可观测性,可以合理地假设这种观测误差服从正态分布。于是我们得到一元线性回归模型为:

  47. y = b1 + b2* x +e 正态分布 E(e) =0,D(e)=s2 其中s 未知,固定的未知参数 b1 、 b2 称为回归系数,自变量x称为回归变量。 式两边同时取期望得: Y = b1 + b2*x ,称为y 对x的回归直线方程。 在该模型下,第i个观测值可以看作样本Yi = b1 + b2*xi +ei (这些样本相互独立但不同分布, i = 1,2,,n)的实际抽样值,即样本值。 一元线性回归分析的主要任务是:用实验值(样本值)对b1 、b2和s 作点估计;对回归系数b1 、b2 作假设检验;在x = x0 处对y 作预测,并对y 作区间估计。

  48. 2 模型的假设、预测、控制 1回归方程的显著性检验 在实际问题中,因变量y 与自变量x之间是否有线性关系(1.1)只是一种假设,在求出回归方程之后,还必须对这种回归方程同实际观测数据拟合的效果进行检验。由(10.10)可知,| b2 |越大, y 随x变化的趋势就越明显;反之,| b2 |越小, y 随x变化的趋势就越不明显。特别当b2 =0时,则认为y 与x之间不存在线性关系,当 b2 ≠0时,则认为y与x之间有线性关(1.1)。因此,问题归结为对假设 H0:b2=0;H1:b2 ≠0 进行检验。假设: H0:b2=0被拒绝,则回归显著,认为y 与x之间存在线性关系,所求的线性回归方程有意义;否则回归不显著, y 与x的关系不能用一元线性回归模型来描述,所得的回归方程也无意义。此时,可能有如下几种情况: (i)x对y 没有显著影响,此时应丢掉变量x; (ii)x对y 有显著影响,但这种影响不能用线性关系来表示,应该用非线性 回归; (iii)除x之外,还有其他不可忽略的变量对y 有显著影响,从而削弱了x对 y 的影响。 此时应用多元线性回归模型。因此,在接受H 0的同时,需要进一步查明原因以便分别处理。

More Related