当前位置:首页 >> 初中教育 >>

MATLAB技术论坛电子期刊第三期


MATLAB 技术论坛
电子期刊

编辑:xiezhh 制作:MATLAB 技术论坛 版权:MatlabSky?版权所有 网址:http://www.matlabsky.com

第3期

2010.11 No.3

中国权威 MATLAB 论坛核心期刊

MATLAB

技术论坛简介
Matlab 技术论坛 | Simulink 仿真论坛 | ——打造优秀、专业和权威的 Matlab 技术交流平台!

http://www.matlabsky.com /cn/net/org
论坛拥有 40 多个专业版块,内容涉及资料下载、视频教学、数学建模、科学计算、程序设计、GUI 开发、simulink 仿真、统计概率、拟合优化、扩展编程、算法研究、控制系统、信号通信、图像处理、经 济金融、生物化学、航空航天、人工智能、汽车设计、机械自动化、毕业设计等几十个方面!

论坛特色:
1.拥有强大的技术人员、热情严谨的管理团队 2.丰富多彩的Matlab电子资源免费共享 3.免费提供技术交流和在线解答 4.首个推出MATLAB函数百科中文帮助系统 5.国内唯一提供Matlab汉化包的团队

联系方式:
客服 QQ:1341692017 服务邮箱:matlabsky@gmail.com 支付宝账号:yuthreestone@163.com 郁磊 互动QQ群:http://www.matlabsky.com/thread-3-1-1.html

开通业务:
技术论坛:http://www.matlabsky.com 函数百科:http://wiki.matlabsky.com 官方博客:http://blog.matlabsky.com 读书频道:http://book.matlabsky.com 视频教程:http://video.matlabsky.com 有偿编程:http://paid.matlabsky.com





1 一元线性回归分析 ........................................................................................................................1 1.1 数据的散点图........................................................................................................................................ 2 1.2 调用 REGRESS 函数作一元线性回归分析 ............................................................................................ 3
1.2.1.regress 函数的用法................................................................................................................................................. 3 1.2.2.调用 regress 函数作一元线性回归分析 ................................................................................................................ 4 1.2.3.残差分析 ................................................................................................................................................................ 6

1.3 调用 REGSTATS 函数作一元线性回归分析........................................................................................... 7
1.3.1.regstats 函数的用法................................................................................................................................................ 7 1.3.2.调用 regstats 函数作一元线性回归分析 ............................................................................................................... 8

1.4 调用 ROBUSTFIT 函数作稳健回归....................................................................................................... 11
1.4.1.robustfit 函数的用法 ............................................................................................................................................ 11 1.4.2.调用 robustfit 函数作稳健回归............................................................................................................................ 13

2 一元非线性回归分析 ...................................................................................................................14 2.1 数据的散点图...................................................................................................................................... 15 2.2 调用 NLINFIT 函数作一元非线性回归分析........................................................................................ 16
2.2.1.nlinfit 函数的用法 ................................................................................................................................................ 16 2.2.2.选择合适的理论回归方程................................................................................................................................... 17 2.2.3.调用 nlinfit 函数作一元非线性回归分析............................................................................................................ 18 2.2.4.绘制一元非线性回归曲线................................................................................................................................... 18 2.2.5.参数估计值的置信区间....................................................................................................................................... 19 2.2.6.头围平均值的置信区间和观测值的预测区间 ................................................................................................... 19

2.3 利用曲线拟合工具 CFTOOL 作一元非线性拟合................................................................................ 21
2.3.1.cftool 函数的调用方式......................................................................................................................................... 22 2.3.2.导入数据 .............................................................................................................................................................. 22 2.3.3.数据的平滑处理................................................................................................................................................... 23 2.3.4.数据筛选 .............................................................................................................................................................. 23 2.3.5.数据拟合 .............................................................................................................................................................. 24 2.3.6.绘图控制 .............................................................................................................................................................. 26 2.3.7.后处理 .................................................................................................................................................................. 27

3 多重回归分析..............................................................................................................................27 3.1 调用自编 REGLM 函数作多重回归分析 ............................................................................................. 28
3.1.1.调用 reglm 函数作多重线性回归分析 ................................................................................................................ 28 3.1.2.调用 reglm 函数作二次回归分析........................................................................................................................ 30 3.1.3.拟合效果图........................................................................................................................................................... 31

3.2 调用 STEPWISE 函数作逐步回归......................................................................................................... 31
3.2.1.stepwise 函数的用法 ............................................................................................................................................ 31 3.2.2.调用 stepwise 函数作逐步回归 ........................................................................................................................... 32

MATLAB 技术论坛电子期刊 第 3 期

·1·

MATLAB 统计概率版块
资料来源:《MATLAB 统计分析与应用:40 个案例分析》,谢中华编著,北 京航空航天大学出版社,2010 年 6 月。

数据拟合——回归分析
在很多实际问题中,通常会涉及很多变量,需要研究变量之间的关系,很多时候变量之间的关系是 不确定的,需要用一个函数来近似表示这种关系。数据拟合就是根据变量的观测数据研究某些变量之间 的近似函数关系,用来帮助我们认识事物的内在规律和本质属性。本章将结合具体案例介绍用回归分析 方法进行数据拟合。

1 一元线性回归分析
现有全国 31 个主要城市 2007 年的气候情况观测数据,如表 1-1 所示。数据来源:中华人民共和国国 家统计局网站 2007 年环境统计数据。
表 1-1 全国 31 个主要城市 2007 年的气候情况观测数据 城 北 天 太 沈 长 上 南 杭 合 福 南 济 郑 武 长 广 南 海 重 成 贵 市 京 津 原 阳 春 海 京 州 肥 州 昌 南 州 汉 沙 州 宁 口 庆 都 阳 年平均气温 单位:℃ 14.0 13.6 14.9 11.4 9.0 9.0 7.7 6.6 18.5 17.4 18.4 17.4 21.0 19.2 15.0 16.0 18.6 18.8 23.2 21.7 24.1 19.0 16.8 14.9 年极端最高气温 单位:℃ 37.3 38.5 39.7 35.8 35.6 33.9 35.8 35.8 39.6 38.2 39.5 37.2 39.8 38.5 38.5 39.7 37.2 38.8 37.4 37.7 37.9 37.9 34.9 31.0 年极端最低气温 单位:℃ -11.7 -10.6 -7.4 -13.2 -17.6 -23.1 -21.7 -22.6 -1.1 -4.5 -1.9 -3.5 3.6 0.5 -7.9 -5.0 -1.5 -0.5 5.7 0.7 10.7 3.0 -1.6 -1.7 年均相对湿度 单位:% 54 61 59 55 47 68 58 58 73 70 71 79 68 68 61 60 67 70 71 76 80 81 77 75 全年日照时数 单位:小时 2351.1 2165.4 2167.7 2174.6 2647.8 2360.9 2533.6 2359.2 1522.2 1680.3 1472.9 1814.6 1543.8 2102.0 1819.8 1747.2 1934.2 1742.2 1616.0 1614.0 1669.1 856.2 935.6 1014.8 全年降水量 单位:毫米 483.9 389.7 430.4 535.4 261.2 672.3 534.2 444.1 1254.5 1070.9 1378.5 929.7 1109.6 1118.5 797.1 596.4 1023.2 9364.0 1370.3 1008.1 1419.3 1439.2 624.5 884.9 序 号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

石家庄 呼和浩特

哈尔滨

2010 年 11 月
昆 拉 西 兰 西 银 明 萨 安 州 宁 川 15.6 9.8 15.6 11.1 6.1 10.4 8.5 30.0 29.0 39.8 34.3 30.7 35.0 37.6

Matlab 数据拟合——回归分析
0.7 -9.8 -5.9 -11.9 -21.8 -15.4 -24.0 72 34 58 53 57 52 56 2038.6 3181.0 1893.6 2214.1 2364.7 2529.8 2853.4

·2·
932.7 477.3 698.5 407.9 523.1 214.7 419.5 25 26 27 28 29 30 31

乌鲁木齐

以上数据保存在文件 examp01_01.xls 中,下面根据以上 31 组观测数据研究年平均气温和全年日照时 数之间的关系。

1.1 数据的散点图
令 x 表示年平均气温, y 表示全年日照时数。由于 x 和 y 均为一维变量,可以先从 x 和 y 的散点图 上直观的观察它们之间的关系,然后再作进一步的分析。 通过以下命令从文件 examp01_01.xls 中读取变量 x 和 y 的数据,然后作出 x 和 y 的观测数据的散点 图(如图 1-1 所示),并求出 x 和 y 的线性相关系数。
>> >> >> >> >> >> ClimateData = xlsread('examp01_01.xls'); % 从 Excel 文件读取数据 x = ClimateData(:, 1); % 提取 ClimateData 的第 1 列,即年平均气温数据 y = ClimateData(:, 5); % 提取 ClimateData 的第 5 列,即全年日照时数数据 plot(x, y, 'k.', 'Markersize', 15) % 绘制 x 和 y 的散点图 xlabel('年平均气温(x)') % 给 X 轴加标签 ylabel('全年日照时数(y)') % 给 Y 轴加标签 %计算 x 和 y 的线性相关系数矩阵 R

>> R = corrcoef(x, y) R = 1.0000 -0.7095 -0.7095 1.0000

3500

3000

2500

全年日照时数(y)

2000

1500

1000

500

6

8

10

12

14

16

18

20

22

24

26

年平均气温 (x)

图 1-1 年平均气温 x 与全年日照时数 y 的散点图

从散点图上看,有 4 组数据有些异常,它们分别是拉萨(9.8, 3181)、重庆(19, 856.2)、成都(16.8, 935.6)

MATLAB 技术论坛电子期刊 第 3 期

·3·

和贵阳(14.9, 1014.8). 其中拉萨的全年日照时数最多,重庆、成都和贵阳的全年日照时数较少。除这 4 组 数据外,散点图表明 x 和 y 的线性趋势比较明显,可以用直线 y = β 0 + β1 x 进行拟合。 从相关系数来看, x 和 y 的线性相关系数为 ?0.7095 ,表明 x 和 y 负相关,这是一个非常有意思的 现象,全年日照时数越多的地方其年平均气温倒越低。

1.2 调用 regress 函数作一元线性回归分析
1.2.1.regress 函数的用法 MATLAB 统计工具箱中提供了 regress 函数,用来作多重线性或广义线性回归分析。对于可控变量 x1 , x2 , , x p 和随机变量 y 的 n 次独立的观测 ( xi1 , xi 2 , , xip ; yi ), i = 1, 2, , n , y 关于 x1 , x2 , , x p 的 p 重广义线性回归模型如下:

? y1 ? ? f1 ( x11 ) ? ? ? f (x ) ? y2 ? = ? 1 21 ? ? ? ? ? ? ? ? yn ? ? f1 ( xn1 )
y

f 2 ( x12 ) f 2 ( x22 ) f 2 ( xn 2 )
X

f p ( x1 p ) ? ? β1 ? ? ε1 ? ?? ? f p ( x2 p ) ? ? β 2 ? ? ε 2 ? + ? ?. ?? ? ? ? ?? ? ? ? f p ( xnp ) ? ? β p ? ? ε n ? ?? ?
β ε

(1.1)

通常假定 ε1 , ε 2 , 量。

, ε n ~ N (0, σ 2 ) ,这里 iid 表示独立同分布。(1.1)式中的 y 为因变量观测值向量, X 为 设计矩阵, f1 , f 2 , , f p 为 p 个函数,对应模型中的 p 项, β 为需要估计的系数向量, ε 为随机误差向
不 同 的 函 数 f1 , f 2 ,

iid

, f p 对 应 不 同 类 型 的 回 归 模 型 , 特 别 地 , 当 f1 ( xi1 ) = xi1 , f 2 ( xi 2 ) = xi 2 ,

,

f p ( xip ) = xip , i = 1,

, n 时,(1.1)式称为 p 重线性回归模型。一元线性回归是多重线性回归的特殊情

况。当模型中需要常数项时,设计矩阵 X 中应有 1 列 1 元素(通常是 X 的第一列)。 regress 函数的调用方式如下:
b = regress(y,X) [b,bint] = regress(y,X) [b,bint,r] = regress(y,X) [b,bint,r,rint] = regress(y,X) [b,bint,r,rint,stats] = regress(y,X) […] = regress(y,X,alpha)

