基于MATLAB的多元非线性回归模型.pdf

第 2 9卷第 2期 2 0 0 9年 3月 云南师范大学学报 J o u r n a l o f Y u n n a n N o r ma l Un i v e r s i t y Vo 1 2 9 No 2 Ma r 2 0 09 基 于 MA T L A B的多元非线性 回归模型 董大校(临沧师范高等专科学校，云南 临沧 6 7 7 0 0 0 0)摘要：M A T L AB是源于矩阵运算的一种高度集成的计算机语言。它提供了强大的科学运算、灵活的 程序设计流程、高质量的图形可视化与界面设计、便捷的与其他程序和语言接口的功能。文章充分利用 MA T L A B统计工具箱的优势，通过程序的实现，对多元非线性回归模型的未知参数的估计方法以及对估 计后的模型预报做 出研究，并 以实例验证了该方法的有效性。关键词：MA T L A B；多元非线性 回归；最小二乘法；统计工具箱 中图分类号：T P 3 0 1 文献标识码：A 文章编号：1 0 0 79 7 9 3(2 0 0 9)0 2 0 0 4 5 0 4 1 预备知识 非线性回归最小二乘法拟合的基本原理 。对给定数据(，Y )(i=0，1，m)，在取定的函数类 中，求 P()E ，使误差=P()一Y (i =0，1，m)的平方和最小，即 r =p()一 Y i 最小，从几何意义上讲，就是寻求 与给定点(，Y )(i=0，1，m)的距离平方和为最小的曲线 Y=P()(图1)。函数P()称为拟合 函数或最小二乘解，求拟合 函数P()的方法称为曲线拟合的最小二乘 法 2 M A T L A B非线性 曲线拟合命令介绍 2 1 n l i n fi t 函数 用 n l i n fi t 函数进行非线性最小二乘数据拟合。该函数使用高斯一牛顿算法，调用格式如下：b e t a=n l i n fi t(X，Y，f u n，b e t a 0)用最小二乘法估计非线性函数系数。Y为响应值(因变量)矢量。一般地，为 自变量值组成的设计矩阵，每一行对应与 Y中的一发个值。但是，可以是 f u n参数能接受的 任何数组。f u n参数为一函数，该函数具有下面的形式 y h a t=m y f u n(b e t a，X)其中 b e t a 为系数矢量，x为设计矩阵。f u n为参数返 回一个拟合 Y 值的 y h a t 矢量。b e t a 0为包含系 数初始值的矢量。b e t a，r，J =n l i n fi t(X，Y，f u n，b e t a 0)返 回拟合系数(b e t a)、残差(r)、和雅可 比矩阵 J，这些参数 可以用于 n l i n t o o l 函数，生成预测值的误差估计；或用于 n l p a r c i 函数生成系数的误差估计。2 2 n l i n t o o l 函数 使用 n l i n t o o l 函数可以对数据进行非线性方程拟合并交互图形显示，其调用格式如下：收稿 日期：2 0 0 8 1 0 2 3 作者简介：董大校(1 9 6 5 一)，男，云南省临沧市人，副教授，主要从 事数学教学科研工作 4 6。云南师范大学学报(自然科学版)第 2 9卷 n l i n fi t(x，Y，f u n，b e t a 0)为一预测图，它提供数据(x，Y)的非线性曲线拟合。它用两条红色曲线 来表示预测值的9 5 置信区间。b e t a 0为一矢量，包含参数的初值。f u n参数为一函数，该 函数具有 y h a t =m y fun(b e t a，x)的形式：其中 b e t a 为系数矢量，x为设计矩阵。n l i n t o o l 函数显示一个图形“矢量”，其中的每个 图形对应于输入矩阵 X的每一列。响应变量 y为一 列矢量，与 x的行数相匹配。n l i n t o o l(x，Y，F U N，b e t a 0，a l p h a)用图形表示预测值的 1 0 0(1 一a l p h a)置信 区间。当 x为矩阵 时，将为每一列单独 生成 图形，应变量 Y为一列变量，对应于 x的行。a l p h a()的默认值 为 0 0 5生成 9 5 置信区间。2 3 n l p r e d c i 函数 可利用 n l p r e d c i 函数计算非线性模型预测值 的置信区间，其调用格式如下：y p r e d=n l p r e d c i(F U N，i n p u t s，b e t a，r，J)给定拟合参数(b e t a)、残差(r)和和雅可 比矩阵(J)，返 回预测响应。输入是非线性 函数中独立变量 的数值矩阵。2 4 n l p a r c i 函数 可利用 n l p a r c i 函数计算非线性模型中参数估计值 的置信区间。n l p a r c i(b e t a，r，J)给定拟合参数(b e t a)、残差(r)、和解处的雅可比矩阵(J)，返回非线性最小二 乘参数估计(b e t a)的 9 5 置信区间 c i。n l p a r c i 函数使用 n l i n fi t 函数的输 出作为输入。2 5 r e g r e s s 函数 用 r e g r e s s函数进行多元线性回归。b，b i n t，r，r i n t，s t a t s=r e gre s s(Y，x，a l p h a)给出 b i n t 和 r i n t 的 1 0 0(1一a l p h a)置信区间。2 6 s t e p w i s e函数 用 s t e p w i s e函数进行逐步回归，它使用交互环境进行分析。其调用格式如下：s t e p w i s e(x，Y，i n m o d e l，p e n t e r，p r e m o v e)指定模型的初始状态 和要使用 的置信区间。i n m o d e l 为 长度为 x中列数的逻辑矢量。i n m o d e l 指定包含在初始模型中的 自变量。默认时不包含 x中的列。p e n t e r 指定引入 自变量的最大 P值，默认时为 0 0 5。p r e m o v e 指定剔除 自变量的最小 P值，默认时为 0 1 0。3实例验证 3 1 反应 动 力学模 型 反应动力学中的 H o u g e nWa t s o n模型是非线性模型的一个典型例子。其模型如下：8 L x 2一X 3 8 s r a t e 其 中：，卢 ，为未知参数，。，：和，为三个输入变量。三个输入为氢(h y d r o g e n)、n一戊烷(n p e n t a n e)、异戊烷(i s o p e n t a n e)。上述模型显然为多元非线性模型。文件 r e a c t i o n m a t 包含了反应的仿真 数据。通过计算可知，反应动力学 H o u g e nWa t s o n模型中的 5个未知参数为 l=1 2 5 2 6，卢 2=0 0 6 2 8，卢 3=0 0 4，卢 4=0 1 1 2 4，卢 5=1 1 9 1 4 1 2 5 2 6 ，一 1 1 91 4 最终模型为：r。r _ _ 3 2 财政 收入预 测 财政收入与国民收入、工农业总产值、人 口、就业人 口、固定资产投资等因素有关。采用逐步回归分 析构造预测模 型。以财政 收入作为 因变量 Y，自变量为国 民收入()、工业 总产值()、农业总产值 第 2期 董大校，等：基于 MA T L A B的多元非线性 回归模型 4 7。()、总人口()、就业人 口()、固定资产投资()。其样本数 7,=3 0，自变量数P=6。首先编写非线性函数的 M函数文件 mo d e 1 I n代码。然后在命令窗 口输入由、，、样本数 据构成的矩阵 和财政收人样本数据 y，取 b e t a O=0 5 0 0 0 3-0 6 0 0 O 1 0 0 2 0 3 5 ；用交互式 非线性拟合图形工具来显示，在 M A T L A B命令窗口中输入如下命令：n l i n t o o l(X，Y，i n o d e 1 b e t a O，0 O 1)；得到 n l i n t o o 1 分析窗 口。由此得到财政收入与各因素之间关系的回归模型 1为：Y=0 5 2 4 5 x l一0 0 2 9 3 x 20 6 3 1 x 3+0 O 1 1 3 x 40 0 2 3 1 x 5+0 3 6 5 2 x 6 (1)对该模型进行多元线性回归检验，输人命令：b，b i n t，r，r i n t，s t a t s =r e g r e s s(Y，x)s t a t s=0 9 8 0 5 3 01 0 3 0 5 0 9 51。5 53 0 结果分析：h i n t 为各系数的置信区间。s t a t s 矢量是值分别为相关系数的平方、，值和显著性概率P。相关系数平方值 R。=0 9 8 0 5，说明模型拟合程度较高。显著性概率 P=0，小于 0 0 5，故拒绝零假设，认 为模型 1中至少有一个 自变量的系数不为零，因而从总体上模型 1 是有意义 E l 的。结果 6给出了回归系数的估计值，检查它们的置信区间发现，系数的置信区间包含了零点，表明 回归变量 对变量 l，的影响不是太显著，因此从模型中移 出次变量，对该模 型进行逐步 回归检验，在 M A T L A B命令窗口输入并运行命令：s t e p w i s e(X，Y)除了。，外，其他 自变量的回归系数置信区间都包含零点，保 留，将其他变量移出回 归模型，R s q u a r e=0 9 7 9 6 4 5，是指因变量的9 7 9 6 可由模型确定，F=4 1 7 1 0 4远远超过了F检验的 临界值，P=0远小于显著性水平，这些说明用包含，的线性模型来描述 y是合适的。的回归系 数为 0 3 8 2 8 6 9，的回归系数为 一0 4 9 6 9 7 5，的回归系数为0 4 0 8 6 2 4，常数项由 I n t e r e e p 给出，即 6=1 7 7 0 7，利用逐步回归得到的回归模型 2为：y=0 3 8 2 8 6 9 x l一0 4 9 6 8 7 5 x 3+0 4 0 8 6 2 4 x 6+1 7 7 0 7 (2)由于，对 Y具有较明显的显著性，用更为精确的拟合曲线来拟合表示。，和 Y的关系。做，分别对 Y的散点图，具体步骤在 MA T L A B命令窗 口输入命令：x l=X(：，1)；x 3=X(：，3)；x 6=X(：，6)；s u b p l o t(1，3，1)；p l o t(x l，Y，+，)；s u b p l o t (1，3，2)；p l o t(x 3，Y，+；s u b p o t(1，3，3)；p l o t(x 6，Y，+再做经变换后的散点图，对 比其线性效果较好的，得出如下拟合 曲线：Y=一7 1 3 9 1+0 4 7 3 x l；，=1 6 2 0 4 6 6 5 0 9 x 3+0 0 1 x；y=1 4 1 9 2 1+1 1 5 8 x 6；在 M A T L A B命令窗口输入如下命令：s u b p l o t(1，3，1)；p l o t(x l，一 7 1 3 9 1+0 4 7 3木 x l，+；s u b p l o t(1，3，2)；p l o t(x 3，1 6 2 0 4 6 6 5 0 9术 x 3+0 0 1术 x 3 ，+s u b p l o t(1，3，3)；p l o t(x 6，1 4 1 9 2 1+1 1 5 8半 x 6，+对应的RS q u a r e 分别为 0 9 5 1，0 8 8 7，0 9 4 4。从相关系数的值可以看出，对应的拟合方程的拟合 程度较高。所 以，用拟合后的曲线来代替原来的，得到如下 回归模型 3：Y=0 1 8 1 l x l+3 2 3 4 2 x 3 0 0 0 5 x；+0 4 7 3 2 x 6 5 9 7 4 3 7 2 (3)作为对比，在引入交互项，在 MA T L AB命令窗 口输入下面命令：x l=x(：，1)；x 2=X(：，2)；x 3=x(：，3)；x 7=x 1 木 x 2；x 8=x 1 木x 3；x=x x 7 x 8 ；s t e p w i s e(x，Y)加入交互项后，利用逐步 回归分析，从运算结果可得 RS q u a r e=0 9 8 6 6 5，F=4 6 1 9 2 3，R MS E：2 7 9 2 8 3，P=0 0 0 0 0。各项参数值较模型(2)、(3)更好，所以，建立模型(4)Y：0 8 45 3 61 xI一 0 01 2 0 2 52 x 5+0 0 0 0 0 0 45 x 700 0 05 7 x 8 即：Y=0 8 4 5 3 6 1 x l一0 0 1 2 0 2 5 2 x 5+0 0 o 0 0 o 4 5 1 20 0 0 0 5 7 x l 3 (4)从各统计量的数值看，模型(2)比模型(1)差一点，但模型(2)比模型(1)更为简单，模型(3)较模型 4 8 云南 师范 大学学报(自然科学版)第 2 9卷(1)、(2)虽然是非线性的回归模型，但实际上它的拟合效果 并不好。模型(4)是在加入交互项后建立 的，其系数中有的较小，几乎为零，代人数据进行预测，其效果不是很好。综合上述情形，建立 的线性 回 归模型(2)更能真实的反映各 自变量与因变量之间的关系。4 结语 本文从两方面人手，介绍了基于 M A T L A B的多元非线性 回归模 型参数的求解方法，同时也尝试 了 通过实验建立多元非线性 回归模型的方法，虽然对于 3 2财政收入预测模型并没有表现 的像多元线性 回归模型那样优秀，但是可以提供一种思想与方法，在以后 的学习与实践中来完善这种思想与方法。参考文献：1 袁建问 经济计量学实验 M 北京：科学出版社，2 0 0 2 8 7 8 9 2 苏金明，阮沈勇，王永利 MA T L A B工程数学 M 北京：电子工业出版社，2 0 0 5 5 0 8 3 The M u l t i v a r i a t e No n l i ne a r Re g r e s s i o n M o de l Ba s e d o n M ATLAB DONG Dax i a o (L i n C a n g T e a c h e r s C o l l e g e)Ab s t r ac t：MATLAB i s t h e。dg i n o f t h e I n a t r i x。p e r a t i。n s i n a h i g h l y i nt e g r a t e d c。mp ut e r l a n g ua g e I t p r o v i de d a p owe r f u l s c i e n-t i fi c c。mp u t i n g，n e x i b l e p r o g r a m d e s i g n p r o c e s s，h i g h q u a l i t y g r a p h i c s a n d V i s u a l i n t e r f a c e d e s i g n a n d c。n V e n i e n t P r 0 c e d“。a n d 1 a ng u a g e a n d。t h e r i nt e r f a c e f e a t u r e s I n t h i s p a p e r，we【a k e f u l l a d v a n t a g e o f MATLAB s t a t i s t i c t 0 0 1 b。x t。a c hi e V e t h e mu l t i v a r i a t e n。n l i ne a r r e g r e s s i【)n m。de 1 Es t i ma t e。f t h e un k no wn pa r a me t e r s e x a mp l e s a r e u s e d t。V a l i d a t e t h e me t h。d i s e f f e t i。Ke y wo r(1 s：MATL AB；Mu l t i v a r i a t e No n l i n e a r r e g r e s s i o n；L ea s t s q u a r e s；s t a t i s t i c t o o l b o x (上接第 4 4页)Ent i r e So l ut i o n s o f Fun c t i o na l Equ a t i o n f 6(z)+g (z)+h ()=1 SU Mi nL I Yuh u a f D e p a r t m e n t o f M a t h e m a t i e s Y u n n a n N o r ma l U n i v e r s i t y，Y u n n a n K u n m i n g 6 5 0 0 9 2)A b s t r a c t：T h i s p a p e r s t u d i e s t h e e n t r i e s o l u t i o n s。f t h e f u n c t i。n a l e q u a t i o n f ()+g。()+h 6()=1 I t i s s h o w n t h a t t h e r e i s n o n 0 n c 0 n s t a n t e n t r i e f u n c t i o n s w i t h o r d e r s ma l l e r t h a n o n e s a t i s f y i n g t h e e q u a t i o n z)，g (z)(z)w i t h o r d e r s m a l l e r t h a n o n e s a t i s f y i n g t h e e q u a t i o n f (z)+g ()+h ()=1 Me a n w h i l e a n e w b rie f p r o o f i s g i v e n f o r t h e kn o wn r e s u l t s t h a t i f n7 t h e r e i s n o n o n c o ns t a n t e n t ne f u n c t i o n s s o l u t o n s o t t he f u n c t i o n a l e q u a t i o n f ()+g ()+h ()=1 Ke v wo r ds：F u nc t i o n a l e q u a t i o n；En t i r e f u nc t i o n；Va l u e d i s t r i b u t i o n t h e o r y