基于MATLAB的地学数据处理与统计分析探讨(共27页)

1、021111班-彭宁俊-20111000713自主学习报告报告(bogo)抬头课题(kt)名称：基于(jy)MATLAB的地学数据处理与统计分析探讨课题组长：姓名 彭宁俊 学号 20111000713 班号 021111 课题成员：（小组协作完成需填写）成员分工：（文中所有MATLAB语句、绘图操作和报告内容系笔者自主完成，实习中提供过部分数据和资料来源的同学已在致谢栏中给出，所有原始数据和图件在附件中给出） 指导教师： 陈志军老师 （可结合你们的毕业实习指导教师的相关研究）起止时间：2014/4/19-2014/4/29目录：0摘要1MATLAB简介2.MATLAB地学数据处理 2.1 MA

2、TLAB数据统计量计算 2.1 MATLAB数据插值和拟合 以Exercises.1-2部分数据、peaks函数为例 2.3 概率函数和分布函数 2.4 参数估计和假设检验 以Exercises.1-1数据为例3MATLAB数据可视化和绘图 3.1数据分布与直方图 以Exercises.1-1数据为例3.2分布检验函数图 以Exercises.1-2、1-3数据为例3.3常见统计图绘制以上届实习-聚类分析原始数据为例3.4其他统计图绘制3.5特殊图形绘制 3.5.1等值线图 3.5.2三维曲线、网格、曲面图4.MATLAB地学统计模型和统计分析 4.1方差分析以某工程师对不同类型金矿预测得率数

3、据为例 4.2回归分析以资勘2班-实习3回归分析原数据为例 4.2.1多元线性回归4.2.2逐步回归 4.3主成分分析以某地土壤(trng)不同类型组成成分数据为例 4.4聚类分析以某地土壤不同类型组成成分(chng fn)数据为例 4.5判别分析以资源1班-实习(shx)3判别分析原数据为例 4.6自相关分析以某海域测得的年平均海平面高度数据为例5.结语与感想6问题探讨7致谢8参考文献9 附件基于MATLAB的地学数据处理与统计分析探讨0摘要：数学在任何领域的应用都主要发挥两种功能：一是作为实验或者观察数据的整理手段；二是构造假设、建立模型和发展理论。无论哪种功能的发挥，都涉及数据的处理和分