(1)b = regress(y,X) 返回多重线性回归方程中系数向量 β 的估计值 b,这里的 b 为一个 p × 1 的向量。输入参数 y 为因变 量的观测值向量,是 n ×1 的列向量。X 为 n × p 的设计矩阵。regress 函数把 y 或 X 中的不确定数据 NaN 作为缺失数据而忽略它们。 (2)[b,bint] = regress(y,X) 还返回系数估计值的 95%置信区间 bint,它是一个 p × 2 的矩阵,第 1 列为置信下限,第 2 列为置信 上限。 (3)[b,bint,r] = regress(y,X) ? 还返回残差(因变量的真实值 yi 减去估计值 yi )向量 r,它是一个 n ×1 的列向量。 (4)[b,bint,r,rint] = regress(y,X) 还返回残差的 95%置信区间 rint,它是一个 n × 2 的矩阵,第 1 列为置信下限,第 2 列为置信上限。 rint 可用于异常值(或离群值)的诊断,若第 i 组观测的残差的置信区间不包括 0,则可认为第 i 组观测值 为异常值。 (5)[b,bint,r,rint,stats] = regress(y,X) 还返回一个 1× 4 的向量 stats,其元素依次为判定系数 R 、 F 统计量的观测值、检验的 p 值和误差
2

2010 年 11 月
2 2

Matlab 数据拟合——回归分析

·4·

? 方差 σ 的估计值 σ 。 (6)[…] = regress(y,X,alpha) 用 alpha 指定计算 bint 和 rint 时的置信水平为 100(1-alpha)% .

注意:当回归模型中需要常数项时,矩阵 X 中应有 1 列 1 元素。当需要计算判定系数 R2 、 F 统
计量的观测值、 p 值时,模型中应包含常数项。若模型中不包含常数项,regress 函数输出的判定系数

R 2 、 F 统计量的观测值、 p 值是不正确的。在不考虑常数项的情况下,计算出的判定系数 R 2 的值可能
是负的,说明所用模型不适合用户的数据。 1.2.2.调用 regress 函数作一元线性回归分析 下面调用 regress 函数对 y 和 x 作一元线性回归分析。
>> >> >> >> >> >> ClimateData = xlsread('examp01_01.xls'); % 从 Excel 文件读取数据 x = ClimateData(:, 1); % 提取 ClimateData 的第 1 列,即年平均气温数据 y = ClimateData(:, 5); % 提取 ClimateData 的第 5 列,即全年日照时数数据 xdata = [ones(size(x, 1), 1), x]; % 在原始数据 x 的左边加一列 1,即模型包含常数项 [b, bint, r, rint, s] = regress(y, xdata); % 调用 regress 函数作一元线性回归 yhat = xdata*b; % 计算 y 的估计值

% 定义元胞数组,以元胞数组形式显示系数的估计值和估计值的 95%置信区间 >> head1 = {'系数的估计值','估计值的 95%置信下限','估计值的 95%置信上限'}; >> [head1; num2cell([b, bint])] ans = '系数的估计值' [ 3.1154e+003] [ -76.9616] '估计值的 95%置信下限' [ 2.6592e+003] [ -105.9970] '估计值的 95%置信上限' [ 3.5716e+003] [ -47.9261]

% 定义元胞数组,以元胞数组形式显示 y 的真实值、y 的估计值、残差和残差的 95%置信区间 >> head2 = {'y 的真实值','y 的估计值','残差','残差的 95%置信下限','残差的 95%置信上限'}; % 同时显示 y 的真实值、y 的估计值、残差和残差的 95%置信区间 >> [head2; num2cell([y, yhat, r, rint])] ans = 'y 的真实值' [2.3511e+003] [2.1654e+003] [2.1677e+003] [2.1746e+003] [2.6478e+003] [2.3609e+003] [2.5336e+003] [2.3592e+003] [1.5222e+003] [1.6803e+003] [1.4729e+003] [1.8146e+003] [1.5438e+003] [ 2102] [1.8198e+003] [1.7472e+003] [1.9342e+003] [1.7422e+003]

'y 的估计值' [2.0379e+003] [2.0687e+003] [1.9686e+003] [2.2380e+003] [2.4227e+003] [2.4227e+003] [2.5228e+003] [2.6074e+003] [1.6916e+003] [1.7762e+003] [1.6993e+003] [1.7762e+003] [1.4992e+003] [1.6377e+003] [1.9610e+003] [1.8840e+003] [1.6839e+003] [1.6685e+003]

'残差' '残差的 95%置信下限' [ 313.1847] [ -461.6917] [ 96.7001] [ -686.1726] [ 199.0501] [ -581.9400] [ -63.4154] [ -840.7773] [ 225.0768] [ -534.8155] [ -61.8232] [ -826.3056] [ 10.8268] [ -744.1662] [-248.2309] [ -987.0492] [-169.3882] [ -944.3372] [ -95.9459] [ -876.4774] [-226.3844] [ -999.5520] [ 38.3541] [ -742.9172] [ 44.6157] [ -719.2939] [ 464.2849] [ -289.2774] [-141.1537] [ -924.0249] [-136.7921] [ -919.1600] [ 250.3079] [ -520.9526] [ 73.7003] [ -702.2379]

'残差的 95%置信上限' [1.0881e+003] [ 879.5727] [ 980.0403] [ 713.9466] [ 984.9692] [ 702.6593] [ 765.8198] [ 490.5873] [ 605.5608] [ 684.5855] [ 546.7832] [ 819.6253] [ 808.5253] [1.2178e+003] [ 641.7175] [ 645.5757] [1.0216e+003] [ 849.6384]

MATLAB 技术论坛电子期刊 第 3 期
[ 1616] [ 1614] [1.6691e+003] [ 856.2000] [ 935.6000] [1.0148e+003] [2.0386e+003] [ 3181] [1.8936e+003] [2.2141e+003] [2.3647e+003] [2.5298e+003] [2.8534e+003] [1.3299e+003] [1.4453e+003] [1.2606e+003] [1.6531e+003] [1.8224e+003] [1.9686e+003] [1.9148e+003] [2.3612e+003] [1.9148e+003] [2.2611e+003] [2.6459e+003] [2.3150e+003] [2.4612e+003] [ 286.1312] [ 168.6888] [ 408.4966] [-796.9074] [-886.8229] [-953.8499] [ 123.8232] [ 819.8461] [ -21.1768] [ -47.0039] [-281.2117] [ 214.8230] [ 392.1961] [ -451.5246] [ -587.4691] [ -311.0565] [-1.5087e+003] [-1.5907e+003] [-1.6466e+003] [ -659.2486] [ 118.1779] [ -805.6671] [ -823.2940] [-1.0132e+003] [ -553.8992] [ -353.8710]

·5·
[1.0238e+003] [ 924.8467] [1.1280e+003] [ -85.1229] [ -182.9957] [ -261.0701] [ 906.8951] [1.5215e+003] [ 763.3135] [ 729.2863] [ 450.7301] [ 983.5453] [1.1383e+003]

% 定义元胞数组,以元胞数组形式显示判定系数、F 统计量的观测值、检验的 p 值和误差方差的估计值 >> head3 = {'判定系数','F 统计量的观测值','检验的 p 值','误差方差的估计值'}; >> [head3; num2cell(s)] '判定系数' [ 0.5033] 'F 统计量的观测值' [ 29.3884] '检验的 p 值' [7.8739e-006] '误差方差的估计值' [ 1.4689e+005]

从输出的结果看,常数项和回归系数的估计值分别为 3.1154e+003 和 -7.6962e+001,从而可以写出线 性回归方程为

? y = 3115.4 ? 76.962 x . 回归系数估计值的 95%置信区间为 [ ?105.997, ? 47.926] .
对回归直线进行显著性检验,原假设和对立假设分别为

H 0 : β1 = 0,

H1 : β1 ≠ 0.

检验的 p 值为 7.8739e-006 < 0.01,可知在显著性水平 α = 0.01 下应拒绝原假设 H 0 ,可认为 y (全年日 照时数)与 x (年平均气温)的线性关系是显著的。 也可以通过 F 统计量的观测值与临界值 Fα (1, n ? 2) 作比较得出结论,当 F ≥ Fα (1, n ? 2) 时,拒绝 原假设,认为 y 与 x 的线性关系是显著的,否则接受原假设,认为 y 与 x 的线性关系是不显著的。对于 本例, F 统计量的观测值为 29.388,临界值 F0.01 (1, 29) = 7.5977 ,显然在显著性水平 α = 0.01 下应拒 绝原假设 H 0 . 利用下面命令画出原始数据散点与回归直线图,如图 1-2 所示。
>> >> >> >> >> >> plot(x, y, 'k.', 'Markersize', 15) % 画原始数据散点 hold on plot(x, yhat, 'linewidth', 3) % 画回归直线 xlabel('年平均气温(x)') % 给 X 轴加标签 ylabel('全年日照时数(y)') % 给 Y 轴加标签 legend('原始散点','回归直线'); % 加标注框

2010 年 11 月
3500

Matlab 数据拟合——回归分析
Residual Case Order Plot 原始散点 回归直线 1500 1000 500 Residuals 0 -500 -1000

·6·

3000

全年日照时数(y)

2500

2000

1500

1000 -1500 500 6 8 10 12 14 16 18 年平均气温(x) 20 22 24 26 5 10 15 20 Case Number 25 30

图 1-2 原始数据散点与回归直线图

图 1-3 原始数据对应残差图

1.2.3.残差分析 通过对残差和残差的置信区间进行分析,可以看出原始数据中是否存在异常点,若残差的置信区间 不包括 0 点,可认为该组观测为异常点。为了直观,这里用 rcoplot 函数按顺序画出各组观测对应的残差 和残差的置信区间。作出的图如图 1-3 所示。
>> figure % 新建一个图形窗口 >> rcoplot(r,rint) %按顺序画出各组观测对应的残差和残差的置信区间

残差图的横坐标表示观测序号,纵坐标表示残差值的大小。图中的每条竖直线段对应一组观测的残 差和残差的置信区间,线段的中点处的圆圈所对应纵坐标为残差值的大小,线段的上端点的纵坐标为置 信上限,下端点的纵坐标为置信下限。从残差图可以看出有 4 条线段(红色虚线)与水平线 y = 0 没有交 点,它们对应的观测序号分别为 22、23、24 和 26,也就是说这 4 组观测对应的残差的置信区间不包含 0 点,可认为这四组观测数据为异常数据。它们分别是拉萨(9.8, 3181)、重庆(19, 856.2)、成都(16.8, 935.6) 和贵阳(14.9, 1014.8),这和从散点图上直接观察的结果相吻合。 将以上四组异常数据剔除后重新作散点图(如图 1-4 所示),重新计算线性相关系数如下:
>> >> >> >> >> >> >> xt = x(y<3000 & y>1250); % 根据条件 y<3000 & y>1250 剔除异常数据 yt = y(y<3000 & y>1250); figure % 新建一个空的图形窗口 plot(xt, yt, 'ko') % 画剔除异常数据后的散点图 xlabel('年平均气温(x)') % 为 X 轴加标签 ylabel('全年日照时数(y)') % 为 Y 轴加标签 Rt = corrcoef(xt, yt) % 重新计算相关系数矩阵

Rt = 1.0000 -0.8554 -0.8554 1.0000

剔除异常数据后, x 和 y 的线性相关系数变为 ?0.8554 ,线性相关性进一步增强。调用 regress 函数 重新计算,命令及结果如下:
>> xtdata = [ones(size(xt, 1), 1), xt]; % 在数据 xt 的左边加一列 1 % 调用 regress 函数对处理后数据作一元线性回归 >> [b, bint, r, rint, s] = regress(yt, xtdata); >> b % 显示常数项和回归系数的估计值 b = 2.9838e+003 -6.3628e+001

MATLAB 技术论坛电子期刊 第 3 期
>> s s = 7.3178e-001 6.8206e+001 1.3088e-008 4.0508e+004 % 显示判定系数 、 F 统计量的观测值、 p 值和误差方差的估计值

·7·

>> ythat = xtdata*b;

% 重新计算 y 的估计值

? 此时的回归方程为 y = 2983.8 ? 63.628 x . 对回归直线进行显著性检验的 p 值为 1.3088e-008,可知 y (全年日照时数)与 x (年平均气温)的线性关系更为显著。作出原始数据散点、原始数据对应的回 归直线和剔除异常数据后的回归直线图,如图 1-5 所示。
>> >> >> >> >> >> >> >> >> >>
3000 2800 2600 全年日照时数(y) 2400 2200 2000 1800 1600 1400 6 8 10 12 14 16 18 年平均气温(x) 20 22 24 26 1000 全年日照时数(y) 2500

