《23 MATLAB数理统计.doc》由会员分享，可在线阅读，更多相关《23 MATLAB数理统计.doc（19页珍藏版）》请在三一办公上搜索。
1、23 MATLAB数理统计本章要点： MATLAB中多元回归分析的实现方法 MATLAB中逐步回归分析的实现方法 MATLAB中聚类分析的实现方法 MATLAB中主因子分析的实现方法 MATLAB中对应分析的实现方法 MATLAB中时间序列分析的实现方法23.1 多元回归分析% 回归分析的MATLAB代码% 0. 清除工作区间所有变量clear all;% 1. 线性回归X1=2:11;% 输入变量X1X=ones(size(X1);X1;% 输入变量X(在X1前加一列1并转置)y=6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59; % 输入变量y并
2、转置b1=regress(y,X);% 回归结果: 系数b1b1,b2=regress(y,X);% 回归结果: b1,95%置信区间b2b1,b2,b3=regress(y,X);% 回归结果: b1,b2,残差b3b1,b2,b3,b4=regress(y,X);% 回归结果: b1,b2,b3,残差95%置信区间b4b1,b2,b3,b4,b5=regress(y,X);% 回归结果: b1,b2,b3,b4,b5(R2,F,p)alp=0.85; b1=regress(y,X,alp);% 回归结果为 b1(左端同上),输入置信度alp% 2. 非线性回归clear all; % 清除
3、工作区间所有变量x=2:11; % 输入变量xy=6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59; % 输入变量ybeta0=0.1 0.1; % 设置参数初始值beta0beta,R,J=nlinfit(x,y,f1,beta0); % 回归结果: 系数beta,残差R,Jacobi矩阵Jbeta=nlinfit(x,y,f1,beta0); % 回归结果: 系数beta% 2-1. 回归函数f1的代码(另存为工作目录下)% function y=f1(beta,x)% a=beta(1);% b=beta(2);% y=x./(a*x+b);2
4、3.2逐步回归分析23.2.1 理论基础逐步回归分析是从一个自变量开始，视自变量Y作用的显著程度，从大到地依次逐个引入回归方程。 当引入的自变量由于后面变量的引入而变得不显著时，要将其剔除掉；引入一个自变量或从回归方程中剔除一个自变量，为逐步回归的一步；对于每一步都要进行Y值检验，以确保每次引入新的显著性变量前回归方程中只包含对Y作用显著的变量；这个过程反复进行，直至既无不显著的变量从回归方程中剔除，又无显著变量可引入回归方程时为止。 MATLAB中逐步回归的命令是stepwise。该命令提供了一个交互式画面，通过此工具可以自由地选择变量进行统计分析。调用格式为：stepwise(x，y，in
5、model，alpha)其中x是自变量数据，y是因变量数据，分别为nm和nl矩阵，inmodel是矩阵的列数指标(缺省时为全部自变量)，alpha,为显著性水平(缺省时为0.5)。结果产生三个图形窗口，在stepwise plot窗口，虚线表示该变量的拟合系数与0无显著差异，实线表示有显著差异，红色线表示从模型中移去的变量；绿色线表明存在模型中的变量，点击一条会改变其状态。在stepwise Table窗口中列出一个统计表，包括回归系数及其置信区间，以及模型的统计量：剩余标准差(RMSE)、相关系数 (R-square)、F值和P值。23.2.2 算例23.2.2.1 问题描述试验原始数据取自
6、于某工程16处的样品(i=26)，描述样品特性的自变量指标为7个、因变量指标为1个，即分别为x(i,1), x(i,2), , x(i,7)和y(1), y(2), , y(7)，试验结果见表23-1。 试确定x和y之间的逐步回归关系。表23-1 逐步回归分析原始数据表样品编号ix(i,1)x(i,2)x(i,3)x(i,4)x(i,5)x(i,6)x(i,7)y(i)10.704.060.0150.0195.00.9550.9101.1821.105.064.0152.5217.50.9500.8801.3631.205.582.5187.5242.50.9150.8551.5841.204
7、.557.5130.0177.50.9700.9251.2050.855.062.5150.0217.50.9500.8801.3761.005.092.5160.0242.50.9450.8551.5871.155.597.5175.0250.00.9300.8351.6281.253.555.0145.0177.50.9550.9251.0890.605.065.0152.5217.50.9450.8801.22100.905.077.5160.0192.50.9400.8551.60111.305.595.0175.0250.00.9300.8351.65120.653.570.0135
8、.0182.50.9700.9201.10131.104.574.0150.0217.50.9500.8751.18141.155.577.5160.0255.00.9450.8351.32151.255.592.5165.0257.50.9400.8251.56161.305.0105.0180.0262.50.9250.8151.5423.2.2.2 程序代码% 1 输入原始数据clear all;x1=0.7,1.1,1.2,1.2,0.85,1.0,1.15,1.25,0.6,0.9,1.3,0.65,1.1,1.15,1.25,1.3; x2=4.0,5.0,5.5,4.5,5.0,
9、5.0,5.5,3.5,5.0,5.0,5.5,3.5,4.5,5.5,5.5,5.0; x3=60.0,64.0,82.5,57.5,62.5,92.5,97.5,55.0,65.0,77.5,95.0,70.0,74.0,77.5,92.5,105.0;x4=150,152.5,187.5,130,150,160,175,145,152.5,160,175,135,150,160,165,180;x5=195,217.5,242.5,177.5,217.5,242.5,250,177.5,217.5,192.5,250,182.5,217.5,255,257.5,262.5;x6=0.95
10、5,0.950,0.915,0.970,0.950,0.945,0.930,0.955,0.945,0.940,0.930,0.970,0.950,0.945,0.940,0.925;x7=0.910,0.880,0.855,0.925,0.880,0.855,0.835,0.925,0.880,0.855,0.835,0.920,0.875,0.835,0.825,0.815;Y=1.18,1.36,1.58,1.20,1.37,1.58,1.62,1.08,1.22,1.60,1.65,1.10,1.18,1.32,1.56,1.54; % 2 回归分析X=x1,x2,x3,x4,x5,x
11、6,x7;stepwise(X,Y,0.1); = x6,x7可去掉X=x1,x2,x3,x4,x5;stepwise(X,Y,0.1); = x1,x7可去掉X=x2,x3,x4,x5;stepwise(X,Y,0.1); = 得到回归系数23.2.2.1 结果分析采用逐步回归分析时不同步骤的统计判断情况见表23-2。表23-2中，RMSE为多元回归分析均方差平方根，其值越小、则拟合效果越好；R2为复相关系数的平方，其值越大拟合效果越好。RMSE和R2分别根据公式(23-1)和(23-2)进行计算。 (23-1) (23-2)式中，m为各步回归分析的变量数；n为样本数，n=16；为第i个样本
12、的估计值；yi为第i个样本的实测值；xik为第i个样本第k个参数的实测值；i为待估参数。由表23-2可以看出，第六、第七参数在第一步即可剔除；剩余五个参数重新计算后，第一参数在第二步可以剔除；此后，再不能剔除各个参数了。因此，第二、第三、第四和第五参数被保留了下来。由表23-5可知， Y与x2、x3、x4和x5之间的关系为 (23-3)由式(23-3)可以看出， x2和x3为主要参数，x4、x5为次要参数；x2越大则Y越大。表 23-2 不同步骤的判断结果步骤回归系数置信区间RMSER2列数值下限值上限值010.004979-0.17340.18330.075610.927120.25030.
13、12610.374530.0081350.0016830.0145940.01635-0.0048840.037585-0.005943-0.009655-0.00223611.79-8.54532.137-0.1135-4.163.933110.003521-0.16290.170.072530.916120.22430.13220.316430.0092940.0049840.013640.0042530.000094280.0084115-0.005156-0.008034-0.002277220.22430.13730.31330.069160.916130.0093030.00525
14、20.0133540.0042620.00035880.0081665-0.005152-0.007867-0.00243723.3聚类分析23.3.1 理论基础聚类分析是将数据分为若干类的一种多变量统计分析方法，通常有层次聚类和非层次聚类两种类型。层次聚类使用不同的样本距离和类间距离对样品进行聚类，而非层次聚类主要包括C-均值聚类和模糊C-均值聚类。聚类分析可用于分析样品的类型和特性、类别之间的相似性和差异性。假设共有n个样品，第i土个样品属于类别Y=Y1,Y2,YK之一（K为类别总数），每一样品使用p个指标x(i,1), x(i,2), , x(i,p)(i=1,2,n)来描述。根据试验可
15、得观测数据矩阵X(=x(i,j)。下面使用层次聚类法(简作HC)、模糊C-均值聚类(简作FCM)对样品类别进行聚类分析与评判。HC首先将n个样本中的每一样本看作一类，再把两个相距最近的样本合成为一个新的类别，然后计算新的(n-1)个样本或样本组之间的距离、将相距最近的样本或样本组合并为一个新类。重复这一过程，直到所有样本组合为一类。FCM是根据隶属度大小把样本归属于某类的数据聚类技术。聚类开始时，假定初始聚类中心并依据隶属度大小将任一样本赋予一类。此后，将聚类中心移动到使目标函数取最小值的合适位置。样本距离和类间距离可以有不同定义。样本距离通常有欧氏距离、城块距离、闵氏距离、余弦距离、相关距离
16、、海明距离、切比契夫距离等，类间距离计算方法通常有类平均法、全距离法、加权均值法、重心法、中间值法。23.3.2 算例23.3.2.1 问题描述下面以某处试验资料为例，说明使用聚类分析方法确定样品类型、不同聚类方法确定结果的比较。试验原始数据取自于某工程的26个样品(i=26)，描述样品特性的重要指标选为10个，即x(i,1), x(i,2), , x(i,10)，试验结果见表23-1。 由于分析的样品类别有三类(见表23-3)，聚类数反映了各类别的物理力学性质差异，因此，聚类分析时聚类数取为3。为了消除不同指标使用不同量纲的影响，使用式(23-4)对原始数据矩阵X进行无量刚化处理: (23-
17、4)式中，xjM和xjN分别为第j个指标的最大值和最小值。表23-3 试验成果简表样品编号实际类别x(i,1)x(i,2)x(i,3)x(i,4)x(i,5)x(i,6)x(i,7)x(i,8)x(i,9)x(i,10)112.7333.618.60.96135.014.51711.30.414.54212.7325.319.50.75438.316.92011.30.453.73312.7223.619.80.69833.913.81611.60.285.84412.7325.919.70.74533.613.61812.40.305.62512.7325.819.60.75237.216.
18、02012.10.295.80612.7332.618.60.94637.716.41511.00.444.24712.7231.518.80.90333.913.71511.90.364.99812.7336.618.61.00538.116.91912.40.523.70912.7325.919.50.76337.616.42011.60.315.431012.7325.919.60.75438.016.61812.40.325.311112.7327.719.00.83537.716.42011.60.443.881222.7335.218.60.98438.116.71511.30.4
19、54.191322.7335.018.60.98138.416.81711.00.444.361422.7335.118.60.98338.416.91513.50.483.871522.7336.018.11.05137.316.11214.00.395.041622.7325.519.60.74836.815.72011.30.354.791732.7436.118.41.02740.718.5129.90.613.121832.7438.918.11.10342.619.8139.10.682.921932.7436.518.01.07839.417.6108.50.782.482032
20、.7437.818.11.08639.817.9128.80.623.142132.7441.217.41.22340.918.61010.20.812.532232.7437.518.11.08139.817.9128.30.653.062332.7438.418.31.07240.718.5109.90.682.812432.7437.218.21.06639.617.81510.50.583.352532.7437.718.31.06240.718.51510.20.643.052632.7441.617.51.22542.920.11610.80.802.6123.3.2.2 程序代码
21、% 1 输入原始数据clear all;X0= 1, 2.73, 33.6, 18.6, 0.961, 35.0, 14.5, 17, 11.3, 0.41, 4.54;1, 2.73, 25.3, 19.5, 0.754, 38.3, 16.9, 20, 11.3, 0.45, 3.73;1, 2.72, 23.6, 19.8, 0.698, 33.9, 13.8, 16, 11.6, 0.28, 5.84;1, 2.73, 25.9, 19.7, 0.745, 33.6, 13.6, 18, 12.4, 0.30, 5.62;1, 2.73, 25.8, 19.6, 0.752, 37.2
22、, 16.0, 20, 12.1, 0.29, 5.80;1, 2.73, 32.6, 18.6, 0.946, 37.7, 16.4, 15, 11.0, 0.44, 4.24;1, 2.72, 31.5, 18.8, 0.903, 33.9, 13.7, 15, 11.9, 0.36, 4.99;1, 2.73, 36.6, 18.6, 1.005, 38.1, 16.9, 19, 12.4, 0.52, 3.70;1, 2.73, 25.9, 19.5, 0.763, 37.6, 16.4, 20, 11.6, 0.31, 5.43;1, 2.73, 25.9, 19.6, 0.754,
23、 38.0, 16.6, 18, 12.4, 0.32, 5.31;1, 2.73, 27.7, 19.0, 0.835, 37.7, 16.4, 20, 11.6, 0.44, 3.88;2, 2.73, 35.2, 18.6, 0.984, 38.1, 16.7, 15, 11.3, 0.45, 4.19;2, 2.73, 35.0, 18.6, 0.981, 38.4, 16.8, 17, 11.0, 0.44, 4.36;2, 2.73, 35.1, 18.6, 0.983, 38.4, 16.9, 15, 13.5, 0.48, 3.87;2, 2.73, 36.0, 18.1, 1
24、.051, 37.3, 16.1, 12, 14.0, 0.39, 5.04;2, 2.73, 25.5, 19.6, 0.748, 36.8, 15.7, 20, 11.3, 0.35, 4.79;3, 2.74, 36.1, 18.4, 1.027, 40.7, 18.5, 12, 9.9, 0.61, 3.12;3, 2.74, 38.9, 18.1, 1.103, 42.6, 19.8, 13, 9.1, 0.68, 2.92;3, 2.74, 36.5, 18.0, 1.078, 39.4, 17.6, 10, 8.5, 0.78, 2.48;3, 2.74, 37.8, 18.1,
25、 1.086, 39.8, 17.9, 12, 8.8, 0.62, 3.14;3, 2.74, 41.2, 17.4, 1.223, 40.9, 18.6, 10, 10.2, 0.81, 2.53;3, 2.74, 37.5, 18.1, 1.081, 39.8, 17.9, 12, 8.3, 0.65, 3.06;3, 2.74, 38.4, 18.3, 1.072, 40.7, 18.5, 10, 9.9, 0.68, 2.81;3, 2.74, 37.2, 18.2, 1.066, 39.6, 17.8, 15, 10.5, 0.58, 3.35;3, 2.74, 37.7, 18.
26、3, 1.062, 40.7, 18.5, 15, 10.2, 0.64, 3.05;3, 2.74, 41.6, 17.5, 1.225, 42.9, 20.1, 16, 10.8, 0.80, 2.61;for i=1:10, X(:,i)=X0(:,i+1);end;% 提取分类参数NM=size(X0);for i=1:NM(1), Y0(i)=X0(i,1);end;% 样品实际类型% 2 使用模糊C均值聚类方法进行聚类center,U,obj_fcn=fcm(X,3);% 模糊C均值聚类Y=sort(U);Y,I=sort(U);% 排序format short g% 实数输出格式
27、center; % 显示模糊聚类中心max(U); % 样品隶属度I(3,:) % 类别代码% 制作聚类图形maxU=max(U);% 样品隶属度index1=find(U(1,:)=maxU); % 找到属于第1类的点index2=find(U(2,:)=maxU); % 找到属于第2类的点index3=find(U(3,:)=maxU); % 找到属于第3类的点% 作图line(X(index1,1),X(index1,2),linestyle,none,marker,o,color,g);line(X(index2,1),X(index2,2),linestyle,none,marker
28、,x,color,r);line(X(index3,1),X(index3,2),linestyle,none,marker,+,color,k);% 显示每一类中心hold on;plot(center(1,1),center(1,2),kpentagram,markersize,7,LineWidth,2);plot(center(2,2),center(2,2),ksquare,markersize,7,LineWidth,2);plot(center(3,2),center(3,2),kv,markersize,7,LineWidth,2);% 3 使用层次聚类方法进行聚类Z = li
29、nkage(X,single,euclidean); % 样本数据,距离定义,聚合方法% single可为下述之一:% single/complete/average/weighted/% centroid/median/ward % 此时，第3项只能是euclidean% euclidean可为下述之一:% euclidean/seuclidean/cityblock/mahalanobis/minkowski/cosine/% correlation/spearman/hamming/jaccard/chebychevdendrogram(Z);% 画出谱图T = cluster(Z,3)
30、;% 将X聚合为3类xlabel(Observations,FontSize,20); hold on;ylabel(Distances of Observations,FontSize,20); hold on;set(gca,FontSize,16,LineWidth,3); hold onT=T;% 分类参数取“全部参数”、距离为欧氏距离下不同层次聚类方法的聚类结果Z = linkage(X,complete,euclidean); % 样本数据,距离定义,聚合方法T21= cluster(Z,3);T21=T21% 将X聚合为3类Z = linkage(X,average,euclid
31、ean); % 样本数据,距离定义,聚合方法T22= cluster(Z,3);T22=T22% 将X聚合为3类Z = linkage(X,weighted,euclidean); % 样本数据,距离定义,聚合方法T23= cluster(Z,3);T23=T23% 将X聚合为3类Z = linkage(X,centroid,euclidean); % 样本数据,距离定义,聚合方法T24= cluster(Z,3);T24=T24% 将X聚合为3类Z = linkage(X,median,euclidean); % 样本数据,距离定义,聚合方法T25= cluster(Z,3);T25=T25
32、% 将X聚合为3类Z = linkage(X,median,euclidean); % 样本数据,距离定义,聚合方法T26= cluster(Z,3);T26=T26% 将X聚合为3类% 分类参数取“全部参数”、聚类方法取类平均法、不同距离定义下的层次聚类结果Z = linkage(X,average,seuclidean); % 样本数据,距离定义,聚合方法T31= cluster(Z,3);T31=T31% 将X聚合为3类Z = linkage(X,average,cityblock); % 样本数据,距离定义,聚合方法T32= cluster(Z,3);T32=T32% 将X聚合为3类Z
33、 = linkage(X,average,seuclidean); % 样本数据,距离定义,聚合方法T33= cluster(Z,3);T33=T33% 将X聚合为3类Z = linkage(X,average,mahalanobis); % 样本数据,距离定义,聚合方法T34= cluster(Z,3);T34=T34% 将X聚合为3类Z = linkage(X,average,minkowski); % 样本数据,距离定义,聚合方法T35= cluster(Z,3);T35=T35% 将X聚合为3类Z = linkage(X,average,cosine); % 样本数据,距离定义,聚合方
34、法T36= cluster(Z,3);T36=T36% 将X聚合为3类Z = linkage(X,average,correlation); % 样本数据,距离定义,聚合方法T37= cluster(Z,3);T37=T37% 将X聚合为3类Z = linkage(X,average,spearman); % 样本数据,距离定义,聚合方法T38= cluster(Z,3);T38=T38% 将X聚合为3类Z = linkage(X,average,hamming); % 样本数据,距离定义,聚合方法T39= cluster(Z,3);T39=T39% 将X聚合为3类Z = linkage(X,
35、average,jaccard); % 样本数据,距离定义,聚合方法T310= cluster(Z,3); T310=T310% 将X聚合为3类Z = linkage(X,average,chebychev); % 样本数据,距离定义,聚合方法T311= cluster(Z,3); T311=T311% 将X聚合为3类25.3.1.3 结果分析图23-1是样本距离和类间距离分别使用欧氏距离法和类平均法得到的26个样本聚类直方图。表23-4是使用不同变量组合、不同样本距离、不同类间距离的层次聚类结果。表23-5为FCM的聚类结果。这两张表中，聚类条件表示选用不同变量的情况：情况1选用x(i,4)
36、和x(i,6)；情况2选用x(i,4)、x(i,5)和x(i,6)；情况3选用x(i,2)、x(i,4)、x(i,5)和x(i,6)；情况4选用x(i,1)、x(i,2)、x(i,4)、x(i,5)和x(i,6)；情况5选用全部10个变量；类别序列为类别编号(见表23-2)的排列；土类序列中的黑斜体表示误判土类。从表23-4可以看出，层次聚类时如果只选用一部分变量，误判率高达42.3%和53.8%；如果选用全部变量、样本距离不选用海明距离、类间距离采用类平均法，误判率可降至19.2%。由表23-5可以看出，如果FCM中所选变量很少（比如选用x(i,4)和x(i,6)或三个变量x(i,4)、x(
37、i,5)和x(i,6)），误判率为30%左右，这一误判率过大。然而，如果此时选择三个以上的变量，误判率则为19.2%左右。从实用观点出发，上述聚类分析方法可用于实际工程中类别的初步确定，对亚类划分也有一定的参考价值。如果聚类分析时使用合适的样本距离、类间距离并且较多选用指标，会获得比较好的聚类结果，可为样品类别鉴定和性质评判提供许多有用信息。 样品编号图23-1 样品层次聚类直方图表23-4 样品的层次聚类结果聚类条件样本距离类间距离类别序列误判率实际情况1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3情况1欧氏距离类平均法1 3 1 1
38、3 3 1 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 214/26=53.8%情况21 2 1 1 3 3 1 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 214/26=53.8%情况33 1 2 2 1 3 3 3 1 1 1 3 3 3 3 1 3 3 3 3 3 3 3 3 3 311/26=42.3%情况43 1 2 2 1 3 3 3 1 1 1 3 3 3 3 1 3 3 3 3 3 3 3 3 3 311/26=42.3%情况5全距离法2 1 1 1 1 2 2 3 1 1 1 2 2 2 2 1 3 3 3 3 3 3
39、3 3 3 35/26=19.2%加权均值法2 1 1 1 1 2 2 2 1 1 1 2 2 2 2 1 3 3 3 3 3 3 3 3 3 35/26=19.2%重心法2 1 1 1 1 2 2 2 1 1 1 2 2 2 2 1 3 3 3 3 3 3 3 3 3 35/26=19.2%中间值法2 1 1 1 1 2 2 2 1 1 1 2 2 2 2 1 3 3 3 3 3 3 3 3 3 35/26=19.2%城块距离类平均法2 1 1 1 1 2 2 2 1 1 1 2 2 2 2 1 3 3 3 3 3 3 3 3 3 35/26=19.2%闵氏距离2 1 1 1 1 2 2 2
40、 1 1 1 2 2 2 2 1 3 3 3 3 3 3 3 3 3 35/26=19.2%余弦距离2 1 1 1 1 2 2 2 1 1 1 2 2 2 2 1 3 3 3 3 3 3 3 3 3 35/26=19.2%相关距离2 1 1 1 1 2 2 2 1 1 1 2 2 2 2 1 3 3 3 3 3 3 3 3 3 35/26=19.2%海明距离1 1 2 1 1 1 2 1 1 1 1 1 1 1 3 1 3 3 3 3 3 3 3 3 3 37/26=26.9%切比契夫距离2 1 1 1 1 2 2 2 1 1 1 2 2 2 2 13 3 3 3 3 3 3 3 3 35/26=19.2%表 23-5 样品模糊C-均值聚类结果聚类条件类别序列总数正确错误误判率%实际情况1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3情况11 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 32618830.8情况21 2 1 1 2 2 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 32619726.9情况32 1