4、析。MATLAB很好地利用了这两种功能，其在数学地质领域的作用自不待言。笔者就是在掌握MATLAB基本操作的基础上，以具体地学数据为出发点，以数学方法为支撑点探讨MATLAB在地学数据处理以及统计分析中的具体应用。MATLAB简介 MATLAB是MATrix LABoratory(/矩阵实验室的缩写，是由美国MathWorks公司开发的集数值计算、符号计算和图形可视化三大基本功能于一体、功能强大、操作简单的语言，是国际公认的优秀数学应用软件之一。强大的科技应用软件和编程语言， 易学易用，应用十分广泛，以备受广大学子的青睐。在数学地质领域，MATLAB基础数据处理显得轻而易举。我们还可以灵活的利

5、用这些MATLAB多变量统计分析及类似方法对给定的特定的地质对象进行定性描述、识别、预测、成因研究；根据系统的数学模型借助强大的数值计算的功能对地质过程进行数值模拟；通过各种地质变量的空间分布运用MATLAB随机过程的理论与方法进行研究，以探索空间数据的变化规律并进行预测等等。在矿床统计预测的数次实习中我们对MATLAB早已屡见不鲜，如实习1.混合总体筛分：复制matlab code.txt代码计算各分组的累计概率；运行MF4SPOD_MixDist_Preproc.m文件调用SPOD_mixdist_TFe_Raw中的数据进行特异值剔除，显示剔除特异值后的数据并绘制出直方图；实习2.证据权法

6、：把4处缺失的语句补充完整即可获得后验概率和系列结果；实习3.聚类分析:导入需进行聚类分析的数据文件： ht_ores_TextID，运行m文件MPHX_HierarchicalCluster，得到方阵 v.pdist 、 x.pdist并绘制谱系图；实习四.趋势面分析:导入数据后调用函数trend234.m文件，算出一系列的统计分析值，运行脚本文件trendFig.m绘出等值线图、趋势面图、剩余图等等。这些都是老师编好的文件，只需稍微修改就可运行，快捷准确、极大的提高了计算的效率。受到课堂(ktng)启发，课下探讨学会(xuhu)了MATLAB的基本操作，研究了其在多方面的应用，并总结出了一

7、些规律。此次探索性学习(xux)笔者就是在这基础上通过探讨MATLAB重点在 = 1 * GB3 基础数据处理 = 2 * GB3 数据可视化及统计绘图 = 3 * GB3 统计模型与统计分析等三方面的具体应用，切实以地学数据为基础、数学和计算机为手段研究客观地质问题。2.MATLAB与地学数据处理 数值计算是MATLAB最强大的功能之一，其矩阵计算能力强、编程简单，计算代码通过精心优化，计算效率高。在数值计算方面，涉及的数学学科面广，例如矩阵代数计算、微积分计算、快速傅立叶变换、常微分方程计算、偏微分方程计算、数据插值和拟合等，MATLAB提供了非常丰富的解决这类问题的命令。地学研究设计众多

8、的工程实验和工程测量，常常需对数据进行离散数据分析处理，找出测量数据的数学规律。数据插值与拟合、数理统计分析都属于这个范畴。本章笔者将采用MATLAB中一些常用的数理统计方法和技巧分析地学数据可能存在的数理规律，认识基础基础统计量、熟悉基本统计方法、探讨数据缺失与异常以及插值和拟合等系列基础问题。2.1 MATLAB数据统计量计算我们常需对我们的实验数据进行统计分析，统计量的计算必不可少。常见的统计量有表示位置的统计量、表示变异程度的统计量、表示分布形状的统计量等等，MATLAB统计量计算命令常含：最大值Max(),最小值min(),中位数median(),平均值Mean(),按指定维数求平均

9、值mean(A,n),算术平均nanmean(),几何平均geomean(),和谐平均harmmean();调整平均trimmean();方差var();标准差std(),协方差Cov(),相关系数coffcoef();由小到大排序sort();极差range();绝对差分平均值mad();按首行排序sortrows();求和sum();求元素此前的元素和cumsum();求梯形累和cumtrapz(),求元素此前所有元素的积cumprod();这些语句简单明了、快速便捷，所得的结果可直接显示是其他非专业软件难以比拟的。2.1 MATLAB插值和拟合MATLAB 软件以其强大的数据分析处理功能

10、以及简洁易懂的编程语言， 在很多领域都有了应用，把试验数据绘制成图是点状分布，将这些点练成曲线，然后转化为有意义的数学函数，才能对数据做对比和分析，这一过程过程涉及插值(interpolation)和曲线拟合(curve-fitting)。插值过程中认为数据(shj)是准确的，依据采样点之间的几何关系和关联信息(xnx)进行补值或求取(qi q)描述点之间的数据， ，MATLAB一维插值函数为interp1(a，c，b)和interp1(a，c，b，method)。a，c为已知数据，而b为要插值的数据点，method为预先设定的插值方法，分别为线性(linear)、三次多项式(cubic)和s

11、pline(样条插值)；曲线拟合中，假定已知曲线的规律，做曲线的最佳逼近，不需要经过所有数据点，如果拟合的曲线限定为多项式就称为多项式最小二乘法曲线拟合。假设假设以下数据为在某测线上等距性采集的11个火山样品铁元素含量（单位：wt%）;5.19，5.25 ，3.87，4.68 ，5.06 ，5.48，4.39，5.17，4.54，4.61，4.81，通过MATLAB-basic-fiting可以方便快捷的绘出不同类型的插值和拟合曲线，不用输入公式便可直接插值或拟合，并可显示多项式回归方程和残差图。根据残差分布可以判断拟合程度。图1.Basic Fiting数据插值和拟合图 高维插值和一维插值原

12、理一样，以二维插值为例，其语法为：ZI=interp2（X,Y,Z,XI,YI,method）,其中的参数和一维的类似，为方便理解和运算，笔者以MATLAB自带的peaks函数为例，首先取较少的数据，然后再求插值函数比较不同方法的区别： clear x,y=meshgrid(-3:1:3); z=peaks(x,y); subplot(2,2,1),surf(x,y,z);title(较少点的表面(biomin)图) axis(-3,3,-3,3,-5,5) xi,yi=meshgrid(-3:0.25:3); zi1=interp2(x,y,z,xi,yi,nearest); subplot

13、(2,2,2),surf(xi,yi,zi1);title(nearest法插值结果(ji gu) axis(-3,3,-3,3,-5,5) zi2=interp2(x,y,z,xi,yi,bilinear); subplot(2,2,3),surf(xi,yi,zi2);title(bilinear法插值结果(ji gu) axis(-3,3,-3,3,-5,5) zi3=interp2(x,y,z,xi,yi,bicubict); subplot(2,2,4),surf(xi,yi,zi3);title(bicubict法插值结果) axis(-3,3,-3,3,-5,5) xlabel(

14、绘制：彭宁俊)图2.Peaks函数二维插值结果图二维插值在地学领域应用也是非常广泛的，比如给出了一系列采样点的X,Y坐标和相应点的的采样元素含量，通过插值可以推算某些未采样点的元素含量值。此外缺失数据的处理是一个非常困难的问题，随根据具体问题的不同，处理方法也各异。为了数据分析的方便，将缺失数据用NaN表示是一个非常便利的方法。对异常数据可以采用与缺失数据相似的处理方法，即去除异常数据。至于异常数据的标准，将视具体问题而定，实际中经常使用的一种标准是：与平均值的偏差大于3倍标准差。2.2密度(md)函数和分布函数 常用(chn yn)的概率密度函数含normpdf()、chi2pdf()、tp

15、df()、fpdf()，通用(tngyng)概率密度公式为Y=pdf(name,x,a,b,c),其中MATLAB支持的name种类含beta、bino、chi2、exp、logn、nbin、ncf、nct、ncxz、norm、poiss、rayl、t、unif、unid、wbl. a、b、c为常数根据具体函数类型选择个数。常用的分布函数有normcdf()、chi2cdf()、tcdf()、fcdf(),通用分布函数为cdf(name,k,a,b,c) ,a、b、c为常数根据具体函数类型选择个数，函数表示的是xK时的概率值之和。此外应用随机数发生器也是MATLAB经常用到的操作，通用格式为Y

16、=random（name，A,B,C,M,N。具体的函数详细应用方法请参照相关资料，这里不详尽叙述，以下为案例：x=-4:0.1:4;p1=normpdf(x,0,0.6);p2=normpdf(x,0,0.1); p3=normpdf(x,-1,0.25);p4=normpdf(x,1,0.4); subplot(2,2,1),plot(x,p1,x,p2,x,p3,x,p4) p5=fpdf(x2,5,10);p6=fpdf(x2,10,10); p7=fpdf(x2,10,5); subplot(2,2,2),plot(x2,p5,x2,p6,x2,p7)%图3、图4类似，图标在Plot

17、 Editor中编辑图3.MATLAB概率密度图和函数分布图2.3参数估计和假设检验 参数估计和假设检验常常(chngchng)是统计假设需运到的问题。参数估计就是用某一函数值作为总体未知参数的估计值，又称为点估计。点估计分为矩估计和最大释然估计。另一类估计为区间(q jin)估计，是对于未知参数给一定范围在一定可靠度下求其真实值。如正态分布的参数估计语句(yj)为：mu,sigma,muci,sigmaci=normfit(x,alpha),mu、sigma分别为均值和标准差，muci和sigmaci为置信区间，为样本矩阵，alpha为显著性水平。假设检验即在总体分布未知或虽知其类型但含未知

18、参数时，提出关于总体的假设，对所做出的假设作出接受还是拒绝的推断，假设检验常遇到的有单个正太总体均值的假设检验；单个正太总体方差的假设检验；两个正太总体均值的假设检验。如单个正太整体的假设检验语句为：h,p,ci=ztest(x,mu,sigma,alpha,tail)(方差已知)。下面以SPOD自主学习网站研究性学习文件夹中的Exercise.1-1.txt数据（The Fe content (in wt%) of 120 volcanic samples），注意由于排版和数据量大，文章中设计的原始数据均在附件中已给出 ，读者自行查阅： Exercise.1-1.txt原始资料为采自120个

19、火山样品的2616个Fe含量百分比数据，如果 = 1 * GB3 已知这批样品服从正太分布，求在0.95置信水平下样品均值和方差的最大释然估计值？ = 2 * GB3 已知样品服从正太分布，方差未知，检验其均值取51是否可靠？ = 1 * GB3 参数估计：导入数据后更名为a,输入如下语句 mu,sigma,muci,sigmaci=normfit(a,0.95) %参数估计即得mu =50.9108 sigma = 10.2591 muci = 50.8982 50.9234 sigmaci = 10.2515 10.2693在置信水平为0.95时，均值和标准差的最大释然估计分别为50.91

20、和10.26，均值和标准差的知心区间为50.8982 ,50.9234 和10.2515 ,10.2693。 = 2 * GB3 假设检验：导入数据后更名为a,输入如下语句 h,p,ci=ttest(a,60.2) %假设检验即得：h=0 p= 0.6566 ci =50.5175 51.3041 H=0表示不拒绝原假设，且p=0.65660.5概率较大表示接受，即认为均值取51是合理的，且位于ci =50.5175 51.3041之间精度较高。3. MATLAB数据可视化及绘图 作图功能是MATLAB的特长之一，MATLAB提供了高端和低端的作图命令，我们可方便的将数据绘成二维、三维的各种图

21、形。作图分为两大类，一类是数据作图，如给定数组x，y为数据的图形，即数据可视化作图，也是科学研究和生产应用频率最高的作图方法。另一类是给定一个函数，如作出X2+Y2的图形。本部分笔者通过实践主要总结出MATLAB如下类型的图形绘制方法，前四部分可以认为是统计图形制作（2.1-2.3例子均以我们实习接触的数据为例，2.4-2.5重在作图方法的理解）。此外2.5还具有美学欣赏价值，重在兴趣激发。3.1 数据分布与直方图（Data distribution and histograms） 以矿床统计预测(yc)网站研究性学习文件夹中的Exercise.1_1.doc数据(shj)（a dataset

22、 of SiO2 content (wt%) in river sediments from the French Massif Central）为例，首先(shuxin)新建1_1.txt文本文档将数据复制进去，输入如下语句即可快速绘图： a=load(C:UsersAdministratorDesktopMATLAB practice1-1.txt)b=0:10:130;subplot(2，2，1)，hist(a，b) xlabl( x 10% per class) ylabel(y values ) title(直方图1)% 以直方图1绘制为例，其他类似 subplot(2，2，4)，h

23、istfit(a) %拟合图语句 title(拟合图) 图4.数据分布与直方图绘制 由上可以看出每组x不同间隔图形的形状类似，均为中间高峰突起，两边逐渐下降的曲线，跟正态分布很类似，最后作出的拟合图更加形象，曲线与直方图拟合的很好也有利地证明了原始数据的分布特征。3.2 MATLAB分布检验函数图 分布检验我们以正态分布检验和对数正太分布检验为基础。正态分布检验以老师给定的Exercise.1-1.txt数据（The Fe content (in wt%) of 120 volcanic samples）， 对数正太分布分布检验以Exercise.1-2.txt 数据(The Sr conte

24、nt of 60 volcanic rocks from the Banda Arc in Indonesia)为例。值得一提的是由于Exeicise.1-2txt给出的数据为120个，我们在最后一个添加一个0，构成1111的矩阵，数据导入之后经b=a(:)变换再截取c=b(1:120,1)成1201维矩阵，既不改变原数据又可方便MATLAB运算。MATLAB里面含有大量(dling)的分布检验函数，如normplot、cdfplot()，。我们也可以通过分布函数(hnsh)产生随机数如x=normrnd()，y=chi2rnd()，对于(duy)初学者可以通过随机数并进行分布验证可以产生理想

25、的效果。图5. Exeicise.1-2正态分布检验系列图The Fe content (in wt%) of 120 volcanic samples在算术概率格子图上所作出图为直线，说明其满足正态分布；所作出的经验分布函数如图2所示曲线，直方图曲线拟合曲线，均明显显示正态分布的特征；最后一图为Boxplot。图6. Exeicise.1-3对数正太分布(fnb)分布检验 The Sr content of 60 volcanic rocks from the Banda Arc in Indonesia数据直方图显示平均值大于众数，表现明显正偏，正太函数分布检验为曲线；而取对数(du sh

26、)后直方图对称分布，检验函数表现为一条直线，说明原数据是呈对数正态分布的。3.3常见(chn jin)统计图绘制（plot、scatter、stairs、stem、bar、pie） 统计图形在地学领域应用的十分广泛，MATLAB绘制常见的统计图形包括散点图scatter（X,Y,filled，S）、阶梯图stairs(X,Y,S)、离散杆状图stem(X,Y,filled，S)、条形图bar(X,Y,WIDTH,参数)、饼状图pie（X,explode）等等；其语法都类似，易学易懂，这里只做简要介绍，读者可以慢慢探索。下面以10级矿床统计预测-聚类分析实习课程中的原始数据（见附件）为例，数据为

27、14（样品）11（元素）矩阵：图7.MATLAB常见(chn jin)统计图绘制以上统计(tngj)图形对于MATLAB来说可谓得心应手(d xn yng shu)，一键植入即可马上出结果。需要注意的是条形图Bar（）和前面的直方图Hist（）在原理和用法上是有区别的：直方图往往为一维的，以区间范围为横坐标、以频数为纵坐标；而条形图bar(x,y,width,参数)对矩阵mn绘出m组每组n个宽度为width的条形图，以种类为x轴以实际数值为y轴，且参数有grouped和stacked之分。3.4其他统计图形绘制（area、polar、pie3、errorbar、compass、feather）

28、 其他统计图形常包括面积图area（X,Y）、三维饼图pie3(X,explode)、误差棒图errorbar(Y,E)、三维火柴棍图stem3(,filled)、罗盘图compass（X,Y）、羽毛图feather（U,V）、极坐标下柱状图rose()、极坐标下曲线polar（）等等也都是比较容易理解的，相关参数和具体应用读者可自行参阅相关资料弄懂，应结合数据的具体要求进行绘制，以下为部分案例：图8.MATLAB其他统计(tngj)图形绘制3.5特殊(tsh)图形(txng)绘制 数学的美体现在数学理论本身的奇特、微妙、简洁有力以及建立这些理论时人的创造性思维.MATLAB出色的数据可视化和

29、图像处理功能，几乎可以满足一般实际工程和科学计算中所有图形图像的需要。我们可以根据需要选择直角坐标、极坐标、柱坐标和球坐标等坐标系绘制平面曲线、空间曲线、空间曲面的表面图和网面图，还可以绘制直方图、向量图、柱状图等.此外，MATLAB还可以对图形进行标注、添色、变换视角等的加工和色彩控制、局部视图及动画等的操作以实实在在的表达各种理想的图形状态。3.5.1等值线图：在军事、地理等学科中经常会有遇到等值线，对于地质工作者来说更是屡见不鲜。MATLAB中有许多绘制等值线的命令，主要有如下：Contourc（）：该命令可用于计算等值线矩阵C，该举证可用于Contour、Contour3等Contou

30、r（）：二维等值线绘制，可以看作一个三维曲面在xy平面的投影Contour3（）：绘制三维等值线图，如图1Contourf（ ）：填充二维等值线，将相邻等值线用同一颜色填充，如图2Clabel（）：用于在二维等值线中添加高度标志，如图3% 以上命令为同一命令，括号内具体命令请自行参考资料图9.MATLAB等值线图绘制(huzh)3.5.2三维曲线图、曲面(qmin)图、网格图 MATLAB提供(tgng)了丰富的三维作图命令，如三维曲线plot3()、三维网格（mesh）、三维曲面（surf）、三维网格+等值线（meshc）、三维曲面+等值线（surfc）等等，每种类型又分为数据作图和函数作图

31、，不同类型语法格式也是丰富多彩的。在绘制过程中“数据网格化”举足轻重，下面为实习过程中绘制的一些案例： 图10.MATLAB三维曲线、网格图、曲面图绘制三维图形为立体图形，和二维图有了本质的不同。首先三维图有了视角view的概念，不同的视角对于同一物体会产生不同的视觉效果；颜色的变化color,光线的照射light角度也会对视觉效果产生影响。MATLAB提供了丰富的立体图形修饰命令，当然其作图的内容也比二维图形复杂得多，以下案例为部分修饰后的三维立体图形欣赏： 图11.MATLAB三维图形(txng)特殊效果展示通过(tnggu)MATLAB绘图化数学抽象(chuxing)为数学直观，透过直观

32、的数学图形，学习者可以发现数学思维的理性艺术的美，从而激发人们创造数学的激情和启发人们探求真理的思路，同时还可以利用数学的审美标准来检验科学真理。4. MATLAB统计模型与统计分析数理统计作为数学学科中应用最广泛的学科之一，它研究的对象是不确定的随机世界。而数学模型的功能主要包含两个方面：一是解释；二是预测。解释主要是分析相关关系或者因果关系；预测则是基于有限知识对未知领域的某种估计或者推断。就数学模型的功能而言，预测包括两种基本情况：一是样品内的预测，可以叫做内插或者估测；二是样品外的预测，可以称为外推或者推测。内插又可以分为两种情况：一是针对已知点；二是针对未知点。外推也分为两种情况：一

33、是针对过去，二是针对未来。最常用的就是针对已知点的内插预测和针对未来的外推预测。MATLAB提供了一个专门用于统计计算的统计工具箱，其中包括基本统计分析、多元统计分析、回归分析、寿命分析、时间序列分析等。统计工具箱涉及大量的统计知识和统计命令，可以在MATLAB的Help界面的左下栏进行搜索，在该栏中可以看到StatisticsTools，即统计模型工具箱。这部分可以说是MATLAB跟地质统计学乃至数学地质联系最紧密的部分，下面将逐一探讨。4.1方差分析（Analysis of variance）方差分析(Analysis of Variance，简称(jinchng)ANOVA) 分析试验（

34、或观测）数据(shj)的一种统计方法，又称“变异数分析”或“F检验”，是R.A.Fishe发明(fmng)的，用于两个及两个以上样本均数差别的显著性检验，分析各种因素及因素之间的交互作用对研究对象某些指标值的影响。由于各种因素的影响，研究所得的数据呈现波动状。造成波动的原因可分成两类，一是不可控的随机因素，另一是研究中施加的对结果形成影响的可控因素。其中单因素方差分析研究单个因素对观测变量的影响，多因素方差分析用来研究两个及两个以上控制变量是否对观测变量产生显著影响。假设下表为某一工程师根据多年实践总结的对不同金矿类型的预测的判得率，下面以单因素方差分析为例探索MATLAB方差分析的应用：表1

35、. 工程师对不同金矿类型的预测的判得率类型（oC）石英脉型破碎带蚀变岩型细脉浸染型砂金风化壳型得率（%）909796848492939683868892958584（%）平均得分9094958584总平均得分 89.6p，table，stats=anova1（X）%X：为数据矩阵，列为组间误差，行为组内误差%p：输出F统计量的上测概率%Table:输出方差分析表%Stats:输出各种有关统计量运行结果如下：图12.单因素(yn s)方差分析图表2.单因素(yn s)方差分析表（One-way ANOVA Table）Source方差来源SS平方和Df自由度MS 均匀FProbFCOLUMNS系

36、统误差435.84108.9556.846.95139*e-009ERRORS随机误差28.75151.917TOTAL464.5519从方差分析表的最后(zuhu)一项F统计量的上侧概率为6.95139*e-009，表明该统计量落入拒绝域。可以得出结论：原假设不成立，即不同类型的金矿床就有明显影响。事实上F统计量的上侧概率小于0.05表面水平有显著影响，而F统计量的上侧概率小于0.01表面水平有非常显著的影响。此外多因素方差分析显得稍复杂一些，公式为anovan().如双因素分析语句:p=anova2(X,reps)或p,table,stats=anova2()，方差分析图和方差分析表类似，

37、如其3个概率值p分别为因素1、因素2及双因素交互作用项起作用的结果。4.2回归分析（Regression analysis）回归模型是地学领域应用最广泛的统计模型之一，本身的内容十分丰富，回归模型根据实际的用途又可派生出不同的回归，如线性回归、岭回归、稳健回归、非线性回归等。利用MATLAB展开线性回归分析的两种途径： = 1 * GB3 运用矩阵运算规则 = 2 * GB3 利用回归分析程序包（其中一元线性回归还可以利用一次多项式拟合函数）。在对MATLAB比较熟悉的情况下两种情况都很方便。比较而言，调用回归分析程序包更为便捷，因为这种方法可以同时给出多种统计量。因此笔者回归分析乃至后面探讨

38、的各种统计分析均以函数调用的方式进行举例。4.2.1多元线性回归实例本例以资看二班的第三次实习“多元线性回归分析法进行矿床统计预测”为例，为方便计算我们已将矩阵数据的第二列全部设置为1，第一列为因变量向量Y，第2，3列分别为x1，x3. X=a(:，2 3 4)为样本矩阵，对应于回归系数的数据矩阵。当然也可以直接由不加常数列原始文件通过语句X=ones(32，1)，x2，x3直接获得。 a=data; X=a(:，2 3 4); y=a(:，1); b，bint，r，rint，stats=regress(Y，X)%输出b为回归系数的估计值，维数31，bint为b的置信水平维数32% r为残差向

39、量Y-Xb，维数321% rint为r的95%置信区间， 维数322%stats为回归模型的检验统计量，即决定系数R2，F统计量，概率值p.%X后缺省置信水平alpha，默认取0.05得系数：b =1.1904 0.3311 0.0482 统计量：stats = 0.3790 8.8508 0.0010 2.5237即得回归方程模型：y=1.190+0.331*x1+0.0482*x2，实际表明此方程系数与二班实习的手工计算版本是吻合的。统计检验：线性回归分析的主要统计检验类型常含：相关系数检验（含复相关系数R、简单相关系数、偏相关系数、部分相关系数）、标准误差检验（s,变异系数v=s/y-）

40、、F检验、T检验、DW检验、共线性诊断和分析。在具体过程中我们不必面面俱到，选取其中几个进行检验即可。拟合(n h)优度：R2=0.379指因变量Y(总评分(png fn)的37.9%。F检验(jinyn)： F值超过 F检验的临界值finv(0.95,2，29)=3.33, finv(0.99,2,29)均小于8.85，p远小于A=0.05或0.01，F统计量检验通过。利用F累计分布函数fcdf,可以将F统计量转换为概率值即F显著性Sig，公式为Fsig=1-fcdf(F统计量，回归自由度，剩余自由度)只要sig小于某个拟定的显著性水平，如0.05，F检验就可以通过。可见，sig的好处是不用

41、查表就可以判断，非常方便。DW检验：输入语句 rcoplot(r，rint) 作残余图如下图13回归分析残余图分析DW检验的本质是模型残差序列的自相关分析。残差序列的分析有多种角度。通过残差序列的变化区间图示可以看出异常样品。对于本例，从残差样品顺序图中可以看出，第6和第23个单元（单元发生重排，不同于原始给定单元）的数值变化范围超过了0平均线，为异常值。表明这些单元可能受到其他因素的影响较强，可以进行剔除后再进行线性回归。此外还可以运用标准误差检验，t检验，偏相关系数分析来检验，这里就不一一详述。4.2.2 MATLAB逐步回归（Stepwise）变量较多时也可以使用使用MATLAB统计工具

42、箱的逐步回归命令，逐步回归就是一种从众多自变量中有效地选择重要变量的方法。基本思路是，先确定一个包含若干自变量的初始集合，然后每次从集合外的变量中引入一个对因变量影响最大的，再对集合中的变量进行检验，从变得不显著的变量中移出一个影响最小的，依此进行，直到不能引入和移出为止.引入和移出都以给定的显著性水平为标准.逐步回归语句(yj)为Stepwise(x，y，inmodel,Penter,Premove).x为自变量构成(guchng)的矩阵，y为因变量构成(guchng)的向量，inmodel为包含在初始模型中的逻辑向量或者指示变量，Penter为变量引入的临界概率值，默认为0.05；为变量剔

43、除的临界概率值，默认为0.1. 同样以上述实习数据为例如下：图14.MATLAB逐步回归分析图-1stepwise命令产生图形窗口包括 = 1 * GB3 回归系数窗口：显示系数point的大小及回归系数的误差error变化范围bar.右侧可相应显示相关系数coeff.t统计量t-stat.概率值p-val。 = 2 * GB3 统计参数窗口：RMSE:目前模型误差的均方根，模型拟合的标准误差。R-square:回归方程的决定系数。F:回归模型的统计量F;p:与分析相关的显著概率. = 3 * GB3 Model History窗口：一个小窗口和三个按钮，通过小窗口可观察RMSE是否在下降。通

44、过点击鼠标可以决定将变量移进或移除，将三个变量全部移进后如图所示：图15.MATLAB逐步回归分析(fnx)图-2Y = -1.52 + 0.344 x1 + 0.047 x2t-stat -0.47 1.03 2.52p-val 0.64 0.31 0.02可见(kjin)b1.b2系数(xsh)跟前面差不多，但b0差别较大，存在进一步探讨的可能。逐步回归更适合于较多变量，在模型已经确定的情况下与原回归计算的结果可能不一样。另外MATLAB还是很具人性化的，回归结果可通过Export输出方便进一步分析。4.3主成分分析 (Principle component analysis)在统计分析中

45、有另外一类问题，就是提取信息的精华，这就是因子分析和主成份分析。当要研究问题时，通常的做法是提取大量的变量（信息），选择一组个数较少，且互不相关，并带有原样本的大部分的信息的另一组变量。MATLAB编写程序对主成分分析分析比较复杂，下面仍以MATLAB统计工具箱为例进行分析。以某地土壤数据（见附件）为分析目标，采样编号即样品组数，X1-X4为土壤不同的物质（沙、淤泥、粘土、有机物）含量百分数，X5为土壤PH值。下面用主成份分析来提取一组主成份：a=load(C:UsersAdministratorDesktopMATLAB practicecffx.txt)； pc，score，latent，

46、tsquare=princomp(a)% Pc:相对于特征值latent的特征向量，维数（55）% Score：主成份Yi中的元素，维数（205）% Latent：存放从大到小的特征值，维数（51）% Tsquare:t平方统计量，维数（201）计算结果（Pc Score Latent Tsquare略）：latent = 223.1197 8.7388 1.4469 0.2196 0.00214 可以看到最大的特征值已经(y jing)占到全部特征值的95,.5% ，即只抽取一个主成份(chng fn)就可以代表全部Xi的信息了，这显然是不合理的，原因是尺度问题。现在(xinzi)将变量X1

47、去掉，用剩下的个变量来分析。 xx=a(:，2:5); pc，score，latent，tsquare=princomp(xx)计算结果：latent = 86.0017 7.5521 1.4276 0.2192 第一个特征值仍占90%，尺度问题仍没有解决。现在将数据归一化，即对每一列的数据初上它们的标准差，然后进行主成分分析。 stdxx=std(xx); xxx=xx./repmat(stdxx，20，1); pc，score，latent，tsquare=princomp(xxx)计算结果：latent = 1.8023 1.5726 0.4021 0.2230前三个主成份占总和的94%

48、 ，因此我们成功的抽取了三个主成份，他们是：Y1=0.710*x2+0.702*x3+0.026*x4+0.048*x5Y2=0.181*x2+0.245*x3+0.834*x4+0.460*x5Y3=-0.148*x2+0.702*x3+0.026*x4+0.048*x5Y1主要由x2、x3组成，而它们分别表示淤泥、粘土含量，可以认为其代表了土壤的物理组成成分；Y2主要由x4组成，可以认为为土壤肥力的某种成分；而Y3主要由x5组成代表酸性指标，可以认为是某种化学成分。当然进行主成份分析时应灵活应用，有时对数据标准化，有时则不用。利用函数pareto，可以得到特征值的百分比和累计百分比。百分比

49、以柱形图的形式出现。从柱形图可以可以各个主成份的方差贡献率，从折线可以看出提取几个因子后方差累积贡献,如下图所示.此外，借助gname函数，可以得到载荷图上的变量标志和主成份得发图上的样品标签。还可以主成份的残差分析: Resi,Reco=pcares(data,Ndim)，Sumsq=sumsqr(Resi)，对Sumsq和Ndim用plot作图，可以得到主成份残差平方和随主成份提取变量的变化图，这里不深入探讨。图16. 特征值的百分比和累计(li j)百分比作图 4.4聚类分析(Cluster analysis)聚类分析（无监督(jind)或称无导师的Unsupervised）是根据研究对

50、象之间各种特征标志的相似程度或相关程度的大小，可将它们进行分类(fn li)归组。地学研究中的分类问题较多，如岩石分类、矿物分类、构造期次研究、古气候古环境划分等，这些都有可能需要利用聚类分析来研究。利用MATLAB开展聚类分析，关键在于学会调用子程序以及读懂聚类分析的过程和结果。为方便计算，我们仍以前一例中土壤成分的主成分分析数据为例进行层次聚类： clear a=load(C:UsersAdministratorDesktopMATLAB practicejulei.txt) %数据导入 x=a; y=pdist(x，Euclid);%计算样品之间的欧式距离，此语句method含多种距离算

51、法 m=squareform(y);%排出距离矩阵 z=linkage(y，single); %根据最短距离法层次聚类，此语句method含多种分类方法 h，t=dendrogram(z，0);%画出聚类谱系图xlabel(聚类层次);ylabel(土壤样品编号);title(土壤样品聚类分析谱系图); t=cluster(z，3);%指定(zhdng)分类数目和结果 find(t=2)ans = 3 5 10 11 12 13 14 15 16 17 20 图17.MATLAB聚类分析谱系(px)图 对于上面(shng min)的例子，关键是要读懂聚类结果t向量 3 3 2 1 2 3 3

52、1 1 2 2 2 2 2 2 2 2 3 3 2，表明1类样品3个，2类11个，3类6个.此外为满足图形表示的特殊需要，需要采用特殊的语法格式：如图形的局部显示、图形颜色的局部改变、图形的方向调整等等，我们可以根据个人喜爱和偏好进行更改和修饰。4.5判别分析(Discriminant function analysis)判别分析为广义的聚类分析，在一批样品已经明确分类的情况下，主要用于判定样品所属类别。判别分析有多种方法，包括距离判别法、Fisher判别法、Bayes判别法、逐步判别法等。其中最基本的是距离判别法，该判别法是基于马氏距离发展起来的一种线性判别技术。MATLAB提供三种判别方式

53、：线性判别、二次判别、和马氏判别。简单的判别分析命令为:Class=classify(Sample,Training,Group)Sample为待判样品矩阵；Training为训练样品矩阵；Group为训练样品的类别向量。注意：Training和Sample 的列必须相同，Training 与Group 应有相同的行数。更复杂一点的语句为： Class,Error=classify(Sample,Training,Group,Type，Prior)Err为误判误差率=训练样品中归类错误的样品数先验概率/训练样品中的样品数。Type有三种选项，即linearquadraticmahalanobi

54、s. Prior为数值向量，每一个数值对应于group中相应类别的发生概率。我们在矿床统计预测实习3判别分析时，是通过先计算样本均值，计算样本协方差，构造判别函数并计算得分，判别函数检验，待判样品归类的程序进行的。下面探讨MATLAB直接进行MATLAB判别分析，仍以咱们班实习三判别分析原始数据为例：实习数据分为两大类有矿和无矿，为方便计算将其转变成如附件所示的格式成172维矩阵，前7行为含矿单元x1、x2值，后10个单元为无矿单元x1,x2值：导入训练(xnlin)数据a(维数172)和欲判数据s（维数52）后语句如下： clear clf t=a; g=1 1 1 1 1 1 1 0 0

55、0 0 0 0 0 0 0 0； Class,Error=classify(s,t,g,linear,0.41,0.58)；Class = 0 0 1 1 0Error = 0.1251由MATLAB的计算结果误差率为：12.51%，跟手工(shugng)版的11.76%相差无几(xing ch w j)（原因是本例中数值向量0.41,0.58）；欲判数组结果为 00110，第第3第4个预测含矿，第1、2、5个预测无矿，与原手工版的结果完全一样。MATLAB过程方便快捷，可以不用建立模型直接预测，乃统计模型的集大成者。4.6自相关分析 ( Autocorrelation analysis)自相

56、关分析可以看作是相关分析的一种推广，为一种有效的时间、空间序列分析方法。不言而喻自相关系数是相关系数的推广，自相关系数的计算公式和检验统计量都十分简单，利用MATLAB自相关分析函数和偏自相关分析函数可以非常方便的计算函数值并画出相关函数图。有了计算结果和图像，就可以根据有关知识进行时间和空间序列的自相关分析。MATLAB的自相关函数为autocorr，语法结构为ACF,Lags,Bounds=autocorr(Series,nLags,M,nSTDs); （）内Series代表单变量时间或空间序列样本路径生成的向量，适用于系统分析的某种信号，或者某种定期、定点采样的时、空观测值。其后的几个为

57、可选项（nLags为正整数标量，M为非负整数标量，nSTDs代表标准离差范围的正标量）。 笔者以2011版高等教育出版社MATLAB地理数据分析中第八章给定的某海域测得的年平均海平面高度数据（100年）为例，导入数据后将数据矩阵（第一列为时间序列，第二列为海平面高度）更名为H（可直接重命名），输入如下语句代码： figure(1) plot(H,r.); hold on plot(H,b- hold off xlabel(Time); ylabel(Height); figure(2) ACF,Lags,Bounds=autocorr(H); %调用(dioyng)自相关函数 C=Lags;A

58、CF; %将Lags与ACF对齐(du q)并转置 autocorr(H); %绘制(huzh)自相关函数图 hold off C,Bounds 结果所得ACF序列为1，0.8054，0.4353，0.0323，-0.2597.相应的Lags序列的数值为0，1，2，3，4。其中0滞后相应的ACF值ACF0没有任何消息。Bounds给出了基于95%置信区度的正负二倍标准差0.2和-0.2。时间序列变化曲线如下图，可以看到ACF图为Lags为横轴，ACF为纵轴，以离散点线表示出来。图18. 某海域海平面平均高度变化曲线图（100年）图19. 某海域(hiy)海平面高度ACF图（nLags=20）既

59、然海平面高度变化不是不是随机的，而有可能具有(jyu)周期性，因为原始序列的散点图是波动变化的。通过典型时间序列的ACF图可知，海平面百年变化信号(xnho)不是白噪音，也没有趋势性。且可以看出 = 1 * GB3 自相关系数在坐标图上呈波动衰减趋势，每隔11年出现一个峰值和一个谷值，这些峰值和谷值周期性的突破二倍标准误差带。 = 2 * GB3 相当多的ACF值在（+/-）1.96/10=0.196间；第一个谷值对应的时滞为5，第二个谷值对应的时滞为16.两者之差为11；第一个峰值对应的时滞为10，第二个谷值对应的时滞为21.两者之差也为11。由此初步推断：海平面高度变化的时间序列存在一个1

60、1年左右的波动周期，自相关系数图的波动衰减特征暗示着时间序列具有双重周期。借助ACF不仅能反映过去对未来的直接影响，还能反映过去对未来的间接影响。当然此次研究性学习笔者还只是对自相关分析的ACF进行了分析，对ACF的检验和偏自相关系数等系列同等重要的内容还没有涉及，更多的分析和研究还需进一步的工作。5.结语和感想(Conclusion)市场上关于MATLAB的书可谓汗牛充栋，网页上关于其应用的资料数不胜数，但有许多也只是简单的翻译MATLAB帮助文件中的结果，基于MATLAB完整讲述分析案例的图书非常罕见，基于MATLAB全面、系统地分析地理数据的图书更是鲜见。因此本次探索性学习集思广益不断挖