figure; % 新建一个图形窗口 plot(x, y, 'ko'); % 画原始数据散点 hold on; % 图形叠加 [xsort, id1] = sort(x); % 为了画图的需要将 x 从小到大排序 yhatsort = yhat(id1); % 将估计值 yhat 按 x 排序 plot(xsort, yhatsort, 'r--','linewidth',3); % 画原始数据对应的回归直线,红色虚线 plot(xt, ythat, 'linewidth', 3); % 画剔除异常数据后的回归直线,蓝色实线 legend('原始数据散点','原始数据回归直线','剔除异常数据后回归直线') % 为图形加标注框 xlabel('年平均气温(x)'); % 为 X 轴加标签 ylabel('全年日照时数(y)'); % 为 Y 轴加标签
3500

3000

原始数据散点 原始数据回归直线 剔除异常数据后回归直线

2000

1500

500 6

8

10

12

14 16 18 年平均气温(x)

20

22

24

26

图 1-4 剔除异常数据后的散点图

图 1-5 原始数据散点与两条回归直线图

图 1-5 中的圆圈为原始数据点,红色虚线为原始数据对应的回归直线,蓝色实线为剔除异常数据后的 回归直线。由于受异常数据的影响,两次回归结果并不相同。

1.3 调用 regstats 函数作一元线性回归分析
1.3.1.regstats 函数的用法 MATLAB 统计工具箱中提供了 regstats 函数,也可用来作多重线性或广义线性回归分析,它的调用方 式如下:
regstats(y,X,model) stats = regstats(…) stats = regstats(y,X,model,whichstats)

(1)regstats(y,X,model) 作多重线性回归分析。输入参数 X 为自变量观测值矩阵(或设计矩阵),它是 n × p 的矩阵。默认情 况下,regstats 函数自动在 X 第 1 列元素的左边加入一列 1,不需要用户自己添加。输入参数 y 为因变量 的观测值向量,是 n ×1 的列向量。可选的输入参数 model 是一个字符串,用来控制回归模型的类型,其

2010 年 11 月
可用的字符串如表 1-2 所示。

Matlab 数据拟合——回归分析
表 1-2 regstats 函数支持的 model 参数

·8·

model 参数的参数值 'linear' 'interaction' 'quadratic' 'purequadratic'





带有常数项的线性模型(默认情况) 带有常数项、线性项和交叉项的模型 带有常数项、线性项、交叉项和平方项的模型 带有常数项、线性项和平方项的模型

在这种调用方式下,regstats 函数会生成一个交互式图形用户界面(GUI),界面上带有回归诊断统 计量列表,包括系数的估计值、因变量的预测值、残差、判定系数、调整的判定系数、 F 检验和 t 检验的 相关结果等,共 23 个可选项。通过这个界面,用户可以很方便地将回归分析的各种结果导入 MATLAB 工作空间。 (2)stats = regstats(…) 返回一个结构体变量 stats,它有 24 个字段,包括了回归分析的所有诊断统计量。这种调用方式不生 成图形用户界面,stats 的后 23 个字段分别与图形用户界面上的 23 个选项相对应。 (3)stats = regstats(y,X,model,whichstats) 仅返回由 whichstats 参数指定的统计量。whichstats 可以是形如 'leverage' 的单个字符串,也可以是形 如{'leverage' 'standres' 'studres'}的字符串元胞数组。若 whichstats 是字符串 'all' ,则返回所有统计量。

注意:当需要计算 F 统计量的观测值时,模型中应包含常数项。若模型中不包含常数项,regstats
函数输出的 F 统计量的观测值是不正确的。在不考虑常数项的情况下,计算出的判定系数 R 的值可能是 负的,说明所用模型不适合用户的数据。 1.3.2.调用 regstats 函数作一元线性回归分析 为了使结果看起来更加直观,作者编写了 MATLAB 函数 reglm,函数代码如下:
function stats = reglm(y,X,model,varnames) % 多重线性回归分析或广义线性回归分析 % % reglm(y,X),产生线性回归分析的方差分析表和参数估计结果,并以表格形式显示在屏幕上. 参 % 数 X 是自变量观测值矩阵,它是 n 行 p 列的矩阵. y 是因变量观测值向量,它是 n 行 1 列的列向量. % % stats = reglm(y,X),还返回一个包括了回归分析的所有诊断统计量的结构体变量 stats. % % stats = reglm(y,X,model),用可选的 model 参数来控制回归模型的类型. model 是一个字符串, % 其可用的字符串如下 % 'linear' 带有常数项的线性模型(默认情况) % 'interaction' 带有常数项、线性项和交叉项的模型 % 'quadratic' 带有常数项、线性项、交叉项和平方项的模型 % 'purequadratic' 带有常数项、线性项和平方项的模型 % % stats = reglm(y,X,model,varnames),用可选的 varnames 参数指定变量标签. varnames % 可以是字符矩阵或字符串元胞数组,它的每行的字符或每个元胞的字符串是一个变量的标签,它的行 % 数或元胞数应与 X 的列数相同. 默认情况下,用 X1,X2,…作为变量标签. % % 例: % x = [215 250 180 250 180 215 180 215 250 215 215 % 136.5 136.5 136.5 138.5 139.5 138.5 140.5 140.5 140.5 138.5 138.5]'; % y = [6.2 7.5 4.8 5.1 4.6 4.6 2.8 3.1 4.3 4.9 4.1]'; % reglm(y,x,'quadratic') % % ----------------------------------方差分析表---------------------------------2

MATLAB 技术论坛电子期刊 第 3 期
% % % % % % % % % % % % % % % % % % % 方差来源 回归 残差 总计 自由度 5.0000 5.0000 10.0000 平方和 15.0277 1.9742 17.0018 0.6284 4.7273 均方 3.0055 0.3948 F值 7.6122

·9·
p值 0.0219

均方根误差(Root MSE) 因变量均值(Dependent Mean)

判定系数(R-Square) 调整的判定系数(Adj R-Sq)

0.8839 0.7678

-----------------------------------参数估计----------------------------------变量 估计值 标准误 t值 p值 常数项 30.9428 2011.1117 0.0154 0.9883 X1 0.7040 0.6405 1.0992 0.3218 X2 -0.8487 29.1537 -0.0291 0.9779 X1*X2 -0.0058 0.0044 -1.3132 0.2461 X1*X1 0.0003 0.0003 0.8384 0.4400 X2*X2 0.0052 0.1055 0.0492 0.9626 Copyright 2009 - 2010 xiezhh. $Revision: 1.0.0.0 $ $Date: 2009/12/22 21:41:00 $

if nargin < 2 error('至少需要两个输入参数'); end p = size(X,2); % X 的列数,即变量个数 if nargin < 3 || isempty(model) model = 'linear'; % model 参数的默认值 end % 生成变量标签 varnames if nargin < 4 || isempty(varnames) varname1 = strcat({'X'},num2str([1:p]')); varnames = makevarnames(varname1,model); % 默认的变量标签 else if ischar(varnames) varname1 = cellstr(varnames); elseif iscell(varnames) varname1 = varnames(:); else error('varnames 必须是字符矩阵或字符串元胞数组'); end if size(varname1,1) ~= p error('变量标签数与 X 的列数不一致'); else varnames = makevarnames(varname1,model); % 指定的变量标签 end end ST = regstats(y,X,model); % 调用 regstats 函数进行线性回归分析,返回结构体变量 ST f = ST.fstat; % F 检验相关结果 t = ST.tstat; % t 检验相关结果 % 显示方差分析表 fprintf('\n');

2010 年 11 月

Matlab 数据拟合——回归分析

·10·

fprintf('------------------------------方差分析表------------------------------'); fprintf('\n'); fprintf('%s%7s%15s%15s%15s%12s','方差来源','自由度','平方和','均方','F 值','p 值'); fprintf('\n'); fmt = '%s%13.4f%17.4f%17.4f%16.4f%12.4f'; fprintf(fmt,'回归',f.dfr,f.ssr,f.ssr/f.dfr,f.f,f.pval); fprintf('\n'); fmt = '%s%13.4f%17.4f%17.4f'; fprintf(fmt,'残差',f.dfe,f.sse,f.sse/f.dfe); fprintf('\n'); fmt = '%s%13.4f%17.4f'; fprintf(fmt,'总计',f.dfe+f.dfr,f.sse+f.ssr); fprintf('\n'); fprintf('\n'); % 显示判定系数等统计量 fmt = '%22s%15.4f%25s%10.4f'; fprintf(fmt,'均方根误差(Root MSE)',sqrt(ST.mse),'判定系数(R-Square)',ST.rsquare); fprintf('\n'); fprintf(fmt,'因变量均值(Dependent Mean)',mean(y),'调整的判定系数(Adj R-Sq)',... ST.adjrsquare); fprintf('\n'); fprintf('\n'); % 显示参数估计及 t 检验相关结果 fprintf('-------------------------------参数估计-------------------------------'); fprintf('\n'); fprintf('%8s%18s%15s%15s%12s','变量','估计值','标准误','t 值','p 值'); fprintf('\n'); for i = 1:size(t.beta,1) if i == 1 fmt = '%8s%20.4f%17.4f%17.4f%12.4f\n'; fprintf(fmt,'常数项',t.beta(i),t.se(i),t.t(i),t.pval(i)); else fmt = '%10s%20.4f%17.4f%17.4f%12.4f\n'; fprintf(fmt,varnames{i-1},t.beta(i),t.se(i),t.t(i),t.pval(i)); end end if nargout == 1 stats = ST; end

% 返回一个包括了回归分析的所有诊断统计量的结构体变量

% -----------------子函数----------------------function varnames = makevarnames(varname1,model) % 生成指定模型的变量标签 p = size(varname1,1); varname2 = []; for i = 1:p-1 varname2 = [varname2;strcat(varname1(i),'*',varname1(i+1:end))]; end varname3 = strcat(varname1,'*',varname1); switch model

MATLAB 技术论坛电子期刊 第 3 期
case 'linear' varnames = varname1; case 'interaction' varnames = [varname1;varname2]; case 'quadratic' varnames = [varname1;varname2;varname3]; case 'purequadratic' varnames = [varname1;varname3]; end

·11·

该函数调用了 regstats 函数,用来作一元或多重线性回归分析,并在屏幕上显示回归分析的方差分析 表和参数估计表。代码中的注释部分给出了该函数的 4 种调用方法,这里不再重述。下面仍从文件 examp01_01.xls 中读取 x (年平均气温)和 y (全年日照时数)的数据,然后调用 reglm 函数作一元线性 回归分析。
>> >> >> >> >> ClimateData = xlsread('examp01_01.xls'); % 从 Excel 文件读取数据 x = ClimateData(:, 1); % 提取 ClimateData 的第 1 列,即年平均气温数据 y = ClimateData(:, 5); % 提取 ClimateData 的第 5 列,即全年日照时数数据 varname = 'x'; % 变量标签 reglm(y,x,[],varname); % 调用 reglm 函数作一元线性回归分析

------------------------------------方差分析表-----------------------------------方差来源 自由度 平方和 均方 F值 p值 回归 1.0000 4316961.1182 4316961.1182 29.3884 0.0000 残差 29.0000 4259914.3811 146893.5993 总计 30.0000 8576875.4994 均方根误差(Root MSE) 因变量均值(Dependent Mean) 383.2670 1965.1742 判定系数(R-Square) 调整的判定系数(Adj R-Sq) 0.5033 0.4862

-------------------------------参数估计------------------------------变量 估计值 标准误 t值 p值 常数项 3115.3773 223.0588 13.9666 0.0000 x -76.9616 14.1967 -5.4211 0.0000

可以看出 reglm 函数显示的结果与 regress 函数返回的结果是一致的,只是 reglm 函数显示的结果更加 直观。由于方差分析表中的 p 值小于 0.0001,说明整个回归方程是极显著的。根据参数估计结果可以写

? 出所求的回归方程 y = 3115.3773 ? 76.9616 x ,对常数项和线性项进行的 t 检验的 p 值均小于 0.0001, 说明在回归方程中常数项和线性项均是极显著的。

1.4 调用 robustfit 函数作稳健回归
regress 函数和 regstats 函数利用普通最小二乘法估计模型中的参数,参数的估计值受异常值的影响比 较大。robustfit 函数采用加权最小二乘法估计模型中的参数,受异常值的影响就比较小。robustfit 函数用 来作稳健的多重线性或广义线性回归分析,下面介绍 robustfit 函数的用法。 1.4.1.robustfit 函数的用法 robustfit 函数有以下几种调用方式:
b = robustfit(X,y) b = robustfit(X,y,wfun,tune) b = robustfit(X,y,wfun,tune,const) [b,stats] = robustfit(…)

(1)b = robustfit(X,y) 返回多重线性回归方程中系数向量 β 的估计值 b,这里的 b 为一个 p × 1 的向量。输入参数 X 为自变

2010 年 11 月

Matlab 数据拟合——回归分析

·12·

量观测值矩阵(或设计矩阵),它是 n × p 的矩阵。与 regress 函数不同的是,默认情况下,robustfit 函数 自动在 X 第 1 列元素的左边加入一列 1,不需要用户自己添加。输入参数 y 为因变量的观测值向量,是 n ×1 的列向量。robustfit 函数把 y 或 X 中不确定数据 NaN 作为缺失数据而忽略它们。 (2)b = robustfit(X,y,wfun,tune) 用参数 wfun 指定加权函数,用参数 tune 指定调节常数。wfun 为字符串,其可能的取值如表 1-3 所 示。
表 1-3 robustfit 函数支持的加权函数 加权函数(wfun) 'andrews' 函数表达式 默认调节常数值 1.339

w=

'bisquare'(默认值) 'cauchy' 'fair' 'huber' 'logistic' 'ols' 'talwar' 'welsch'

sin r ? I (| r |< π ) r w = (1 ? r 2 )2 ? I (| r |< 1)
2

4.685 2.385 1.400 1.345 1.205 无 2.795 2.985

w = 1 (1 + r ) w = 1 (1+ | r |) w = 1 max(1, | r |) w = tanh(r ) r
普通最小二乘,无加权函数

w = I (| r |< 1)

w=e

?r2

若调用时没有指定调节常数 tune,则用表 1-3 中列出的默认调节常数值进行计算。表 1-3 中加权函数 中的 r 通过下式计算

r=

resid tune × s × 1-h

其中 resid 为上一步迭代的残差向量,tune 为调节常数,h 是由最小二乘拟合得到的中心化杠杆值向量,s 为误差项的标准差的估计。s 的计算公式为:s = MAD/0.6745,其中 MAD 为残差绝对值的中位数,在正 态分布下,这个估计是无偏的。若 X 中有 p 列,计算 MAD 时,将残差绝对值向量的前 p 个最小值舍去。 用户可以定义自己的权重函数,函数的输入必须是残差向量,输出是权重向量。在调用 robustfit 函数 时,把自定义权重函数的句柄(形如@myfun)作为 wfun 参数传递给 robustfit 函数,此时必须指定 tune 参数。 (3)b = robustfit(X,y,wfun,tune,const) 用参数 const 来控制模型中是否包含常数项。若 const 取值为 'on' 或 1,则模型中包含常数项,此时自 动在 X 第 1 列的左边加入一列 1,若 const 取值为 'off' 或 0,则模型中不包含常数项,此时不改变 X 的 值。 (4)[b,stats] = robustfit(…) 返回一个结构体变量 stats,它的字段包含了用于模型诊断的统计量。stats 有以下字段: ? stats.ols_s — 普通最小二乘法得出的 σ 的估计(RMSE); ? stats.robust_s — σ 的稳健估计; ? stats.mad_s — 用残差绝对值的中位数计算 σ 的估计; ? stats.s — σ 的最终估计,是 ols_s 和 robust_s 的加权平均与 robust_s 中的最大值; ? stats.se — 系数估计的标准误差; ? stats.t — b 与 stats.se 的比值; ? stats.p — t 检验的 p 值; ? stats.covb — 系数向量的协方差矩阵的估计; ? stats.coeffcorr — 系数向量的相关系数矩阵的估计; ? stats.w — 稳健拟合的权重向量; ? stats.h — 最小二乘拟合的中心化杠杆值向量; ? stats.dfe — 误差的自由度;

MATLAB 技术论坛电子期刊 第 3 期

·13·

? stats.R

— 矩阵 X 的 QR 分解中的 R 因子。

1.4.2.调用 robustfit 函数作稳健回归 仍从文件 examp01_01.xls 中读取 x (年平均气温)和 y (全年日照时数)的数据,然后利用 robustfit 函数的第 4 种调用方式作稳健回归,命令及结果如下:
>> ClimateData = xlsread('examp01_01.xls'); % 从 Excel 文件读取数据 >> x = ClimateData(:, 1); % 提取 ClimateData 的第 1 列,即年平均气温数据 >> y = ClimateData(:, 5); % 提取 ClimateData 的第 5 列,即全年日照时数数据 % 利用 robustfit 函数作稳健回归,返回系数的估计值 b 和相关统计量 stats >> [b, stats] = robustfit(x, y) b = 1.0e+003 * 3.0348 -0.0683 stats = ols_s: 383.2670 robust_s: 302.4368 mad_s: 290.4914 s: 312.7337 resid: [31x1 double] rstud: [31x1 double] se: [2x1 double] covb: [2x2 double] coeffcorr: [2x2 double] t: [2x1 double] p: [2x1 double] w: [31x1 double] R: [2x2 double] dfe: 29 h: [31x1 double] >> stats.p % 查看常数项和回归系数的 t 检验的 p 值 ans = 1.0e-005 * 0.0000 0.2119

? 稳健回归得出的回归方程为 y = 3034.8 ? 68.3 x . 常数项和回归系数的 t 检验的 p 值均小于 0.0001, 可知线性回归是显著的。 结构体变量 stats 的 resid 字段的值为残差向量,w 字段的值为相应的权重向量,下面作出残差向量和 权重向量的散点图,从直观上了解残差和权重的关系。
>> >> >> >> figure; % 新建一个图形窗口 plot(stats.resid, stats.w, 'o') xlabel('残差') % 为 X 轴加标签 ylabel('权重') % 为 Y 轴加标签 % 绘制残差和权重的散点图

以上命令作出的图如图 1-6 所示。从图上可以看出残差绝对值越大,其权重就越小,残差为零时,相 应的权重为 1,这就保证了拟合的稳健性。也可以通过对比 regress 和 robustfit 函数的拟合效果,看出加权 最小二乘拟合的稳健性。通过下面的命令作出 regress 和 robustfit 函数的拟合效果图,如图 1-7 所示。

2010 年 11 月
>> >> >> >> >>

Matlab 数据拟合——回归分析

·14·

[xsort, id] = sort(x); % 为了后面画图的需要,将 x 从小到大排序 ysort = y(id); % 将 y 按与 x 相同的方式排序 xdata = [ones(size(xsort, 1), 1), xsort]; % 在原始数据 x 的左边加一列 1 b1 = regress(ysort, xdata); % 调用 regress 函数作一元线性回归 yhat1 = xdata*b1; % 求 regress 函数对应的 y 的估计值

>> b2 = robustfit(xsort, ysort); % 调用 robustfit 函数作稳健的一元线性回归 >> yhat2 = xdata*b2; % 求 robustfit 函数对应的 y 的估计值 >> plot(x, y, 'ko') % 画原始数据散点 >> hold on % 图形叠加 >> plot(xsort, yhat1, 'r--','linewidth',3) % 画 regress 函数对应的回归直线,红色虚线 >> plot(xsort, yhat2, 'linewidth', 3) % 画 robustfit 函数对应的回归直线,蓝色实线 % 为图形加标注框 >> legend('原始数据散点','regress 函数对应的回归直线','robustfit 函数对应的回归直线'); >> xlabel('年平均气温(x)') % 为 x 轴加标签 >> ylabel('全年日照时数(y)') % 为 y 轴加标签
1 0.9 3000 0.8 全年日照时数(y) 0.7 权重 0.6 0.5 0.4 0.3 1000 0.2 0.1 -1500 -1000 -500 残差 0 500 1000 500 6 8 10 12 14 16 18 年平均气温(x) 20 22 24 26 2500 3500

原始数据散点 regress函数对应的回归直线 robustfit函数对应的回归直线

2000

1500

图 1-6 残差和权重的散点图

图 1-7 原始数据散点与 regress、robustfit 回归直线图

regress 和 robustfit 函数的主要区别就在于,regress 函数是基于普通最小二乘拟合,而 robustfit 函数是 基于加权最小二乘拟合。从图 1-7 可以看出通过加权可以消除异常值的影响,增强拟合的稳健性。

2 一元非线性回归分析
头围(head circumference)是反映婴幼儿大脑和颅骨发育程度的重要指标之一,对头围的研究具有非常 重要的意义。作者研究了天津地区 1281 位儿童(700 个男孩,581 个女孩)的颅脑发育情况,测量了年 龄、头宽、头长、头宽/头长、头围和颅围等指标,测量方法:读取头颅 CT 图像数据,根据自编程序自 动测量。测量得到 1281 组数据,年龄跨度从 7 个星期到 16 周岁,数据保存在文件 examp01_02.xls 中,数 据格式如表 1-4 所示。
表 1-4 天津地区 1281 位儿童的颅脑发育情况指标数据(只列出部分数据) 序号 性 别 1 2 3 4 5 m m m m m 年龄及 标识 11Y 20M 10Y 3Y 3Y 年 龄 月 龄 头 宽 头 长 头宽/头长 头 围 颅 围

单位:岁 11 1.666667 10 3 3

单位:月 132 20 120 36 36

单位:毫米 136.0476 149.9043 144.4456 145.7053 139.8267

单位:毫米 168.7998 161.2416 156.6227 163.761 153.2635 0.805970149 0.9296875 0.922252011 0.88974359 0.912328767

单位:厘米 50.90952 50.4282 51.35181 50.27417 48.52064

单位:厘米 48.3008 49.01562 48.14725 48.73305 46.925

MATLAB 技术论坛电子期刊 第 3 期

·15·

1277 1278 1279 1280 1281

f f f f f

17M 5Y 3Y 13M 5Y

1.416667 5 3 1.083333 5

17 60 36 13 60

147.8048 144.4456 150.7441 129.3292 146.5451

140.2466 162.0814 145.7053 143.1859 157.8824

1.053892216 0.89119171 1.034582133 0.903225806 0.928191489

46.52105 49.56883 47.0336 44.99825 49.65208

45.54998 48.48535 46.02226 43.32917 47.91818

注:年龄数据中的 Y 表示年,M 表示月,W 表示星期,D 表示天。性别数据中的 m 表示男性,f 表 示女性。 下面根据这 1281 组数据建立头围关于年龄的回归方程。

2.1 数据的散点图
令 x 表示年龄, y 表示头围。 x 和 y 均为一维变量,同样可以先从 x 和 y 的散点图上直观的观察它 们之间的关系,然后再作进一步的分析。 通过以下命令从文件 examp01_02.xls 中读取变量 x 和 y 的数据,然后作出 x 和 y 的观测数据的散点 图,如图 1-8 所示。
>> >> >> >> >> >> HeadData = xlsread('examp01_02.xls'); %从 Excel 文件读取数据 x = HeadData(:, 4); % 提取 HeadData 矩阵的第 4 列数据,即年龄数据 y = HeadData(:, 9); % 提取 HeadData 矩阵的第 9 列数据,即头围数据 plot(x, y, 'k.') % 绘制 x 和 y 的散点图 xlabel('年龄(x)') % 为 X 轴加标签 ylabel('头围(y)') % 为 Y 轴加标签
60

55

50 头围(y) 45 40 35 0

2

4

6

8 年龄(x)

10

12

14

16

图 1-8 年龄与头围数据的散点图

从图 1-8 可以看出 y (头围)和 x (年龄)之间呈现非线性相关关系,可以考虑作非线性回归。根据 散点图的走势,可以选取以下函数作为理论回归方程:

? 负指数函数: y = β1e , x + β1 ? 双曲线函数: y = , β 2 x + β3 ? 幂函数: y = β1 ( x + β 2 ) β3 ,

x + β3

β2

2010 年 11 月
? Logistic 曲线函数: y =

Matlab 数据拟合——回归分析

·16·

1 + β 2 e ? ( x + β3 ) ? 对数函数: y = β1 + β 2 ln( x + β3 ) .

β1



以上函数中都包含有 3 个未知参数 β1 、 β 2 和 β 3 ,需要由观测数据进行估计,根据需要还可以减少 或增加未知参数的个数。以上函数都可以呈现出先急速增加,然后趋于平缓的趋势,比较适合头围和年 龄的观测数据,均可以作为备选的理论回归方程。

2.2 调用 nlinfit 函数作一元非线性回归分析
2.2.1.nlinfit 函数的用法 MATLAB 统计工具箱中的 nlinfit 函数用来作一元或多重非线性回归。多重非线性回归模型的一般形 式如下:

yi = f
通常假定 ε1 , ε 2 , 记
iid

(( β , β ,
1 2

β k ) ; xi1 , xi 2 , , xip ) + ε i , i = 1, 2, β k ) ; x1 , x2 , , x p ) .

n.

(1.2)

, ε n ~ N (0, σ 2 ) . 由(1.2)式可得理论回归方程为

? y= f

(( β , β ,
1 2

(1.3)

x1 p ? ? x11 x12 ? y1 ? ? β1 ? ? ? ? ? ? ? x x22 x2 p ? β ? y2 ? , X = ? 21 (1.4) Y= , β = ? 2 ?. ? ? ? ? ? ? ? ? ? ? ? ? ?x xnp ? ? yn ? ? βk ? ? n1 xn 2 ? Y 为因变量观测值向量。矩阵 X 为自变量观测值矩阵, X 的每一列对应一个变量,每一行对应一组观 测。 β 为模型的未知参数向量(可以为列向量或行向量)。
nlinfit 函数的调用方法如下:
beta = nlinfit(X,y,fun,beta0) [beta,r,J,COVB,mse] = nlinfit(X,y,fun,beta0) […] = nlinfit(X,y,fun,beta0,options)

(1)beta = nlinfit(X,y,fun,beta0) 返回非线性回归方程(1.3)式中未知参数的估计值 beta,若回归方程中有 k 个未知参数,则 beta 为包 含 k 个元素的向量。输入参数 X 为自变量的观测值矩阵,是 n × p 的矩阵,形如(1.4)式。输入参数 y 为因 变量的观测值向量,是 n ×1 的列向量,形如(1.4)式。输入参数 fun 是回归方程(1.3)式所对应函数的句柄或 函数名,应根据回归方程(1.3)式编写如下形式的函数:
yhat = modelfun(b, X) 函数体

其中 modelfun 为函数名,b 为未知参数向量。nlinfit 函数的输入参数 beta0 为用户设定的未知参数的初 值,不同的初值可能会有不同的估计结果,设定初值时应慎重,最好能根据实际问题提前有个预判。 nlinfit 函数把 y 或 modelfun(beta0,X)中的不确定数据 NaN 作为缺失数据而忽略。 (2)[beta,r,J,COVB,mse] = nlinfit(X,y,fun,beta0) 还返回残差向量 r,雅可比矩阵 J,未知参数的协方差矩阵的估计 COVB,误差方差 σ 的估计 mse (均方误差平方和)。 这里的输出可作为 nlparci 函数的输入,用来计算参数估计值的置信区间。也可作为 nlpredci 函数的 输入,用来计算给定 x 处的预测值及预测值的置信区间。结合案例,后面会介绍 nlparci 和 nlpredci 函数 的用法。 (3)[…] = nlinfit(X,y,fun,beta0,options) nlinfit 函数利用麦夸特(Levenberg-Marquardt)算法求解非线性回归方程中的未知参数,它是一种迭 代算法,这里用参数 options 设定该算法的控制参数。options 是由 statset 函数生成的结构体变量,它可用
2

MATLAB 技术论坛电子期刊 第 3 期 的字段及字段值如表 1-5 所示。
表 1-5 参数 options 的可用字段及字段值 字段名 Display 字段值 'off' (the default) 'iter' 'final' MaxIter TolFun 正整数 正实数 说 明

·17·

设置显示中间结果的方式。默认值是'off',表示不显示, 'iter'表示分步显示,'final'表示只显示最终结果 设置最大迭代次数,默认值为 100. 当迭代次数超过最大迭代次数时, 停止迭代 设置残差平方和的终止容限,默认值为 1e8. 当相邻两步迭代的残差平 方和差的绝对值小于或等于终止容限与当前残差平方和的乘积时,停 止迭代

TolX DerivStep FunValCheck Robust WgtFun

正实数 正的标量或与 beta0 等长的向量 'off' 或 'on' 'off' 或 'on' 'bisquare' 'andrews' 'cauchy' 'fair' 'huber' 'logistic' 'talwar' 'welsch' 自定义函数句柄

设置参数估计值的终止容限,默认值为 1e8. 设置有限差分梯度计算的步长,默认值是 eps^(1/3) 默认值是 'on' ,检查目标函数的无效值,例如 NaN 或 Inf 默认值是 'off',若为 'on',表示调用稳健拟合方法 当 Robust 字段值为 'on' 时,WgtFun 字段才是有效的。用来设置加权函 数,默认值是 'bisquare'. 若用户自定义加权函数,函数的输入是标准化残差向量,输出是权重 向量

Tune

正实数

设置稳健拟合的调节常数。默认值与加权函数有关,如表 1-3 所示。 若用户自定义加权函数,此时必须指定 tune 参数

注意:稳健的非线性拟合与稳健的线性拟合所用的最小二乘算法是不同的,但是调整权重的方式
是相同的。特别地,它们所用的加权函数和相应的调节常数都是相同的。 2.2.2.选择合适的理论回归方程 在调用 nlinfit 函数作一元非线性回归分析之前,应根据观测数据的特点选择合适的理论回归方程, 理论回归方程往往是不唯一的,可以有多种选择。2.1 节中列出了 5 个备选的理论回归方程,这里选择负 指数函数 y = β1e
x + β3

β2

作为理论回归方程,当然用户也可以选择其他函数。

有了理论回归方程之后,首先编写理论回归方程所对应的 M 函数。函数应有两个输入参数,一个输 出参数。第 1 个输入为未知参数向量,对于一元回归,第 2 个输入为自变量观测值向量,而对于多重回 归,第 2 个输入为自变量观测值矩阵。函数的输出为因变量观测值向量。针对所选择的负指数函数,编 写 M 函数如下:
function y = HeadCir1(beta, x) y = beta(1) * exp(beta(2) ./ (x + beta(3)));

将以上 2 行代码写到 MATLAB 程序编辑窗口,保存为 HeadCir1.m 文件,保存路径默认即可。在这 种定义方式下,可以将@HeadCir1 或 'HeadCir1' (注意有引号)传递给 nlinfit 函数,作为 fun 参数。

注意:对于比较简单的理论回归方程,还可以用@符号定义匿名函数。例如可以如下定义负指数
函数:
HeadCir2 = @(beta, x)beta(1)*exp(beta(2)./(x+beta(3)));

与前面的定义不同,这里返回的 HeadCir2 直接就是函数句柄,可以把它作为 fun 参数传递给 nlinfit

2010 年 11 月
函数。

Matlab 数据拟合——回归分析

·18·

2.2.3.调用 nlinfit 函数作一元非线性回归分析 下面调用 nlinfit 函数对 y (头围)和 x (年龄)作稳健的一元非线性回归分析。
>> HeadData = xlsread('examp01_02.xls'); % 从 Excel 文件读取数据 >> HeadData = sortrows(HeadData, 4); % 为作图需要,将 HeadData 按第 4 列从小到大排序 >> x = HeadData(:, 4); % 提取排序后年龄数据 >> y = HeadData(:, 9); % 提取排序后头围数据 >> options = statset; % 创建一个结构体变量,用来设定迭代算法的控制参数 >> options.Robust = 'on'; % 调用稳健拟合方法 % 调用 nlinfit 函数作非线性回归 >> [beta,r,J,COVB,mse] = nlinfit(x,y,@HeadCir1,[53,-0.2604,0.6276], options); %或[beta,r,J,COVB,mse] = nlinfit(x,y,'HeadCir1',[53,-0.2604,0.6276],options); %或[beta,r,J,COVB,mse] = nlinfit(x,y,HeadCir2,[53,-0.2604,0.6276],options); >> beta %查看未知参数估计值 beta = 52.3766 >> mse mse = 2.7498 -0.2595 0.7604

% 查看均方残差平方和

为了程序的完整性,上面再次给出了读取数据的命令,并且为了后面作图的需要,将数据按年龄从 小到大进行了排序。上面调用 statset 函数定义了一个结构体变量 options,通过设置 options 的 Robust 字段 值为 'on' 来调用稳健拟合方法。 未知参数初值的选取是一个难点。从散点图 1-8 上可以看到,随着年龄的增长,人的头围也在增长, 但是头围不会一直增长,到了一定年龄之后,头围就稳定在 50 至 55 之间,注意到 lim
β2 β3
x →+∞

β1e

x + β3

β2

= β1 ,可

以选取 β1 的初值为 50 至 55 之间的一个数,不妨选为 53. 再注意到初生婴儿的头围应在 35cm 左右,可得

β 53e = 35 ,从而 2 = ?0.4149 . 从图 1-8 还可看到 2 岁儿童的头围在 48cm 左右,可得 53e 2+ β = 48 , β3 β2 从而 = ?0.0991 . 于是可得 β 2 = ?0.2604, β3 = 0.6276 ,故选取未知参数向量 ( β1 , β 2 , β 3 ) 的初 2 + β3 值为[53, -0.2604, 0.6276]。实际上,在确定 β1 的初值在 50 至 55 之间后, β 2 和 β 3 可以尝试随意指定,例
3

β2

如[50, 1, 1]、[50, -1, 1]都是可以的,对估计结果影响非常小,也就是说初值在一定范围内都是稳定的。 由未知参数的估计值可以写出头围关于年龄的一元非线性回归方程为

? y = 52.3766e
2.2.4.绘制一元非线性回归曲线

?

0.2595 x + 0.7604

.

调用下面的命令作出年龄与头围的散点和头围关于年龄的回归曲线图。
>> >> >> >> yhat = HeadCir1(beta, x); % 求 y 的估计值 figure % 新建一个空的图形窗口 plot(x, y, 'k.') % 绘制 x 和 y 的散点图 hold on

% 绘制 robustfit 函数对应的回归直线,蓝色实线,线宽为 3

MATLAB 技术论坛电子期刊 第 3 期
>> >> >> >> plot(x, yhat, 'linewidth', 3) xlabel('年龄(x)') % 给 X 轴加标签 ylabel('头围(y)') % 给 Y 轴加标签 legend('原始数据散点','非线性回归曲线')

·19·

% 为图形加标注框

以上命令作出的图形如图 1-9 所示,从图上可以看出拟合效果还是很不错的。有些遗憾的是,由于是 非线性回归,这里不能像线性回归那样对回归方程作显著性检验。
60

原始数据散点 非线性回归曲线
55

50

头围(y)
45 40 35 0

2

4

6

8

10

12

14

16

年龄 (x)

图 1-9 年龄与头围的散点和回归曲线图

2.2.5.参数估计值的置信区间 在调用 nlinfit 函数得到未知参数的估计值 beta、残差向量 r、雅可比矩阵 J、未知参数的协方差矩阵 的估计 COVB 和误差方差 σ 的估计 mse 后,还可以调用 nlparci 函数计算参数估计值的置信区间。下面 给出 nlparci 函数的两种调用方式,计算参数估计值的 95%置信区间。
2

% 求参数估计值的 95%置信区间的第 1 种方式 >> ci1 = nlparci(beta, r, 'covar', COVB, 'alpha', 0.05) ci1 = 52.0923 -0.2912 0.6173 52.6609 -0.2278 0.9035

% 求参数估计值的 95%置信区间的第 2 种方式 >> ci1 = nlparci(beta, r, 'jacobian', J, 'alpha', 0.05) ci1 = 52.0935 -0.2911 0.6178 52.6597 -0.2279 0.9029

在 'alpha' 参数及其参数值缺省的情况下,将返回参数估计值的 95%置信区间。需要注意的是,上面 两种方式求出的置信区间稍有不同,当 nlinfit 函数调用稳健拟合方法时,应利用上面的第 1 种方式求参数 估计值的置信区间。 2.2.6.头围平均值的置信区间和观测值的预测区间 对于 x (年龄)的一个给定值 x0 ,相应的 y (头围)是一个随机变量,具有一定的分布。 x 给定时

2010 年 11 月

Matlab 数据拟合——回归分析

·20·

y 的总体均值的区间估计称为平均值(或预测值)的置信区间, y 的观测值的区间估计称为观测值的预
测区间。 求出头围关于年龄的回归曲线后,对于给定的年龄,可以调用 nlpredci 函数求出头围的预测值(即 x 给定时 y 的总体均值)、预测值的置信区间和观测值的预测区间。nlpredci 函数的调用方式如下:
[ypred,delta] = nlpredci(modelfun,x,beta,r,'covar',COVB) [ypred,delta] = nlpredci(modelfun,x,beta,r,'jacobian',J) […] = nlpredci(…,param1,val1,param2,val2,…)

(1)[ypred,delta] = nlpredci(modelfun,x,beta,r,'covar',COVB) 返回 x 处的预测值 ypred 和预测值的 95%置信区间宽度的一半 delta. 输入参数 modelfun 为回归方程 所对应函数的句柄。x 为指定的自变量的取值,可以为标量或数组。beta 为 nlinfit 函数返回的参数估计 值,r 为残差,COVB 为未知参数的协方差矩阵的估计。 (2)[ypred,delta] = nlpredci(modelfun,x,beta,r,'jacobian',J) 求预测值及预测值的 95%置信区间的另外一种方式。J 为 nlinfit 函数返回的雅可比矩阵,其余参数同 上。

注意:当 nlinfit 函数调用稳健拟合方法时,应利用上面的第 1 种调用方式求预测值及预测值的 95%置信区间。 (3)[…] = nlpredci(…,param1,val1,param2,val2,…) 接受可选的参数名/参数值对,可用来求观测值的预测区间。可用的参数名及参数值如表 1-6 所示。
表 1-6 nlpredci 函数支持的参数名与参数值列表 参数名 'alpha' 'mse' 'predopt' 参数值 0 到 1 之间的数 nlinfit 函数返回的均方误差平方和 'curve' 或 'observation' 说 明

设定置信水平为 100(1-alpha)%. 默认值为 0.05 若 nlinfit 函数调用了稳健拟合方法,在求预测区间时,需设定此参数 若取值为'curve',则计算预测值及预测值的置信区间。若为'observation ', 则计算预测值及观测值的预测区间,若此时 nlinfit 函数调用了稳健拟合方 法,需设定 'mse' 参数

'simopt'

'on' 或 'off'

若取值为 'on',表示联合置信限,若为'off',表示非联合置信限

nlpredci 函数会把残差向量 r 或雅可比矩阵 J 中的不确定数据 NaN 作为缺失数据而忽略。 下面调用 nlpredci 函数求 y (头围)的 95%预测区间,并作出回归曲线和预测区间图。
>> xx =[0:0.1:16]'; % 设定年龄向量

% 计算给定年龄处头围预测值和预测区间 >> [ypred,delta] = nlpredci(@HeadCir1,xx,beta,r,'covar',COVB,'mse',mse, ... 'predopt','observation'); >> yup = ypred + delta; >> ydown = ypred - delta; >> >> >> >> % 预测区间上限(线) % 预测区间下限(线)

figure % 新建一个空的图形窗口 hold on h1 = fill([xx; flipud(xx)],[yup; flipud(ydown)],[0.5,0.5,0.5]); % 填充预测区间 set(h1,'EdgeColor','none','FaceAlpha',0.5); % 设置填充区域边界线条颜色和面板透明度 % 画预测区间上限曲线,红色虚线 % 画预测区间下限曲线,蓝色点划线 % 画回归曲线,黑色实线

>> plot(xx,yup,'r--','LineWidth',2); >> plot(xx,ydown,'b-.','LineWidth',2); >> plot(xx, ypred, 'k','linewidth', 2)

>> grid on % 添加辅助网格 >> ylim([32, 57]) % 设置 y 轴的显示范围为 32 至 57 >> xlabel('年龄(x)') % 给 X 轴加标签

MATLAB 技术论坛电子期刊 第 3 期

·21·

>> ylabel('头围(y)') % 给 Y 轴加标签 >> h2 = legend('预测区间','预测区间上限','预测区间下限','回归曲线'); % 为图形加标注框 >> set(h2, 'Location', 'SouthEast') % 设置标注框的放置位置为图形窗口右下角

以上命令作出的图如图 1-10 所示。

图 1-10 头围关于年龄的回归曲线和 95%预测区间

图 1-10 可以作为评价儿童颅脑发育情况的参考图,参照图中给出的各年龄段头围的 95%预测区间, 可以评价儿童颅脑发育是否正常。读者可以尝试对男孩和女孩的数据分别作一元非线性回归,得出更具 实际意义和参考价值的结果。

2.3 利用曲线拟合工具 cftool 作一元非线性拟合
MATLAB 中有一个功能强大的曲线拟合工具箱(Curve Fitting Toolbox),其中提供了 cftool 函数, 用来通过界面操作的方式进行一元数据拟合。在 MATLAB 命令窗口运行 cftool 命令将打开如图 1-11 所示 的曲线拟合主界面,通过该界面可以实现以下功能: · 从 MATLAB 工作空间导入数据; · 进行数据的预处理,如数据筛选和数据平滑处理; · 利用内置的不同类型的模型进行数据拟合,允许用户自定义模型; · 生成相关结果,如参数估计值,置信区间,相关统计量等; · 进行插值、外推、差分和积分等后处理; · 导出结果到 MATLAB 工作空间,以便进行后续的分析和可视化处理。

2010 年 11 月

Matlab 数据拟合——回归分析

·22·

图 1-11

曲线拟合主界面

下面结合头围与年龄数据的拟合介绍曲线拟合界面的用法。 2.3.1.cftool 函数的调用方式 cftool 函数的调用方式比较简单,有以下 3 种:
cftool cftool(xdata,ydata) cftool(xdata,ydata,w)

以上 3 种方式均可打开曲线拟合主界面,其中输入参数 xdata 为自变量观测值向量,ydata 为因变量 观测值向量,w 为权重向量,它们应为等长向量。 2.3.2.导入数据 如果利用 cftool 函数的第 1 种方式打开曲线拟合主界面,则此时曲线拟合主界面的坐标系还是一片空 白,还没有可以分析的数据,应该先从 MATLAB 工作空间导入变量数据。 首先运行下面的命令将变量数据从文件读入 MATLAB 工作空间。
>> >> >> >> HeadData = xlsread('examp01_02.xls'); % 从 Excel 文件读取数据 HeadData = sortrows(HeadData, 4); % 为作图需要,将 HeadData 按第 4 列从小到大排序 x = HeadData(:, 4); % 提取排序后年龄数据 y = HeadData(:, 9); % 提取排序后头围数据

现在 y (头围)和 x (年龄)的数据已经导入 MATLAB 工作空间,此时单击曲线拟合主界面上的 “Data”按钮,打开数据导入界面如图 1-12 所示。

图 1-12 数据导入界面

MATLAB 技术论坛电子期刊 第 3 期

·23·

数据导入界面的上方有两个按钮:“Data Sets”和“Smooth”。 默认情况下“Data Sets”按钮处于 按下状态,此时用户可以从 MATLAB 工作空间选择变量,创建数据集。 单击“X Data:”后的下拉菜单,从 MATLAB 工作空间选择自变量 x,同样的方式选择因变量和权 重向量。选择完成后,在“Data set name:”后面的文本编辑框中输入数据集的名称“HeadAndAge”, 然后单击“Create data set”按钮,就创建了一个名称为“HeadAndAge”的数据集,界面右侧出现数据的 预览。界面下方的“View”、“Rename”和“Delete”按钮分别用来查看、重命名和删除数据集。 2.3.3.数据的平滑处理 单击数据导入界面上方的“Smooth”按钮,出现如图 1-13 所示界面。

图 1-13 数据平滑界面

数据平滑界面调用了 smooth 函数对数据进行平滑处理,读者可以参阅 3.1 节(数据的平滑处理)中 smooth 函数的用法。界面上“Original data set:”后的下拉菜单用来选择已经定义的数据集;“Smoothed data set:”后面的编辑框用来输入平滑处理后数据集的名称,默认情况下是在原始数据集的名称后加上 “(smooth)”;“Method:”后面的下拉菜单用来选择 smooth 函数进行数据平滑所用的方法;“Span:”后 面的编辑框用来输入 smooth 函数的 span 参数的参数值。以上操作完成后,单击“Create smoothed data set”按钮,即创建了一个经过平滑处理后的数据集,此时单击界面右下角的“Close”按钮或界面右上角 的关闭图标,关闭导入数据界面即可。 2.3.4.数据筛选 所谓的数据筛选就是按照用户设定的条件去除不需要的数据。单击数据拟合主界面上的“Exclude” 按钮,弹出数据筛选界面如图 1-14 所示。

图 1-14 数据筛选界面

图 1-15 手动选点去除数据界面

数据筛选界面用来创建数据筛选规则。在界面左下方的“Exclude Sections”区域,有 4 个下拉菜单和 4 个编辑框,用来设定去除数据的条件,可以根据这些条件创建一个数据筛选规则,并且所创建的数据筛

2010 年 11 月

Matlab 数据拟合——回归分析

·24·

选规则可用于已定义的数据集。用户可在“Exclusion rule name”后面的编辑框中输入数据筛选规则的名 称,然后单击界面下方的“Create exclusion rule”按钮,就可创建一个数据筛选规则。 “Select data set:”后面的下拉菜单用来选择已定义的数据集。一旦选择了某个数据集,在“Check to exclude point:”下方的区域将出现该数据集的预览,可以查看数据集中的每一组观测,并且每一组观测数 据的前面都有一个复选框,通过勾选某些复选框,可以去除相应的观测数据。 选择已定义的 HeadAndAge 数据集,单击“Exclude graphically”按钮,出现手动选点去除数据界 面,如图 1-15 所示。这个界面的左侧区域显示了观测数据的散点图,默认情况下散点为蓝色的圆圈,表 示全部观测均为选入状态,没有被去除的数据。界面右侧有 1 个下拉菜单,2 个单选按钮,3 个普通按 钮。其中下拉菜单用来选择因变量。2 个单选按钮分别为“Excludes Therm”和“Includes Therm”。若选 择“Excludes Therm”,表示当前操作为去除数据,此时在图中散点上单击鼠标左键,选中的散点变为红 色的叉号,相应的观测被去除。若选择“Includes Therm”,则表示当前操作为选入数据,此时在图中散 点上单击鼠标左键,选中的散点变为蓝色的圆圈,相应的观测被选入。3 个普通按钮分别为“Exclude All”、“Include All”和“Close”。若单击“Exclude All”按钮,则表示去除所有观测数据,此时图中散 点全部变为红色的叉号。若单击“Include All”按钮,则表示选入所有观测数据,此时图中散点全部变为 蓝色的圆圈。若单击“Close”按钮,则关闭手动选点去除数据界面。 2.3.5.数据拟合 导入头围和年龄的数据之后,曲线拟合主界面的坐标系里出现了相应的散点图。单击曲线拟合主界 面上的“Fitting”按钮,在弹出的数据拟合界面上单击“New fit”按钮,将创建一个新的拟合,出现如图 1-16 所示界面。

图 1-16 数据拟合界面

界面上“Fit name:”后面的编辑框用来输入拟合的名称;“Data set:”后面的下拉菜单用来选择数据 集;“Exclusion rule:”后面的下拉菜单用来选择数据筛选规则;“Center and scale X data”前面有一个复 选框,勾选该复选框,将对自变量观测数据进行中心化和归一化处理;“Type of fit:”后面的下拉菜单用 来选择拟合类型。可选的拟合类型如表 1-7 所示。

MATLAB 技术论坛电子期刊 第 3 期
表 1-7 可选的拟合类型列表 拟合类型 Custom Equations 说 明 基本模型表达式 广义线性方程: a sin( x ? π ) + c 一般方程: ae
? bx

·25·

自定义函数类型,可在 基本模型表达式基础上 修改 指数函数

+c

Exponential

Fourier

傅立叶级数

aebx aebx + ce dx a0 + a1 cos( xw) + b1 sin( xw)

a0 + a1 cos( xw) + b1 sin( xw) +
Gaussian 高斯函数

+ a8 cos(8 xw) + b8 sin(8 xw)

a1e? (( x ?b1 ) / c1 )

2

a1e? (( x ?b1 ) / c1 ) +
2

+ a8e ? (( x ?b8 ) / c8 )

2

Interpolant Polynomial Power Rational Smoothing Spline Sum of Sin Functions

插值 多项式函数 幂函数 有理分式函数 光滑样条 正弦函数之和

linear、nearest neighbor、cubic spline、shape-preserving 1 至 9 次多项式

axb , axb + c
分子为常数、1 至 5 次多项式,分母为 1 至 5 次多项式 无

a1 sin(b1 x + c1 ) a1 sin(b1 x + c1 ) + + a8 sin(b8 x + c8 )

Weibull

威布尔函数

abx b ?1e ? ax

b

当选中某种拟合类型后,“Type of fit:”下面的空白区域(模型预览区)将出现该拟合类型下的基本 模型表达式,可以通过鼠标选择模型表达式。特别地,当选择自定义函数类型时,模型预览区的右边将 出现 4 个按钮。单击“New”按钮,可在弹出的界面中输入自定义模型表达式。对于自定义函数类型下的 一般方程,还可以设定参数初值和取值范围。单击“Edit”按钮,可对自定义模型表达式进行编辑。在模 型预览区中选中某个模型,然后单击“Copy Edit”按钮,可复制该模型表达式并进行编辑;单击 “Delete”按钮,将删除选中的模型表达式。 单击模型预览区下面的“Fit options”按钮,在弹出的界面中可以设定拟合算法的控制参数,当然也 可以不用设定,直接使用参数的默认值。勾选模型预览区下面的“Immediate apply”复选框或单击 “Apply”按钮,将启动数据拟合程序,数据拟合结果在“Results”下方的结果预览区显示。主要显示模 型表达式、参数估计值与估计值的 95%置信区间和模型的拟合优度。其中模型的拟合优度包括残差平方 和(SSE)、判定系数(R-square)、调整的判定系数(Adjusted R-square)和均方根误差(RMSE)。 在结果预览区的下方有一个拟合列表,如果用户创建了多个拟合,将在拟合列表中显示所有拟合相 关结果。通过拟合列表下方的“Table options”按钮可以控制列表中显示的内容,默认情况下,列表中只 显示拟合的名称(Fit name)、数据集(Data set)、方程名称(Equation name)、残差平方和、判定系 数。在创建多个拟合的情况下,可通过拟合列表对比拟合效果的优劣,可以用残差平方和、调整的判定 系数和均方根误差作为对比的依据。残差平方和越小,均方根误差也越小,调整的判定系数则越大,可 认为拟合的效果越好。对于效果不好的拟合,可以通过拟合列表下方的“Delete fit”按钮删除该拟合。在 拟合列表中选中某个拟合,然后单击拟合列表下方的“Save to workspace”按钮,可以将拟合的相关结果 导入 MATLAB 工作空间。 从图 1-16 可以看出,对于前面给出的 1281 组头围和年龄的观测数据,至少可以用 5 种函数(见 2.1 节)进行拟合,得到的非线性回归方程分别为

2010 年 11 月

Matlab 数据拟合——回归分析
? y = 52.43e x + 0.7906 , x + 0.6644 ? y= , 0.01907 x + 0.01779 ? y = 45.22( x + 0.05)0.05907 , 50.3 ? y= , 1 + 32.96e ? ( x + 4.634) ? y = 45.18 + 2.859 ln( x + 0.05).
? 0.2676

·26·

其中负指数函数和双曲线函数的拟合效果较好。拟合效果如图 1-17 所示。通过图 1-17 所示界面的 File 菜单的“Print to Figure”选项可以把拟合效果图导入单独的图形窗口,而“Generate M-file”选项则 可以生成数据拟合的 M 文件,文件中包含了数据拟合相关操作的 M 代码。

图 1-17 头围和年龄数据的拟合效果图

2.3.6.绘图控制 单击曲线拟合主界面上的“Plotting”按钮,弹出绘图控制界面,如图 1-18 所示。

图 1-18 绘图控制界面

从图 1-18 可以看到,绘图控制界面被分成了左右两个区域,左侧“Plot data sets”区域显示数据集列 表,右侧“Plot fits”区域显示拟合列表。每个数据集和拟合的前面都有一个复选框,通过勾选某些复选

MATLAB 技术论坛电子期刊 第 3 期

·27·

框,可以选择用于绘图的数据集和拟合曲线。若勾选了界面下方的“Clear associated fits when clearing data sets”复选框,则当取消某个数据集的选中状态后,与该数据集相关的拟合也都将处于取消选中状态,图 1-17 中相应的散点和拟合曲线也会被删除。 2.3.7.后处理 拟合结束后,单击曲线拟合主界面上的“Analysis”按钮,将弹出后处理界面,如图 1-19 所示。利用 该界面可进行插值、外推、差分和积分等后处理,这里所谓的插值是在自变量的样本数据的取值范围内 选定一些新的坐标点,计算这些点处因变量的估计值、置信区间和预测区间等。而外推是指新设定的自 变量的取值超出了自变量的样本数据的取值范围。

图 1-19 后处理界面

后处理界面上“Fit to analyze:”后面的下拉菜单用来选择拟合,“Analyze at Xi = ”后面的文本编辑 框用来输入插值或外推点处横坐标的取值,可以输入标量或向量。 勾选“Evaluate fit at Xi”前面的复选框,然后在“Prediction or confidence bounds:”下面的 3 个单选 按钮中选择其一。若选择“None”,则不计算置信区间或预测区间;若选择“For function”,则计算预 测值的置信区间;若选择“For new observation”,则计算观测值的预测区间。通过“Level”后的文本编 辑框设置置信水平,默认情况下文本编辑框中的数字为 95,表示区间估计的置信水平为 95%. “1st derivative at Xi”和“2nd derivative at Xi”复选框分别表示计算回归函数在指定自变量取值处的 一阶导数和二阶导数。 若勾选“Integrate to Xi”前面的复选框,则计算回归函数的积分值,积分上限为自变量的当前指定 值,积分下限通过其下面的 2 个单选按钮来确定。若选择“Start from min(Xi)”,则表示积分下限为自变 量的样本观测值的最小值;若选择“Start from”,则需在“Start from”后面的文本编辑框中输入积分下 限。 以上操作完成后,单击界面下方的“Apply”按钮,将按照用户的选择进行相关的计算,然后在界面 的右侧区域以表格形式显示计算结果。可通过单击“Save to workspace”按钮将计算结果导入 MATLAB 工作空间。 界面左下角的“Plot results”复选框表示根据计算结果绘制图形,“Plot data set:HeadAndAge”则表 示绘制所选数据集 HeadAndAge 的散点图。

3 多重回归分析
在有氧锻炼中,人的耗氧能力 y (ml/(min ? kg)) 是衡量身体状况的重要指标,它可能与以下因素有 关:年龄 x1 (岁),体重 x2 (kg),1500 米跑所用的时间 x3 (min),静止时心速 x4 (次/min),跑步 后心速 x5 (次/min). 对 24 名 40 至 57 岁的志愿者进行了测试,结果如表 1-8 所示。表 1-8 中的数据保存

2010 年 11 月

Matlab 数据拟合——回归分析

·28·

在文件 examp01_03.xls 中,试根据这些数据建立耗氧能力 y 与诸因素之间的回归模型。
表 1-8 人体耗氧能力测试相关数据 序 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 号

y
44.6 45.3 54.3 59.6 49.9 44.8 45.7 49.1 39.4 60.1 50.5 37.4 44.8 47.2 51.9 49.2 40.9 46.7 46.8 50.4 39.4 46.1 45.4 54.7

x1
44 40 44 42 38 47 40 43 44 38 44 45 45 47 54 49 51 51 48 47 57 54 52 50

x2
89.5 75.1 85.8 68.2 89 77.5 76 81.2 81.4 81.9 73 87.7 66.5 79.2 83.1 81.4 69.6 77.9 91.6 73.4 73.4 79.4 76.3 70.9

x3
6.82 6.04 5.19 4.9 5.53 6.98 7.17 6.51 7.85 5.18 6.08 8.42 6.67 6.36 6.2 5.37 6.57 6 6.15 6.05 7.58 6.7 5.78 5.35

x4
62 62 45 40 55 58 70 64 63 48 45 56 51 47 50 44 57 48 48 67 58 62 48 48

x5
178 185 156 166 178 176 176 162 174 170 168 186 176 162 166 180 168 162 162 168 174 156 164 146

3.1 调用自编 reglm 函数作多重回归分析
3.1.1.调用 reglm 函数作多重线性回归分析 对于多重回归分析,由于自变量较多,理论回归方程的选择是比较困难的。这里先尝试作 5 重线性 回归分析,假设 y 关于 x1 , x2 , , x5 的理论回归方程为

? y = b0 + b1 x1 + b2 x2 + b3 x3 + b4 x4 + b5 x5 .
下面调用 1.3.2 节中的自编函数 reglm,求未知参数 b0 , b1 , 著性检验的方差分析表。
% 从 Excel 文件 examp01_03.xls 中读取数值型数据 >> xydata = xlsread('examp01_03.xls'); >> y = xydata(:, 2); % 提取矩阵 xydata 的第 2 列数据,即耗氧能力数据 y >> X = xydata(:, 3:7); % 提取矩阵 xydata 的第 3 至 7 列数据,即自变量观测值矩阵 X % 调用 reglm 函数作 5 重线性回归,显示回归分析的方差分析表和参数估计表 >> reglm(y,X)

(1.5)

, b5 的估计值,在屏幕上显示参数估计表和显

MATLAB 技术论坛电子期刊 第 3 期

·29·

------------------------------------方差分析表-----------------------------------方差来源 自由度 平方和 均方 F值 p值 回归 5.0000 625.3110 125.0622 16.0069 0.0000 残差 18.0000 140.6340 7.8130 总计 23.0000 765.9450 均方根误差(Root MSE) 因变量均值(Dependent Mean) 2.7952 47.6750 判定系数(R-Square) 调整的判定系数(Adj R-Sq) 0.8164 0.7654

-------------------------------参数估计------------------------------变量 估计值 标准误 t值 p值 常数项 121.1655 17.4064 6.9610 0.0000 X1 -0.3471 0.1435 -2.4185 0.0264 X2 -0.0167 0.0874 -0.1914 0.8504 X3 -4.2903 1.0268 -4.1784 0.0006 X4 -0.0399 0.0942 -0.4236 0.6769 X5 -0.1587 0.0788 -2.0122 0.0594

根据上面的参数估计表可以写出经验回归方程如下:

? y = 121.1655 ? 0.3471x1 ? 0.0167 x2 ? 4.2903 x3 ? 0.0399 x4 ? 0.1587 x5 . (1.6) 由于方差分析表中的 p 值小于 0.0001,说明整个回归方程(1.6)式是极显著的,但是并不能说明方程中的 每一项都是显著的。参数估计表中列出了对(1.5)式中常数项和各线性项进行的 t 检验的 p 值,可以看出, x2 , x4 和 x5 所对应的 p 值均大于 0.05,说明在显著性水平 0.05 下,回归方程中的线性项 x2 , x4 和 x5 都是
不显著的,其中 x2 最不显著,其次是 x4 ,然后是 x5 . 下面将(1.5)式中最不显著的线性项 x2 , x4 去掉,重新写出理论回归方程

? y = b0 + b1 x1 + b3 x3 + b5 x5 ,
然后重新调用自编函数 reglm 作 3 重线性回归分析,相应的命令和结果如下:
>> X135 = X(:,[1 3 5]); % 提取矩阵 X 的第 1,3,5 列 >> varnames = {'X1','X3','X5'}; % 定义变量标签 >> reglm(y, X135, [], varnames); % 调用 reglm 函数重新作 3 重线性回归 ------------------------------------方差分析表-----------------------------------方差来源 自由度 平方和 均方 F值 p值 回归 3.0000 623.7206 207.9069 29.2364 0.0000 残差 20.0000 142.2244 7.1112 总计 23.0000 765.9450 均方根误差(Root MSE) 因变量均值(Dependent Mean) 2.6667 47.6750 判定系数(R-Square) 调整的判定系数(Adj R-Sq) 0.8143 0.7865

-------------------------------参数估计------------------------------变量 估计值 标准误 t值 p值 常数项 118.0135 14.3399 8.2297 0.0000 X1 -0.3254 0.1288 -2.5274 0.0200 X3 -4.5694 0.7741 -5.9026 0.0000 X5 -0.1561 0.0750 -2.0809 0.0505

从以上结果可以看出,剔除线性项 x2 , x4 后的经验回归方程为

? y = 118.0135 ? 0.3254 x1 ? 4.5694 x3 ? 0.1561x5 .

(1.7)

该方程也是极显著的,并且常数项和线性项 x1 , x3 也是显著的,只是 x5 所对应的 p 值略大于 0.05,即在显

2010 年 11 月

Matlab 数据拟合——回归分析

·30·

著性水平 0.05 下, x5 接近于显著。结合 x5 所代表的实际意义,还是应将 x5 保留在回归方程中。 3.1.2.调用 reglm 函数作二次回归分析 虽然(1.7)式中已经剔除了最不显著的线性项 x2 , x4 ,并且整个方程是极显著的,但是不能认为(1.7)式 就是最好的回归方程,还应尝试作非线性回归分析,例如二次回归。假设 y 关于 x1 , x2 , 方程为

, x5 的理论回归
(1.8)

? y = b0 + ∑ bi xi + ∑
i =1

5

4

i =1 j =i +1

∑ b x x + ∑b x .
ij i j i =1 2 ii i

5

5

这是一个完全二次方程(包括常数项、线性项、交叉乘积项和平方项)。可调用 reglm 函数求方程(1.8)式 中未知参数 b0 , b1 , , b5 , b12 , b13 , , b45 , b11 , , b55 的估计值,并进行显著性检验。
>> reglm(y,X,'quadratic') % 调用 reglm 函数作完全二次回归

------------------------------------方差分析表-----------------------------------方差来源 自由度 平方和 均方 F值 p值 回归 20.0000 765.0129 38.2506 123.1051 0.0010 残差 3.0000 0.9321 0.3107 总计 23.0000 765.9450 均方根误差(Root MSE) 因变量均值(Dependent Mean) 0.5574 47.6750 判定系数(R-Square) 调整的判定系数(Adj R-Sq) 0.9988 0.9907

-----------------------------------参数估计----------------------------------变量 估计值 标准误 t值 p值 常数项 1804.0553 176.6722 10.2113 0.0020 X1 -26.7679 3.3174 -8.0690 0.0040 X2 -16.4220 1.4725 -11.1526 0.0015 X3 -7.2417 17.3281 -0.4179 0.7041 X4 1.7071 1.5284 1.1169 0.3454 X5 -5.5878 1.2082 -4.6248 0.0190 X1*X2 0.1885 0.0148 12.7017 0.0011 X1*X3 0.2381 0.2163 1.1006 0.3514 X1*X4 0.0168 0.0158 1.0649 0.3650 X1*X5 0.0350 0.0115 3.0359 0.0560 X2*X3 -0.5616 0.0879 -6.3874 0.0078 X2*X4 0.0031 0.0058 0.5294 0.6332 X2*X5 0.0679 0.0064 10.6823 0.0018 X3*X4 -0.0656 0.0713 -0.9206 0.4251 X3*X5 0.1751 0.0639 2.7408 0.0713 X4*X5 -0.0017 0.0056 -0.2968 0.7860 X1*X1 0.0340 0.0223 1.5240 0.2249 X2*X2 -0.0024 0.0031 -0.7907 0.4868 X3*X3 0.6882 0.6357 1.0826 0.3583 X4*X4 -0.0164 0.0048 -3.4342 0.0414 X5*X5 -0.0077 0.0027 -2.8577 0.0647

由于方差分析表中的 p 值为 0.001,说明在显著性水平 0.05 下, y 关于 x1 , x2 ,

, x5 的完全二次回归

方程是显著的。由参数估计值列表可写出经验回归方程,这里从略。从参数估计值列表中的显著性检验
2 的 p 值可以看出,常数项、 x1 、 x2 、 x5 、 x1 x2 、 x2 x3 、 x2 x5 和 x4 所对应的 p 值均小于 0.05,说明回归

方程中的这些项是显著的。读者可以尝试去除不显著项,重新作二次回归分析。

MATLAB 技术论坛电子期刊 第 3 期 3.1.3.拟合效果图

·31·

上面调用 reglm 函数作了 5 重线性回归拟合、3 重线性回归拟合和完全二次拟合,得出了 3 个经验回 归方程。从误差标准差 σ 的估计值(即均方根误差)可以看出 3 种拟合的准确性,均方根误差越小,说 明残差越小,拟合也就越准确。当然也可以从拟合效果图上直观地看出拟合的准确性,下面作出 3 种拟 合的拟合效果图,相关 MATLAB 命令如下:
>> s1 = reglm(y,X); % 调用 reglm 函数作 5 重线性回归拟合

>> X135 = X(:,[1 3 5]); % 提取矩阵 X 的第 1,3,5 列 >> varnames = {'X1','X3','X5'}; % 定义变量标签 >> s2 = reglm(y, X135, [], varnames); % 调用 reglm 函数作 3 重线性回归拟合 >> s3 = reglm(y,X,'quadratic'); % 调用 reglm 函数作完全二次回归拟合

% 绘制拟合效果图 >> figure; % 新建图形窗口 >> plot(y,'ko'); % 绘制因变量 y 与观测序号的散点 >> hold on >> plot(s1.yhat,':','linewidth',2); % 绘制 5 重线性回归的拟合效果图,蓝色虚线 >> plot(s2.yhat,'r-.','linewidth',2); % 绘制 3 重线性回归的拟合效果图,红色点划线 >> plot(s3.yhat,'k','linewidth',2); % 绘制完全二次回归的拟合效果图,黑色实线 >> legend('y 的原始散点','5 重线性回归拟合','3 重线性回归拟合','完全二次回归拟合'); % 标注框 >> xlabel('y 的观测序号'); % 为 X 轴加标签 >> ylabel('y'); % 为 Y 轴加标签
65 y的原始散点 5重线性回归拟合 3重线性回归拟合 完全二次回归拟合

60

55

50

y

45

40

35

0

5

10 15 y的观测序号

20

25

图 1-20 拟合效果图

以上命令作出的拟合效果图如图 1-20 所示,横坐标是因变量的观测序号,纵坐标是因变量的取值。 单纯从拟合的准确性来看,完全二次回归拟合的拟合效果较好,5 重和 3 重线性回归拟合的拟合效果差不 多,相对都比较差。

3.2 调用 stepwise 函数作逐步回归
3.2.1.stepwise 函数的用法 MATLAB 统计工具箱中提供了 stepwise 函数,用来作交互式逐步回归分析,寻找最优回归方程,其 调用方式如下:

2010 年 11 月

Matlab 数据拟合——回归分析

·32·

stepwise stepwise(X,y) stepwise(X,y,inmodel,penter,premove)

(1)stepwise 打开一个交互式图形用户界面(GUI),对 MATLAB 自带的数据文件 hald.mat 中的变量 heat 和 ingredients 进行交互式逐步回归分析,其中 heat 是因变量观测值向量,ingredients 是设计矩阵。 (2)stepwise(X,y) 打开交互式图形用户界面,对用户指定的数据进行交互式逐步回归分析。输入参数 X 为 n × p 的设计 矩阵,y 为因变量的观测值向量,是 n ×1 的列向量。stepwise 函数把 y 或 X 中的不确定数据 NaN 作为缺 失数据而忽略它们。

注意:stepwise 函数自动在 X 中加入一列 1 元素(即模型中自动包含常数项),不需要用户自己
添加。 (3)stepwise(X,y,inmodel,penter,premove) 用 inmodel 参数指定初始模型中所包含的项,inmodel 可以是一个长度与 X 的列数相等的逻辑向量, 也可以是一个下标向量(其元素取值介于 1 和 X 的列数之间,表示列序号)。用 penter 参数指定变量进 入模型的最大显著性水平(默认值为 0.05),显著性检验的 p 值小于 penter 的变量才有可能被引入模 型。用 premove 参数指定从模型中剔除变量的最小显著性水平(默认值为 penter 和 0.1 的最大值),显著 性检验的 p 值大于 premove 的变量有可能被剔除出模型。penter 参数的值必须小于或等于 premove 参数的 值。 3.2.2.调用 stepwise 函数作逐步回归 对于表 1-8 中的人体耗氧能力测试相关数据,考虑线性模型,调用 stepwise 函数作交互式逐步回归分 析,命令如下:
% 从 Excel 文件 examp01_03.xls 中读取数值型数据 >> xydata = xlsread('examp01_03.xls'); >> y = xydata(:, 2); % 提取矩阵 xydata 的第 2 列数据,即耗氧能力数据 y >> X = xydata(:, 3:7); % 提取矩阵 xydata 的第 3 至 7 列数据,即自变量观测值矩阵 X >> inmodel = 1:5; % 初始模型中除了常数项,还包含 x1 至 x5 等线性项 >> stepwise(X,y,inmodel); % 交互式逐步回归分析

图 1-21 逐步回归分析的交互式图形用户界面

运行上面的命令,打开一个图形用户界面,如图 1-21 所示。界面的左上方,“Coefficients with Error Bars”的下面有一些水平条,它们显示了所有项(包括已在模型中的项和有待引进模型中的项)的系数的 估计值的 90%(彩色)和 95%(灰色)置信区间,水平线条中点处圆点的横坐标为系数估计值。红色线 条表示该项不在模型中,蓝色线条表示该项已在模型中。用鼠标单击某个水平线条,可以选入该项或剔 除该项。界面的中上方显示了系数估计值、 t 检验统计量的观测值和显著性检验的 p 值,根据 p 值的大

MATLAB 技术论坛电子期刊 第 3 期

·33·

小可以确定将要进入模型的项和将要剔除的项。 p 值最大并且大于 premove 参数值的项最先被剔除出模 型, p 值最小并且小于 premove 参数值的项最先被引进模型。界面的中部显示了与整个模型相关的一些 统计量,包括常数项的估计值“Intercept”、判定系数“R-square”、对模型进行显著性检验的 F 检验统 计量的观测值“F”、均方根误差“RMSE”、调整的判定系数“Adj R-sq”和显著性检验的 p 值。随着 模型中项的调整,这些统计量也会有相应的变化。界面的下方显示了模型的历史“Model History”,即随 着模型的变化均方根误差 RMSE 的变化轨迹,蓝色圆点的横坐标表示逐步回归的步次,纵坐标表示该步 的均方根误差,用鼠标单击某个蓝色圆点,可以查看该步的结果。 单击界面上“Stepwise”菜单的“Scale Inputs”选项可以对输入数据进行标准化变换。通过下面两种 操作可以实现交互式逐步回归分析: (1)单击界面右上方的“Next Step”按钮,单击一次就执行一步回归,引进最显著的项或剔除最不 显著的项,按钮上方的文字信息清楚地显示了将要引进的项或将要剔除的项。一直单击“Next Step”按 钮,直到均方根误差 RMSE 达到一个局部最小值,此时按钮上方的文字信息变为 "Move no terms",表示 没有可以引进的项,也不再有可以剔除的项,逐步回归结束。以上操作的替代操作是单击界面右上方的 “All Steps”按钮,一次完成所有步的回归。 (2)用鼠标单击界面左上方的水平线条,同样可以进行“引进”和“剔除”的操作。若单击红色线 条,线条颜色将变为蓝色,说明该项被引进模型;若单击蓝色线条,线条颜色将变为红色,说明该项被 剔除出模型。 逐步回归结束后,单击界面左边的“Export”按钮,可以将相关结果导入 MATLAB 工作空间。

注意:stepwise 函数得到的最终结果与初始模型有关,初始模型不同,最终结果也可能不同。
对于本例,初始模型中包含了常数项和线性项 x1 , x2 , x3 , x4 , x5 ,通过单击“Next Step”按钮,依次剔 除 x2 和 x4 ,最终得到的经验回归方程如(1.7)式所示。读者还可以尝试改变设计矩阵,在模型中加入一些 非线性项,然后调用 stepwise 函数作交互式逐步回归,寻找更优的回归方程。


相关文章:
微波实验报告 smith原图 matlab
(包括百度知道、 百度文库、 以及 matlab 技术论坛电子期刊等、 学习了 guide 的重要语句的使用, ) 包括 handles 语句、str2double 语句、set 语句、delete 语句...
科学文献检索课程作业-研究生版本
北京交大图书馆中文期刊论文电子资源有: 维普中国科技...总结和参考文献) 前言:结构损伤识别技术近些年取得...Matlab 各种信息,可以交流讨论 [5]中华钢结构论坛 ...
信号的调制与解调(完整版)
机械电子工程系 专业班级:09 应用电子技术 学生姓名...(自己选,一般取常见的信号), 利用 MATLAB 分析幅度...图 13 封装的 DPSK 调制与解调电路图 第一路信号...
信息检索
WEB全文信息检索技术 5页 免费如要投诉违规内容,请到...我在“万方学术期刊全文数据库”中进行了本次课题的...“MATLAB/SIMULINK 在电力电子学基础课程中的应用”...
期刊论文
各类期刊论文库 11页 免费如要投诉违规内容,请到百度...基于 MATLAB/Simulink 的 平面连杆机器人的动力学 ...云南省第一届科学技术论坛集萃 2003-01-10 发表 ...
数据库实验报告
MATLAB技术论坛电子期刊... 41页 免费数​据​库​实​验​报​告...市场职 111 1120020836 何婷老师 成绩 第七、八周周四 上午一、二节 启动 SQL...
MATLAB在电力电子中的应用
≤π。 第五章 MATLAB 在电力系统中应用学习心得这次的 MATLAB 电力电子技术设计给了我理论转化为实际的机会,作为一名工科学校 的学生,我们学习的理论知识不仅仅是...
文献检索报告
中文期刊论文 [数据库名称]万方数据库 [检索时间]:....Matlab 在数字图像处理中的应用[J]:荆门职业技术...当时的电子计算机已经发展到一定水平, 人们开始利用...
信号的采样与恢复
(2.2.3) 2 设计内容及步骤 2.1 用 MATLAB ...[外文期刊] 2000 [2]张厥盛;张会宁;刑静 锁相环...《核电子学与探测技术》 2007 第 5 期 - 维普...
课题综合检索报告
二、文献检索范围及检索策略: 1、超星数字图书馆 2、CNKI 期刊全文数据库 2、...在 MATLAB 中应用 PSB 进行电气传动系统混合仿 真 [A]. 第七届中国电力电子...
更多相关标签:
电子期刊 | 电子期刊制作 | 电子测量技术 期刊 | 电子技术期刊 | 电子科技期刊 | 电子期刊评职称能用吗 | 电子器件期刊 | 现代电子技术 期刊 |