数学建模社区-数学中国

标题: Matlab数学建模学习报告(一) [打印本页]

作者: 杨利霞    时间: 2019-4-10 15:18
标题: Matlab数学建模学习报告(一)
Matlab数学建模学习报告(一)8 i/ J9 s, G2 Y" X
一、学习目标。

(1)了解Matlab与数学建模竞赛的关系。

(2)掌握Matlab数学建模的第一个小实例—评估股票价值与风险。

(3)掌握Matlab数学建模的回归算法。


) ^0 Q6 F1 g, S% v, c二、实例演练。
' Y+ H. Y1 X& C: q) y; o* b  L( r6 Z* e1 {  a% Y5 E
   1、谈谈你对Matlab与数学建模竞赛的了解。
% o! X- Y3 M0 U( d3 n: a
. W2 J5 l  Z" T. ]0 L        Matlab在数学建模中使用广泛:MATLAB 是公认的最优秀的数学模型求解工具,在数学建模竞赛中超过 95% 的参赛队使用 MATLAB 作为求解工具,在国家奖队伍中,MATLAB 的使用率几乎 100%。虽然比较知名的数模软件不只 MATLAB。/ b) {! |1 M/ ^: q, H' {% y

0 D. O8 }' K& O        人们喜欢使用Matlab去数学建模的原因:
0 X) r8 o7 P4 d; o+ b
( q4 C, x: a: d, W3 h1 i2 l. c(1)MATLAB 的数学函数全,包含人类社会的绝大多数数学知识。
: ?3 B) a! h) s/ a% H6 l/ i) s. a) M, q, M7 ]
(2)MATLAB 足够灵活,可以按照问题的需要,自主开发程序,解决问题。5 H' G- O1 I) f% [  q/ H- T

$ ~! v6 j7 F* ?( y4 M: `- H/ G(3)MATLAB易上手,本身很简单,不存在壁垒。掌握正确的 MATLAB 使用方法和实用的小技巧,在半小时内就可以很快地变成 MATLAB 高手了。
; H, R1 F: V5 P, _6 r' Z* e1 n: M4 m$ W8 W* s* A6 s
        正确且高效的 MATLAB 编程理念就是以问题为中心的主动编程。我们传统学习编程的方法是学习变量类型、语法结构、算法以及编程的其他知识,因为学习时候是没有目标的,也不知道学的知识什么时候能用到,收效甚微。而以问题为中心的主动编程,则是先找到问题的解决步骤,然后在 MATLAB 中一步一步地去实现。在每步实现的过程中,遇到问题,查找知识(互联网时代查询知识还是很容易的),定位方法,再根据方法,查询 MATLAB 中的对应函数,学习函数用法,回到程序,解决问题。在这个过程中,知识的获取都是为了解决问题的,也就是说每次学习的目标都是非常明确的,学完之后的应用就会强化对知识的理解和掌握,这样即学即用的学习方式是效率最高,也是最有效的方式。最重要的是,这种主动的编程方式会让学习者体验到学习的成就感的乐趣,有成就感,自然就强化对编程的自信了。这种内心的自信和强大在建模中会发挥意想不到的力量,所为信念的力量。# `3 X" z+ U$ P# F! [6 x1 q
3 `3 ~4 j; \$ V8 O: Z9 |" `
         数学建模竞赛中的 MATLAB 水平要求:
& o* \0 \. j, ]: X& y( x% q( |% h$ g+ R7 K
要想在全国大学生数学建模竞赛中拿到国奖, MATLAB 技能是必备的。 具体的技能水平应达到:
- d9 D2 E" d( L6 y6 L0 ^
% n% v$ j) h, k4 @4 k) D1)了解 MATLAB 的基本用法,包括几个常用的命令,如何获取帮助,脚本结构,程序的分节与注释,矩阵的基本操作,快捷绘图方式;
! l1 i- c- z& f$ d& F7 ~# Q% B3 q+ _3 y4 F7 g- V" e" R) O+ }
2)熟悉 MATLAB 的程序结构,编程模式,能自由地创建和引用函数(包括匿名函数);3 T6 j" s2 E$ v% b
9 ~+ |4 A  y) O) w
3)熟悉常见模型的求解算法和套路,包括连续模型,规划模型,数据建模类的模型;
+ b3 ~; n, i" Y: K0 N8 e  P( O" A# Y5 i$ }; K- }8 m. u2 h3 f3 e
4)能够用 MALTAB 程序将机理建模的过程模拟出来,就是能够建立和求解没有套路的数学模型。
7 E3 r( y: O8 k& P( t" Z( C2 \$ x
/ @' }# R' O1 V: L: t  @4 `5 `要想达到如上要求, 不能按照传统的学习方式一步一步地学习, 而要结合上述提到的学习理念制定科学的训练计划。5 ~: V0 ]+ }6 v! D' z

0 c! d# `+ R! P7 w  2、已知股票的交易数据:日期、开盘价、最高价、最低价、收盘价、成交量和换手率,试用某种方法来评价这只股票的价值和风险。如何用MATLAB去求解该问题?(交易数据:点击此处获取数据)' {/ \; ]/ A9 Q' z$ }
8 F; ]/ x% n  c) u- R
解题步骤:1 }! n* @$ D( L$ G) r  z2 G5 d
) ^) o, Y) Y+ w' j
第一阶段:从外部读取数据
2 @! K( l2 o( P% J
8 b! U* e7 z6 j  B) F, i" RStep1.1:把数据文件sz000004.xls拖曳进‘当前文件夹区’,选中数据文件sz000004.xls,右键,将弹出右键列表,很快可发现有个“导入数据”菜单,如图 1 所示。
  E0 e- [. j, ]0 Y. i3 I6 P3 f  G9 f1 P3 v
7 A6 E: l4 R* a; @

* \( A' f( ?6 d* ]" v' _4 o                                                                  图1. 启动导入数据引擎示意图) M0 J6 [3 `: R+ C+ u

% V$ V& S% l5 K+ T: G# r+ q, JStep1.2:单击“导入数据”这个按钮,则很快发现起到一个导入数据引擎,如图 4 所示。9 w# m- k( ~$ b3 l  }- n% V. Y

4 R: d0 h9 `2 ?1 h* [0 W( R; P7 ^8 A
5 S+ n2 D1 `* s* E& c' \" s) J+ j9 @
                                                                    图2. 导入数据界面
% X/ I( T* L$ }- O. x. w
( j" x" H1 d% j% L$ m$ ^& [$ I, K% zStep1.3:观察图 2,在右上角有个“导入所选内容”按钮,则可直接单击之。马上我们就会发现在 MATLAB 的工作区(当前内存中的变量)就会显示这些导入的数据,并以列向量的方式表示,因为默认的数据类型就是“列向量”,当然您可以可以选择其他的数据类型,大家不妨做几个实验,观察一下选择不同的数据类型后会结果会有什么不同。至此,第一步获取数据的工作的完成。8 I& D% k0 p8 ]: [) s; X' y

, Q( [$ ^3 K% Y" R
- {5 ]8 r8 b& i8 e1 ^4 T  K% e! c5 m" I5 E5 `" A0 Y8 r$ F, ?
第二阶段:数据探索和建模5 E3 F; R5 i' v" D, ^! U3 ]- r6 P5 j

3 s5 `3 e: m; ~$ o# [8 `现在重新回到问题,对于该问题,我们的目标是能够评估股票的价值和风险,但现在我们还不知道该如何去评估,MATLAB 是工具,不能代替我们决策用何种方法来评估,但是可以辅助我们得到合适的方法,这就是数据探索部分的工作。下面我们就来尝试如何在 MATLAB 中进行数据的探索和建模。4 m. n, i4 X! h/ T

+ W% b1 L4 _3 E. p$ kStep2.1:查看数据的统计信息,了解我们的数据。具体操作方式是双击工具区(直接双击这三个字),此时会得到所有变量的详细统计信息。通过查看这些基本的统计信息,有助于快速在第一层面认识我们所正在研究的数据。当然,只要大体浏览即可,除非这些统计信息对某个问题都有很重要的意义。数据的统计信息是认识数据的基础,但不够直观,更直观也更容易发现数据规律的方式就是数据可视化,也就是以图的形式呈现数据的信息。下面我们将尝试用 MATLAB 对这些数据进行可视化。! G+ H* H) U. C' q

6 j" z" j( `3 j' o5 C由于变量比较多,所以还有必要对这些变量进行初步的梳理。对于这个问题,我们一般关心收盘价随时间的变化趋势,这样我们就可以初步选定日期(DateNum)和收盘价(Pclose)作为重点研究对象。也就是说下一步,要对这这两个变量进行可视化。5 {+ J7 {. }4 ~  Y

! x2 Z0 t$ M* |/ Z对于一个新手,我们还不知道如何绘图。但不要紧,新版 MATLAB 提供了更强大的绘图功能——“绘图”面板,这里提供了非常丰富的图形原型,如图 3 所示。) S6 {$ D" j6 w
3 K0 m: j. q1 c2 M( W. C1 P! n

6 D9 [" k2 U. K; J
% n' U9 ~( [$ {+ o2 l3 R; |                                                                                 图3 MATLAB绘图面板中的图例
% ^2 y2 {! k! B3 d
6 i; H' ]+ [" U- i要注意,需要在工作区选中变量后绘图面板中的这些图标才会激活。接下来就可以选中一个中意的图标进行绘图,一般都直接先选第一个(plot)看一下效果,然后再浏览整个面板,看看有没有更合适的。下面我们进行绘图操作。* `9 T; R2 [- M# @  B

+ H7 i$ u6 H" a, m5 HStep2.2:选中变量 DataNum 和 Pclose,在绘图面板中单机 plot 图标,马上可以得到这两个变量的可视化结果,如图 4 所示,同时还可以在命令窗口区看到绘制此图的命令:, Z' ~0 j% i; W1 V% o/ ]- g
& s. Y; L7 |  V
>> plot(DateNum,Pclose)
9 X# i5 Q7 {1 k$ p
: G9 w' W" A+ e: y' d
5 R1 s4 _# g  Y6 E, [- M# C/ b9 n% g
9 _& s1 h9 x* K- j+ Q9 i                                                                                       图4 通过 plot 图标绘制的原图
) H: J/ i/ c& S/ h
# c+ d- E% u$ {! S4 y/ q8 y这样我们就知道了,下次再绘制这样的图直接用 plot 命令就可以了。一般情况下,用这种方式绘图的图往往不能满足我们的要求,比如我们希望更改:& k& e! ^; C7 ]- s

7 H' r4 P6 Z* ]$ S: Q) a# r(1)曲线的颜色、线宽、形状;
9 h9 x" W! f2 e6 k' Y1 \, f& `4 k
(2)坐标轴的线宽、坐标,增加坐标轴描述;
0 L6 o0 p" N& G: m7 \+ \7 J/ G  E6 z# E
(3)在同个坐标轴中绘制多条曲线。2 [4 g6 B7 B" r, b

! V* Y4 M5 v7 j7 I此时我们就需要了解更多关于命令 plot 的用法,这时就可以通过 MATLAB 强大的帮助系统来帮助我们实现期望的结果。最直接获取帮助的两个命令是 doc 和 help,对于新手来说,推荐使用 doc,因为 doc 直接打开的是帮助系统中的某个命令的用法说明,不仅全,而且有应用实例,这样就可以“照猫画虎”,直接参考实例,从而将实例快速转化成自己需要的代码。* C2 [( y; f# I3 ?
) X: U1 ^' \* ]
接下来我们就要考虑如何评估股票的价值和风险呢?
9 n  {7 w! K* K# e: k/ v$ j: e# V9 M5 [2 O( |) \
         对于一只好的股票,我们希望股票的增幅越大越好,体现在数学上,就是曲线的斜率越大越好。
$ z! |+ }6 K1 A7 r0 y/ H4 S! i7 w
         对于风险,则可用最大回撤率来描述更合适,什么是最大回撤率?: B- Y- G: l0 j2 M* |3 b

. \) P  y* l* q  Z  r         最大回撤率的公式可以这样表达:% e7 e5 i9 m- N) p. {7 {* ?

  u- `. {2 `" \2 e, \& a' ?1 \7 nD为某一天的净值,i为某一天,j为i后的某一天,Di为第i天的产品净值,Dj则是Di后面某一天的净值) G* q7 a! Y0 \7 \

  S+ t& K. q! Z: ~. o5 Gdrawdown=max(Di-Dj)/Di,drawdown就是最大回撤率。其实就是对每一个净值进行回撤率求值,然后找出最大的。可以使用程序实现。最大回撤率越大,说明该股票的风险越高。所以最大回撤率越小,股票越好。0 {* X' t5 J) U5 \2 d; [5 [, p

) u( }0 B" \% v- H3 n, G           斜率和最大回撤率不妨一个一个来解决。我们先来看如何计算曲线的斜率。对于这个问题,比较简单,由于从数据的可视化结果来看,数据近似成线性,所以不妨用多项式拟合的方法来拟合该改组数据的方程,这样我们就可以得到斜率。( w! X5 T& y  e3 n5 _$ B$ ^5 A4 V

" j; K9 ^& {4 C2 E  b% ]2 }4 lStep2.3:通过polyfit()多项式拟合的命令,并计算股票的价值,具体代码为:+ r) Z: l+ O* f
$ a/ R. b) j. [  c$ n2 ]6 S
>> p = polyfit(DateNum,Pclose,1); % 多项式拟合
/ f+ T1 W, k: f( E7 c& E5 O+ ]+ M1 G
>> value = p(1) % 将斜率赋值给value,作为股票的价值0 @$ ?8 K" V% t! K* \

( {3 t% k6 b; \# G. \value =, @! [) t8 g! O+ W* a
* b% R1 `- ?1 o( O' Q
    0.1212
9 E" x6 F2 F4 G: m! M7 ]" v; ?6 J0 ?) E& v* T7 ?& k% {0 D/ k3 H* s
代码分析:%后面的内容是注释。polyfit()有三个参数,前两个大家都能明白是什么意思,那第三个参数是什么意思呢?它表示多项式的阶数,也就是最高次数。比如:在本例中,第三个参数为1,说明其为一次项,即一次函数。第三个参数为你要拟合的阶数,一阶直线拟合,二阶抛物线拟合,并非阶次越高越好,看拟合情况而定。polyfit()返回阶数为 n 的多项式 p(x) 的系数,p 中的系数按降幂排列。在本例中的P(1)指的是最高项的系数,即斜率。
9 E) U3 ?- p' Y. H( w' ^* }: B; X5 G* U& v# O  }) O
Step2.4:用相似的方法,可以很快得到计算最大回撤的代码:- }* V; j+ c: @9 Z
1 ]: _8 l; s; `, w5 E+ Q
>> MaxDD = maxdrawdown(Pclose); % 计算最大回撤
0 P& r( C) h# `# D7 |4 ?5 o) R
# C# b" p; c+ [- H; \>> risk = MaxDD  % 将最大回撤赋值给risk,作为股票的风险- A+ o. Q5 o7 i: G# |/ `2 _

% ~+ G2 W9 N9 T/ i/ {risk =5 ~4 r# N: p, }+ @$ ~! |5 A) R
, ~, x. E  z& b- Q
    0.11558 k/ }3 n# f8 _  d9 ]
2 _5 m) S( U4 m( A- h
代码分析:最大回撤率当然计算的是每天收盘时的股价。最大回撤率越大,说明该股票的风险越高。所以最大回撤率越小,股票越好。  u% F/ g. H3 o* a3 Y) C
5 X- [9 P2 {1 x& I
到此处,我们已经找到了评估股票价值和风险的方法,并能用 MALTAB 来实现了。但是,我们都是在命令行中实现的,并不能很方便地修改代码。而 MATLAB 最经典的一种用法就是脚本,因为脚本不仅能够完整地呈现整个问题的解决方法,同时更便于维护、完善、执行,优点很多。所以当我们的探索和开发工作比较成熟后,通常都会将这些有用的程序归纳整理起来,形成脚本。现在我们就来看如何快速开发解决该问题的脚本。
& B3 m% j3 p; R; x1 m2 |
8 A2 {! C! G5 j, ~, @! kStep2.5:像 Step1.1 一样,重新选中数据文件,右键并单击“导入数据”菜单,待启动导入数据引擎后,选择“生成脚本”,然后就会得到导入数据的脚本,并保存该脚本。
9 |0 G" Y' R, I3 h( z4 F8 m; f' E5 M3 l+ X% ~/ r; B! P6 H
脚本源代码中有些地方要注意:
$ W  X; C- U: Y; f9 h* u5 g; r( |8 k, G
       %%在matlab代码中的作用是将代码分块,上下两个%%之间的部分作为一块,在运行代码的时候可以分块运行,查看每一块代码的运行情况。常用于调试程序。%%相当于jupyter notebook中的cell。4 p1 D) B9 K8 [8 I' [
2 f3 p( W" X, C! o- ^1 u6 [) }- {
       %后的内容是注释。  J0 S- a) M+ y3 P! \4 g1 P

& w8 x. {, C0 f1 R  D        每句代码后面的分号作用为不在命令窗口显示执行结果。
$ T) a7 U4 e6 V, q
- _: \, r0 K) F% f' {2 s脚本源代码:/ }, x- j2 t0 i6 q, z! M
. o5 W2 x; T6 [, w0 n+ |1 x4 n
%% 预测股票的价值与风险4 z& i' ~/ {- V- _! e

, }( h% `/ n2 X# [3 C%% 导入数据/ ]) f, }/ O: u. s
clc, clear, close all, ], z: H0 e) X5 X
% clc:清除命令窗口的内容,对工作环境中的全部变量无任何影响
$ p  w! _# a9 L) u% clear:清除工作空间的所有变量 1 u, J" F4 `' t. ?/ W
% close all:关闭所有的Figure窗口
( U% a. ~# _9 m! p# H
/ B- t8 E" W& W+ D$ Y: k% 导入数据7 X7 [& Q2 k; \& i* B0 }
[~, ~, raw] = xlsread('sz000004.xlsx', 'Sheet1', 'A2:H7');7 a* d$ @, o, ]5 ^
% [num,txt,raw],~表示省略该部分的返回值
  e) z* D1 x3 R1 b, P7 t' ~  d7 s0 M% xlsread('filename','sheet', 'range'),第二个参数指数据在sheet1还是其他sheet部分,range表示单元格范围
4 o; ^$ s- ?8 X0 o
1 G8 L7 i. q* q; @% 创建输出变量0 Z6 H2 M- c, t' n% P
data = reshape([raw{:}],size(raw));' ~5 x1 r) N4 g( q7 }
% [raw{:}]指raw里的所有数据,size(raw):6 x 8 ,该语句把6x8的cell类型数据转换为6x8 double类型数据
5 i5 N" [4 J; I$ C
; p; T$ U. `3 O# Q: Z* B% 将导入的数组分配列变量名称
+ }4 f: j. O) J5 CDate = data(:, 1); % 第一个参数表示从第一行到最后一行,第二个参数表示第一列
! E8 N8 L3 f6 w4 C- A! f( F& c- @DateNum = data(:, 2);& ?4 D# L& Y6 z
Popen = data(:, 3);
( y8 ^4 _0 M1 e" OPhigh = data(:, 4);: U, V; q, z# x& O* b, ^# n2 }
Plow = data(:, 5);/ k9 v. [# ]9 I
Pclose = data(:, 6);  ( y/ R! ?! \1 j% @  [0 D5 H& X8 \
Volum = data(:, 7); % Volume 表示股票成交量的意思,成交量=成交股数*成交价格 再加权求和$ s, J3 [" @/ D4 l
Turn = data(:, 8); % turn表示股票周转率,股票周转率越高,意味着该股股性越活泼,也就是投资人所谓的热门股
* |# @  \/ g0 j# t. n, }; S8 G$ l7 R! D; R
% 清除临时变量data和raw
* T- o/ u, d: t- f6 O0 g0 Aclearvars data raw;* }" P' ]- s% H9 s8 L

+ Q7 F  x* d5 c5 x%% 数据探索
9 ^' V/ v7 |0 h- k
+ ^) o" ?9 i* A1 P% cfigure % 创建一个新的图像窗口
5 ]6 W0 Q9 A) @* X" A) xplot(DateNum, Pclose, 'k'); % 'k',曲线是黑色的,打印后不失真: u6 Q0 I4 I! ^) V5 q$ g
datetick('x','mm-dd'); % 更改日期显示类型。参数x表示x轴,mm-dd表示月份和日。yyyy-mm-dd,如2018-10-27
* _; s6 d+ G. S! N0 rxlabel('日期') % x轴1 k" l+ l: H& u7 U1 B$ G
ylabel('收盘价') % y轴2 j& E1 V) S" u  x$ e3 O
figure! ]9 J. s) S2 i0 u! M& I# }& P
bar(Pclose) % 作为对照图形/ D' W/ I, |" H7 ]/ `
7 r/ U0 C. X" y# Z
%% 股票价值的评估
# v6 R. [, e( _8 H3 l- h7 V9 J0 k* q/ h/ I1 _6 j9 |' Y8 f
p = polyfit(DateNum, Pclose, 1); % 多项式拟合) u6 V/ \4 S, {1 [) Y0 _
% polyfit()返回阶数为 n 的多项式 p(x) 的系数,p 中的系数按降幂排列
: Y7 X; R2 ?' o9 nP1 = polyval(p,DateNum); % 得到多项式模型的结果2 E3 i2 h& e, f1 H
figure- d1 Q* O+ z# c& E5 `; i
plot(DateNum,P1,DateNum,Pclose,'*g'); % 模型与原始数据的对照, '*g'表示绿色的*0 _8 @0 G$ F9 T! N7 H" I, S6 ^
value = p(1) % 将斜率赋值给value,作为股票的价值。p(1)最高项的次数
' F+ a9 \# V" H& i8 p& Q
/ O, y* b' ]1 a4 d. L- `* M%% 股票风险的评估1 V) o7 R0 L# O# M0 u
MaxDD = maxdrawdown(Pclose); % 计算最大回撤2 E& K0 r) Q! ~5 b% r$ x' a
risk = MaxDD  % 将最大回撤赋值给risk,作为股票的风险
' {& v; m% n' d9 L2 h  3、回归算法演练。  Y: E. o) y) V( h/ e4 O  L
& W& I; m, E+ N- U. G& u
(1)一元线性回归3 O3 c7 `! Q5 k* f; h$ g9 {8 f
9 z$ u, p' `9 i$ O+ T3 b
[ 例1 ] 近 10 年来,某市社会商品零售总额与职工工资总额(单位:亿元)的数据见表1,请建立社会商品零售总额与职工工资总额数据的回归模型。
% A4 f/ h2 e7 [* V. J7 m$ ~3 f* t: p5 i, H) b, D% e* Q
* D, v  t/ X" H( Q! d$ l. y
% N5 }# f  {( d$ I: _
该问题是典型的一元回归问题,但先要确定是线性还是非线性,然后就可以利用对应的回归方法建立他们之间的回归模型了,具体实现的 MATLAB 代码如下:3 s% |3 }0 l+ z9 U
0 Z, q. k6 @- F% o% o# o
(1)输入数据' H. n. \- c7 z. w7 f$ q# k, k& ]

  |: T! t& Q, p4 E' ]/ `/ |%% 输入数据
" k1 e5 m/ }! _, a4 B. m- g& {clc, clear, close all
+ v3 H" |' @" U% ~% 职工工资总额2 @2 y" o" _- P" Z% z* G
x = [23.8,27.6,31.6,32.4,33.7,34.90,43.2,52.8,63.8,73.4];
0 _0 L0 T' m' W1 ^% y2 ]' Z% 商品零售总额# \4 E1 N, o* q) e1 i! N
y = [41.4,51.8,61.7,67.9,68.7,77.5,95.9,137.4,155.0,175.0];0 u& K1 B, b" E/ l% {0 p
(2)采用最小二乘回归$ w* v/ y4 Z6 h. ]% K
& H6 P/ B3 b: A
%% 采用最小二乘法回归, A# a* O2 v" k( I! @
% 作散点图5 B7 n1 e, v; m. G$ c: h' Y
figure, i! r& `! @% M/ f
plot(x,y,'r*') % 散点图,散点为红色
% E& E5 j( P( o& M# _xlabel('x(职工工资总额)','fontsize',12)
6 }+ g8 s9 p& Z' [6 ]8 ]ylabel('y(商品零售总额)','fontsize',12)
3 C% x. k8 A$ Y0 `* sset(gca, 'linewidth',2) % 坐标轴线宽为2$ [( k! i" ^5 J5 Y" l5 J6 Q  d4 S+ W

5 M0 d; p  X0 m: u8 Q0 R% 采用最小二乘法拟合
  @% v3 G5 O- N7 ]9 M2 v. QLxx = sum((x-mean(x)).^2); %在列表运算中,^与.^不同! d; n! q$ x! |  d0 {
Lxy = sum((x-mean(x)).*(y-mean(y)));
/ {4 F1 L# @9 d! k; p7 zb1 = Lxy/Lxx;
6 T* F0 b' Q) h& I7 @/ n2 C% Y: X" kb0 = mean(y) - b1 * mean(x);
6 c+ k/ E3 E& K1 dy1 = b1 * x + b0;
) }5 w9 T  a7 f/ `& Z) ~: j- v5 z- |' }5 W/ F& l, d% K
hold on % hold on是当前轴及图像保持而不被刷新,准备接受此后将绘制的图形,多图共存* o: A5 {) _3 l+ H7 I( ^
plot(x,y1, 'linewidth',2);
: G+ D0 `+ c/ \4 q4 _+ M: ~3 {运行本节程序,会得到如图5所示的回归图形。在用最小二乘回归之前,先绘制了数据的散点图,这样就可以从图形上判断这些数据是否近似成线性关系。当发现它们的确近似在一条线上后,再用线性回归的方法进行回归,这样也更符合我们分析数据的一般思路。" R& R$ ^6 a: Z" a/ m' ^8 T
! X/ H3 K+ k3 Y' l4 Z

8 c, t- C; D+ {3 B- v; B
+ j0 Y% J) k/ O4 v% M* P3 C' a                                                                                                    图57 p( G4 \. l' k6 V9 [5 Y0 I

/ K1 o- p6 Y* c9 l! M1 V$ N% z(3)采用 LinearModel.fit 函数进行线性回归
7 h# p2 l% B$ J6 a0 |  c
' K4 t/ U# ]$ K' K! t1 G%% 采用 LinearModel.fit 函数进行线性回归! e/ q% r* Q; |8 \9 S' H: }% a
m2 = LinearModel.fit(x, y)
4 l4 T  Z  Q! }8 ^+ I运行结果如下:% C: S. D6 d  y' P: |4 i8 w

$ {( Q3 j' q. K0 l9 p5 Y$ Fm2 =/ A) I3 N% D8 V' x0 Q0 j! n4 u

3 h2 i7 K, W: L. Y- d; n% o" G+ RLinear regression model:
1 ]- t3 O0 X7 O% J: K' z3 Q$ E7 O- {  Y7 I) ?6 ^! D9 Z, o2 e
    y ~ 1 + x1. o1 C& h6 r4 M( _. ?. e
Estimated Coefficients:$ P$ v# i$ u9 j7 a1 p* L

0 u6 x: n$ B& Y, D, \               Estimate      SE       tStat       pValue
; b: Q# q8 j9 M$ \3 G) v( k, y( R7 B# F) a3 E* h2 S, g' S5 d
    (Intercept)    -23.549      5.1028    -4.615     0.0017215( Z2 Z. S7 a2 o4 r% T
, U7 T- B3 X* s! [5 `/ E# w
    x1           2.7991     0.11456    24.435    8.4014e-09
2 t/ r- l, F) L3 Q" j" [9 u( G! @# O( e+ E
R-squared: 0.987,  Adjusted R-Squared 0.985
( B% K" u2 H  ~% j
! g5 ]/ Z- y2 B7 w! L% d9 Z5 qF-statistic vs. constant model: 597, p-value = 8.4e-09
. j% ~* l9 n& h+ x
- g0 R( v4 v- E5 }0 M/ ]如下图,我们只需记住-23.594是一次函数的中x的系数,2.7991是一次函数中的常数项即可,其它的不用理会。- S$ g4 R8 z% p, L( X: z6 Y

2 w  d+ i+ n" x! Z4 u! z$ w/ C2 x  Q: }: w

& s. J( W. C- q: P  Y2 V- b4)采用 regress 函数进行回归
( t+ A" n4 L' h1 y! B/ i3 N8 V4 G  i, s
1 @: A0 w1 d6 L' _7 x. I%% 采用 regress 函数进行回归1 t1 E- y: ^8 f: w/ W. S. s
Y = y', j, h  V' B$ r% N5 {# c
X = [ones(size(x,2),1),x']
% L+ n& W. h/ _/ q4 K  ?/ i0 h[b,bint,r,rint,s] = regress(Y,X)
' D; G1 w# H* F( r$ n0 U运行结果如下:
7 n: [4 ?1 D1 f6 m
- ^! \. {4 [/ V! q& }3 T7 }b =
2 V, t5 t  b2 j- s3 ]$ C1 j+ [5 d: d) w0 O$ r! ^
  -23.5493+ P3 K" p& m+ y& T5 T1 K( c
3 }5 v, e' Q* I6 k9 l+ i
    2.7991
9 M$ K5 i% h& ~+ j- Y2 ~5 C/ f
$ h$ E* A0 V, {% d0 c7 }. ?, S我们只需记住-23.594是一次函数的中x的系数,2.7991是一次函数中的常数项即可,其它的不用理会。
% i0 A$ {" t/ p# Z+ s, b1 T
5 d# l& l7 t9 b5 [# }6 a( P(2)一元非线性回归
0 S9 l. N$ L9 v$ c9 [7 F6 M( n" g- U  C' i# `  f+ v
[ 例2 ] 为了解百货商店销售额 x 与流通率(这是反映商业活动的一个质量指标,指每元商品流转额所分摊的流通费用)y 之间的关系,收集了九个商店的有关数据(见表2)。请建立它们关系的数学模型。0 u! P2 M2 S# H7 m% s* D
' U0 M" t; j0 q

! n$ x' \( f1 }0 i5 `+ @6 ^% a8 f6 M" V2 Z6 b$ P6 V' q4 g
7 q, c$ B) ?* P5 a* u8 R/ t$ o

  P% A1 v5 `* `! @" P9 g        为了得到 x 与 y 之间的关系,先绘制出它们之间的散点图,如图 2 所示的“雪花”点图。由该图可以判断它们之间的关系近似为对数关系或指数关系,为此可以利用这两种函数形式进行非线性拟合,具体实现步骤及每个步骤的结果如下:9 e6 M* g! u3 L+ i1 c
0 U4 z1 [5 A% }; o- i) A- Y
(1)输入数据
! _  ?5 n6 A5 F+ |1 o5 S4 b0 M+ L& D9 p* \
%% 输入数据
. F4 T& T4 x( i9 n) uclc, clear all, close all
; @2 }% ]( m6 Q' R5 F3 }. ~x = [1.5, 4.5, 7.5,10.5,13.5,16.5,19.5,22.5,25.5];
* j$ Q% a7 o/ |9 n1 iy = [7.0,4.8,3.6,3.1,2.7,2.5,2.4,2.3,2.2];. @# b5 D  F0 o
plot(x, y, '*', 'linewidth', 1) % 这里的linewidth指的是散点大小
4 L) k3 _9 ]* T5 b+ K; N! m5 Aset(gca,'linewidth',2) % 设置坐标轴的线宽为2
* k( c* K0 j  W# r* p& |; nxlabel('销售额x/万元','fontsize',12)
; N6 i6 g; h3 }. Lylabel('流通率y/%','fontsize',12)
( E5 e$ O, y: K  V+ P7 }* w(2)对数形式非线性回归) n$ ~* N) ^, W+ M
- N5 a" e, w; L% }
%% 对数形式非线性回归
! N2 q( w2 @; n& s$ U  m7 ~m1 = @(b,x) b(1) + b(2)*log(x);
8 A5 k2 R/ @. T$ D9 b; t+ L& [; lnonlinfit1 = fitnlm(x,y,m1,[0.01;0.01])
) k/ u$ K- S: r1 L/ o/ d- Db = nonlinfit1.Coefficients.Estimate;. {7 u( G3 s* \/ _6 ]3 s5 S5 _" o2 ?
Y1 = b(1,1) + b(2,1)*log(x);) Y! M- @6 m5 ~! p7 M( M
hold on
: D! O) Q! `' {1 G, jplot(x, Y1, '--k', 'linewidth',2)
4 g1 C- W" q* z1 e  V& v' z运行结果如下:
/ ?4 y' ~& n. ]  Y! s1 \1 s0 Y
( `- a4 Y/ A% c3 {7 B" c5 Q2 }nonlinfit1 =& S* G+ ?- k' e: a

8 ?& F+ |5 Z4 y9 dNonlinear regression model:
  u4 \$ R# i( H+ L6 U. k1 v0 u; p7 v- K+ N0 A2 o2 G2 L  L3 M
    y ~ b1 + b2*log(x)) Q* O! H% ]7 p0 ~! a5 s4 @$ p6 D
  F9 B5 F, P9 B, X: ^% s0 P- B/ l
Estimated Coefficients:" h/ T* a  a$ D/ S5 w
$ z( _) ~+ \# D: d; S0 R! {
          Estimate      SE        tStat       pValue - M, C- i6 y3 {. X  ]

1 G  @$ p# U+ _. g$ s: V    b1    7.3979      0.26667     27.742    2.0303e-088 z" N* ^$ v3 n

# ~% k# w5 w4 ^% X! n+ g    b2    -1.713      0.10724    -15.974    9.1465e-07
- [' f" d' b6 \& y1 I7 `3 o7 d2 u( h8 H
R-Squared: 0.973,  Adjusted R-Squared 0.969
2 E/ F6 A3 ~& U& E: [4 G+ D1 r$ M/ y6 J$ ]5 r
F-statistic vs. constant model: 255, p-value = 9.15e-07( t9 E3 {+ ]0 Y! Y
# f  i2 ?& C5 a( k0 [( O
(3)指数形式非线性回归
) F& O/ b% @6 ]8 ?
& y! w& l- v; U%% 指数形式非线性回归
* D) L7 e! [6 r! S* U& U7 _m2 = 'y ~ b1*x^b2';+ J9 K; Y$ D( J) T
nonlinfit2 = fitnlm(x,y,m2, [1;1])
: O, u& Q5 ]* }2 j( ob1 = nonlinfit2.Coefficients.Estimate(1,1);5 m$ z* U: X2 I1 ?$ h& P
b2 = nonlinfit2.Coefficients.Estimate(2,1): G) N* E, _( X. ~3 O
Y2 = b1*x.^b2;4 i- Q1 U; A! o# k( n3 K
hold on;
) d, M5 G( b& x: C* H! x& Zplot(x,Y2,'r','linewidth',2)" ~" C9 S9 J+ k. X
legend('原始数据','a+b*lnx','a*x^b') % 图例0 ?$ y2 M' v9 R! X# y2 m
运行结果如下:) h( w9 x6 k/ j% f. P! u* o
0 ^+ g5 q/ t/ l# Z5 [8 U
nonlinfit2 =
6 Z$ D/ k4 |4 |$ [$ w
+ u" V( j* V( i/ X/ g; CNonlinear regression model:
+ _/ T& y( R( D+ P4 Z9 G$ f
; o" N. w" a. l  u3 {    y ~ b1*x^b2# b+ m6 l9 P8 K- j
6 m! W5 X8 s: ~3 @6 r
Estimated Coefficients:' ?3 Y9 ?9 F5 D# G; O

7 M+ v- n% J1 S' l          Estimate       SE        tStat       pValue 1 c+ `, }5 S: e2 D" R8 G8 ]/ @* G

# D' L7 v9 W2 j6 b    b1      8.4112     0.19176     43.862    8.3606e-10
1 b+ [0 X& h5 y. p- t
# @0 u  O7 i% a' |0 O1 w; ]+ d    b2    -0.41893    0.012382    -33.834    5.1061e-09
" H7 e+ @9 [# d3 z1 i# b: m  u2 M: S
R-Squared: 0.993,  Adjusted R-Squared 0.992
; ^9 F8 y/ @1 Q* q, I) K4 }% V/ t! ~+ T; z2 Z9 g' U- |
F-statistic vs. zero model: 3.05e+03, p-value = 5.1e-11* l0 V9 M) k$ [

- y$ |% d! P; g9 p在该案例中,选择两种函数形式进行非线性回归,从回归结果来看,对数形式的决定系数为 0.973 ,而指数形式的为 0.993 ,优于前者,所以可以认为指数形式的函数形式更符合 y 与 x 之间的关系,这样就可以确定他们之间的函数关系形式了。
1 S& ]& L' B; `/ q9 G% g* v
1 x# S' T" z% {, l2.多元回归8 m4 M8 C* o2 s- v

5 }2 H) f2 G& [8 ]& e1.多元线性回归) V! B1 G" U4 Y: w
% q, Q7 B2 v8 U8 N; R5 \6 [" G, g  j( q
[ 例3 ] 某科学基金会希望估计从事某研究的学者的年薪 Y 与他们的研究成果(论文、著作等)的质量指标 X1、从事研究工作的时间 X2、能成功获得资助的指标 X3 之间的关系,为此按一定的实验设计方法调查了 24 位研究学者,得到如表3 所示的数据( i 为学者序号),试建立 Y 与 X1 , X2 , X3 之间关系的数学模型,并得出有关结论和作统计分析。6 V" X) v( s4 K" ~. a

+ n7 G8 e9 i& t: `" m* i9 U
+ d# s% M- C7 P% Q4 x& \& L" A; }4 b, k5 L
该问题是典型的多元回归问题,但能否应用多元线性回归,最好先通过数据可视化判断他们之间的变化趋势,如果近似满足线性关系,则可以执行利用多元线性回归方法对该问题进行回归。具体步骤如下:
3 w9 X. P+ j7 V6 D* L+ n2 `2 `' @6 `8 Q
. J. I! m) S& v! O  `9 K(1)作出因变量 Y 与各自变量的样本散点图1 O; @  G% r1 t/ `. \

3 e2 H' j# ?, x9 C' p2 {作散点图的目的主要是观察因变量 Y 与各自变量间是否有比较好的线性关系,以便选择恰当的数学模型形式。图3 分别为年薪 Y 与成果质量指标 X1、研究工作时间 X2、获得资助的指标 X3 之间的散点图。从图中可以看出这些点大致分布在一条直线旁边,因此,有比较好的线性关系,可以采用线性回归。绘制图3的代码如下:8 M" V& y: j+ C; |% C! d
, ?+ c. @1 ~" X) b: ^! z  ]
%% 作出因变量Y与各自变量的样本散点图
( M  @, n( `$ k( x3 k% x1,x2,x3,Y的数据- f" o7 C' r2 E" ^! I/ G
x1=[3.5 5.3 5.1 5.8 4.2 6.0 6.8 5.5 3.1 7.2 4.5 4.9 8.0 6.5 6.5 3.7 6.2 7.0 4.0 4.5 5.9 5.6 4.8 3.9];
! x# G8 A% |' Z/ ^% R9 \* g$ {x2=[9 20 18 33 31 13 25 30 5 47 25 11 23 35 39 21 7 40 35 23 33 27 34 15];
* ^% H  ~& Z0 R% m3 [8 L3 L3 Jx3=[6.1 6.4 7.4 6.7 7.5 5.9 6.0 4.0 5.8 8.3 5.0 6.4 7.6 7.0 5.0 4.0 5.5 7.0 6.0 3.5 4.9 4.3 8.0 5.0];
( Y" \" o; l( e  V4 MY=[33.2 40.3 38.7 46.8 41.4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.6 34.2 48.0 38.0 35.9 40.4 36.8 45.2 35.1];
- o8 Z& J5 c/ V$ B. F% 绘图,三幅图横向并排. }% ~* O+ t4 h3 _! M8 a% g  n
subplot(1,3,1),plot(x1,Y,'g*')0 o0 C) T& G- w. l+ e' A% a* K
subplot(1,3,2),plot(x2,Y,'k+')9 S+ n% X& X1 r) x1 z
subplot(1,3,3),plot(x3,Y,'ro')
$ X& J! u) J6 ^+ C4 D' b$ T绘制的图形如下:9 ~5 {! P' b7 y, ?1 P

' g' Z% p  ]: k, G5 h, [! P  k# n/ c3 X% X9 `

/ K( T& }/ j8 h4 ~(2)进行多元线性回归4 T, }4 b, j" X6 x

8 `4 I# l6 H+ s& i" v- J这里可以直接使用 regress 函数执行多元线性回归,注意以下代码模板,以后碰到多元线性问题直接套用代码,具体代码如下:$ |) q! [0 O/ K
* X8 g$ D0 O1 H& t4 W# {  ~
%% 进行多元线性回归9 e: s2 |0 n! G7 s$ T+ [( j
n = 24; m = 3; % 每个变量均有24个数据,共有3个变量
, ?( t" S1 X: }0 y. r2 ~X = [ones(n,1),x1',x2',x3'];
6 x. F( ?3 N% \/ R1 ?1 _+ H[b,bint,r,rint,s]=regress(Y',X,0.05) % 0.05为预定显著水平,判断因变量y与自变量之间是否具有显著的线性相关关系需要用到。2 T# c: ?/ D2 M
运行结果如下:  J" }) A, A% ^. E! h/ d& A4 O& \. ^- w

8 C' K7 i$ ?/ v. Ub =+ n: s# H2 l, Q: T0 _0 S: ~* q
7 v" M8 c% g* {) [
   18.01573 k. w7 p0 D9 o! T! v
    1.0817/ W; i$ l+ f; O: L
    0.32125 m7 P% Z: q  H( w( L
    1.2835& ]/ \; q! j% a

7 W3 Q0 w9 C! U0 q6 x+ W! x; E# X! V
bint =
3 S# d, ]# ?( f/ b3 H( Y, S0 j
" z/ F2 t% Z0 K; y7 o" ~; @   13.9052   22.1262
3 r2 Q* f$ t$ Q, \  S9 S' o0 Y: J    0.3900    1.7733$ x  ?' a. \9 Z1 J& \) O$ p& w8 l9 d
    0.2440    0.3984
# f; m5 E' E' F; \    0.6691    1.8979
! w7 i: M  G, [: k0 F+ K
; Y& m6 a& K$ z3 j  s4 t3 D
4 j9 m& S/ C# q6 K/ Q" Nr =
' s( q- _1 f9 [& ]+ B
5 Q2 t; F6 U" f) F4 N    0.67819 t' E; B+ P; A3 {
    1.91295 V# \- _* o5 [
   -0.1119
2 o4 z+ }" u# L! g( i$ p0 V" r: u    3.3114
  t' _' ]* S1 s- v, i) X   -0.74247 r# n+ N! ?0 r3 X; C) {
    1.2459: w. Z1 _5 [/ b
   -2.1022
+ C! C; f7 \) @( M, y" Z    1.9650: j! s9 q0 e9 x' Q4 ^
   -0.3193
# X& n" r4 |0 v$ H5 S# h" H) e) P    1.3466. I# T" s7 Q+ }, K& X; p
    0.8691& w$ f4 U' U; |9 a' J- G  B$ |
   -3.2637
3 T" M% T" p! Q* t/ H8 }   -0.5115
6 C2 X9 ?  X0 \5 }   -1.1733
3 k, [* u- ~. F: |* r   -1.4910
, {: {# `4 f+ k5 Q; B  M, l: [   -0.2972( j. s) n$ {/ ?2 d
    0.1702
( b, D- ]9 d) u1 e6 }; C5 h    0.5799
' O. m; k* G* W1 H9 q   -3.28568 ^( E# e: G- g2 u( h
    1.1368
; e$ x! z; r( x1 A: C   -0.8864
9 b+ ^; N' Z" _, q, j$ z' ~   -1.4646# E# H: o3 ~4 r. e7 |
    0.8032
4 B8 ^! F3 A! t: N4 [5 e+ i( Q    1.6301, w9 u' {+ z6 n5 ~  z5 `
0 _. T2 b! F2 z
; g; z6 u( @7 C) b' b3 Q! r& m
rint =* h8 X. G, Z1 s: D' O- S! f

& V, ]: l( u, }1 S4 ]7 e. E9 h8 q   -2.7017    4.0580
6 x- r0 R' w- r8 S: ?1 j. \7 p   -1.6203    5.4461
# ]  E( }( }, T7 N% z0 `( H   -3.6190    3.39518 G: |% f- K$ [( Y7 D1 a, C8 H# D
    0.0498    6.5729- i9 s$ M# h3 B7 ^& a2 T# l( M
   -4.0560    2.5712, n4 S( p' e$ y
   -2.1800    4.6717
! ]& @; m+ `7 g   -5.4947    1.2902
. z; r0 E2 r5 X   -1.3231    5.2531
' x" F8 ^- a0 d. o9 c   -3.5894    2.9507- I5 }5 ?( c- z$ l' T4 d7 M
   -1.7678    4.4609
" C% i& Q0 u' l1 W2 k   -2.7146    4.4529
' @, s; f" e9 Q2 J8 ?$ a   -6.4090   -0.11839 D8 |5 ~+ O! x7 i& _" ]% b* o
   -3.6088    2.5859# b; i- ]) h% a8 p  d8 ]  f
   -4.7040    2.3575
2 V( V& h" u+ K2 b% B   -4.8249    1.8429
+ L) G% P7 W3 B. Y% b/ X" S   -3.7129    3.11855 k1 Z1 q% N5 F% |% d0 a# N- _5 w
   -3.0504    3.3907
3 _! Q+ m- l& a2 p/ i" X' N6 D   -2.8855    4.0453$ N! U+ v9 k$ M' k1 u- n8 W$ P& R
   -6.2644   -0.3067
% S2 Q" d$ \, B4 C   -2.1893    4.46302 D- F" e1 t# r% [; i
   -4.4002    2.6273
4 M9 N( }! K7 q* O  Z, B$ j   -4.8991    1.9699* u% ?5 X; f2 }% o
   -2.4872    4.0937# b$ J0 v! k' [4 @$ y* \
   -1.8351    5.0954
* l' H9 o1 `5 g  d  t
/ Z% ~  B. z( X% O$ u
' \% e! I$ T( C/ v9 js =/ n+ n  q9 K9 T) b1 c! S$ a

: c1 h# e3 I/ X7 r. N; f- _    0.9106   67.9195    0.0000    3.07197 Y% M3 v+ p* V/ d5 }" G$ D
看到如此长的运行结果,我们不要害怕,因为里面很多数据是没用的,我们只需提取有用的数据。( X1 i& y7 g3 R+ q: {1 e

! U3 P) g. U& l% g在运行结果中,很多数据我们不需理会,我们真正需要用到的数据如下:
% I7 [9 O; z( u% N! G( c) [8 \6 I; S7 ^2 f3 o2 g
b =( H8 d. \0 d" ~0 `

/ ]" q+ @3 z$ i0 C   18.0157
) \# m+ o2 Z( R5 l$ ?4 \9 J$ G    1.0817
" X9 u6 O, l6 x+ D2 t    0.3212
7 s! B2 d3 V8 F# p7 [    1.2835! [6 f: h( H9 u. g4 {
3 {; C: B6 R7 E+ i5 G% |! A/ \% h6 o
s =7 M+ l8 {3 c& M1 O5 G; m2 G! m0 M

1 K8 f  w1 x7 A& d    0.9106   67.9195    0.0000    3.0719
9 Q; \2 K* I8 t- H% }% D0 T  t回归系数 b = (β0,β1,β2,β3) = (18.0157, 1.0817, 0.3212, 1.2835),回归系数的置信区间,以及统计变量 stats(它包含四个检验统计量:相关系数的平方R^2,假设检验统计量 F,与 F 对应的概率 p,s^2 的值)。观察表4的数据,会发现它来源于运行结果中的b和s:
# G1 G; R  q" r* }8 J: d9 ~9 P- a8 H6 w' n
3 w' Z5 n( E- [- L
* ~* X# z% T# u# c
根据β0,β1,β2,β3,我们初步得出回归方程为:
3 k5 t4 P0 u  h9 X* A8 E% b1 Z2 q" S! H. M( K7 l

9 p# E! O/ r3 I0 V. n$ b" @9 P5 U8 f
如何判断该回归方程是否符合该模型呢?有以下3种方法:7 M  }: B3 \. N2 r/ l: A& M
- G' x- j( Q% P$ I
1)相关系数 R 的评价:本例 R 的绝对值为 0.9542 ,表明线性相关性较强。
: R8 [7 r1 j( c; A1 f' c! U- H  a# g6 X% I
2)F 检验法:当 F > F1-α(m,n-m-1) ,即认为因变量 y 与自变量 x1,x2,...,xm 之间有显著的线性相关关系;否则认为因变量 y 与自变量 x1,x2,...,xm 之间线性相关关系不显著。本例 F=67.919 > F1-0.05( 3,20 ) = 3.10。% o  n% b1 h2 N/ _8 o; S2 i/ w

& m+ M6 W% y+ T0 |3)p 值检验:若 p < α(α 为预定显著水平),则说明因变量 y 与自变量 x1,x2,...,xm之间显著地有线性相关关系。本例输出结果,p<0.0001,显然满足 p<α=0.05。/ x( L2 a  Q" ]3 Q. k

$ Q  \$ ~9 \( s& r以上三种统计推断方法推断的结果是一致的,说明因变量 y 与自变量之间显著地有线性相关关系,所得线性回归模型可用。s^2 当然越小越好,这主要在模型改进时作为参考。; {4 W! q8 U/ x- d  g1 O
* F- q5 J; N- l1 x$ v' q4 U
3. 逐步回归. q. }! D6 g7 I& p6 S3 V, B

: K* _  \+ ~" g+ s# Y- J7 X/ p+ o[ 例4 ] (Hald,1960)Hald 数据是关于水泥生产的数据。某种水泥在凝固时放出的热量 Y(单位:卡/克)与水泥中 4 种化学成品所占的百分比有关:
; P: X4 ^( M3 _$ n; m) E/ K6 M5 V* [/ I, [7 A$ X- N

4 z4 I' X! s- l' f1 ]5 v7 P/ V5 P8 {4 U) Y/ T
在生产中测得 12 组数据,见表5,试建立 Y 关于这些因子的“最优”回归方程。5 D  w1 f) c& _9 C, Y

0 L7 l" F7 m# \" G9 H$ k! w
: x* ~9 \, y* {0 s, n
/ G) E: R% I" [3 ^对于例 4 中的问题,可以使用多元线性回归、多元多项式回归,但也可以考虑使用逐步回归。从逐步回归的原理来看,逐步回归是以上两种回归方法的结合,可以自动使得方程的因子设置最合理。对于该问题,逐步回归的代码如下:
) ~+ e( z7 e, P0 H( L2 D2 W5 i) [6 I
%% 逐步回归- {9 L! E6 h9 i& ^7 p
X=[7,26,6,60;1,29,15,52;11,56,8,20;11,31,8,47;7,52,6,33;11,55,9,22;3,71,17,6;1,31,22,44;2,54,18,22;21,47,4,26;1,40,23,34;11,66,9,12];   %自变量数据6 S( C3 ?" U+ O' K3 S, i' s9 j8 j
Y=[78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3];  %因变量数据
% k/ W- q8 {; z; Y1 {  |" lstepwise(X,Y,[1,2,3,4],0.05,0.10)% in=[1,2,3,4]表示X1、X2、X3、X4均保留在模型中
8 l1 a# _1 i/ Q) z% j) p程序执行后得到下列逐步回归的窗口,如图 4 所示。
3 o& m: L: t, b8 s9 S# P/ }2 ]0 p! l) [2 s
: f5 ?2 Q. G. N) U' G* X% l

; S: u2 X4 W  ]; z; F3 _                                                                                                             图4: T$ @' ]3 ^6 r) Y  L( R$ Q, H. C
* U9 R/ C0 E& X4 W
在图 4 中,用蓝色行显示变量 X1、X2、X3、X4 均保留在模型中,窗口的右侧按钮上方提示:将变量X4剔除回归方程(Move X4 out),单击 Next Step 按钮,即进行下一步运算,将第 4 列数据对应的变量 X4 剔除回归方程。单击 Next Step 按钮后,剔除的变量 X3 所对应的行用红色表示,同时又得到提示:将变量 X3 剔除回归方程(Move X3 out),单击 Next Step 按钮,这样一直重复操作,直到 “Next Step” 按钮变灰,表明逐步回归结束,此时得到的模型即为逐步回归最终的结果。最终结果如下:# y& o- E6 T, N8 `
6 V) y/ @0 ~# Q/ {  ]) a

9 O, j8 ]! H5 i5 z, ~" j
, ^! P% W1 {$ R4. 逻辑回归- ~3 p3 s: [* B+ ]- G; U; u
* @3 B5 S$ c+ x- Y& `4 K0 D" N
[ 例5 ] 企业到金融商业机构贷款,金融商业机构需要对企业进行评估。评估结果为 0 , 1 两种形式,0 表示企业两年后破产,将拒绝贷款,而 1 表示企业 2 年后具备还款能力,可以贷款。在表 6 中,已知前 20 家企业的三项评价指标值和评估结果,试建立模型对其他 5 家企业(企业 21-25)进行评估。
3 e% R  F" k1 B3 [$ @% i
* e& j% Z7 l  W; F& C& i$ {
! l0 H! B! d+ r+ y8 K
' _+ Q: w4 K" m* ~/ ~* v对于该问题,很明显可以用 Logistic 模型来回归,具体求解程序如下:. P/ o4 q4 v" J% ^

# N6 p$ E" p! U- o) z程序中需要用到的数据文件logistic_ex1.xlsx已上传github:https://github.com/xiexupang/mathematical-modeling/tree/master/%E5%9B%9E%E5%BD%92/%E9%80%BB%E8%BE%91%E5%9B%9E%E5%BD%92
6 R6 ]$ B  c. i' S; t( w% _: t' D  q0 L% J
% logistic回归( z6 p' x  R* n4 ?8 R# E: i; n$ b

% o2 W* h5 o8 D9 W9 |%% 导入数据+ ~' \/ M+ a" v
clc,clear,close all
6 ^& q# H4 P% Q- Z% `/ p6 Q( C& M6 NX0 = xlsread('logistic_ex1.xlsx','A2:C21'); % 前20家企业的三项评价指标值,即回归模型的输入
* n5 X' ?2 _! i5 XY0 = xlsread('logistic_ex1.xlsx','D221'); % 前20家企业的评估结果,即回归模型的输出
, H- Q2 J: T/ `' CX1 = xlsread('logistic_ex1.xlsx','A2:C26'); % 预测数据输入
, p3 i8 L  K4 ~, W# y5 s, i9 \' o# D% N5 X
%% 逻辑函数
% Z: z5 s" y* X  H. dGM = fitglm(X0,Y0,'Distribution','binomial');
5 k0 I5 p# J" c- MY1 = predict(GM,X1);
: P) K1 }, N0 x6 K8 u# g, j$ g' s0 N
%% 模型的评估
) e  I" Y' a- a) }N0 = 1:size(Y0,1); % N0 = [1,2,3,4,……,20]
- a$ @0 {' q3 ~& \5 \2 CN1 = 1:size(Y1,1); % N1 = [1,2,3,4,……,25]6 U) m& C) L$ ~7 a) V$ n$ O; x
plot(N0',Y0,'-kd'); % N0'指的是对N0'进行转置,N0'和Y0的形式相同,该行代码绘制的是前20家企业的评估结果  f( s6 A- s2 J1 |' q' k( z; O
% plot()中的参数'-kd'的解析:-代表直线,k代表黑色,d代表菱形符号
/ p) W  O# G; C4 o7 n# xhold on;
3 W! V. E3 _; n, |scatter(N1',Y1,'b'); % N1'指的是对N1'进行转置,N1'和Y1的形式相同
( G8 E: ~: r5 W/ M6 `% pxlabel('企业编号');
  _3 m. q# C- Yylabel('输出值');2 Z' G# t3 h& T9 Q3 k) i
得到的回归结果与原始数据的比较如图5所示。
3 A% u, O( k' V9 h1 @% ~- v3 H  E0 z* l

# ^- e' G3 B# i0 m  @! p( r4 H; d, Z
                                                                   图5" A- _8 x! D3 y, l; Y
3 M/ |1 Q  h8 d! j6 s2 C' \6 |
三、总结与感悟。
9 H7 C* I$ Z' [5 A# I# e/ I- g* m1 h
        总结:通过这次学习,我了解到Matlab在数学建模竞赛中使用广泛;在评估股票价值与风险的小实例中,我掌握了用Matlab去建模的基本方法和步骤;在回归算法的学习过程中,我掌握了一元线性回归、一元非线性回归、多元线性回归、逐步回归、逻辑回归的算法。: T; e+ Y# G( f# b: S8 G& q3 G

5 P& l  V6 v% l2 H' z1 a* i        感悟:正确且高效的 MATLAB 编程理念就是以问题为中心的主动编程。我们传统学习编程的方法是学习变量类型、语法结构、算法以及编程的其他知识,因为学习时候是没有目标的,也不知道学的知识什么时候能用到,收效甚微。而以问题为中心的主动编程,则是先找到问题的解决步骤,然后在 MATLAB 中一步一步地去实现。在每步实现的过程中,遇到问题,查找知识(互联网时代查询知识还是很容易的),定位方法,再根据方法,查询 MATLAB 中的对应函数,学习函数用法,回到程序,解决问题。在这个过程中,知识的获取都是为了解决问题的,也就是说每次学习的目标都是非常明确的,学完之后的应用就会强化对知识的理解和掌握,这样即学即用的学习方式是效率最高,也是最有效的方式。最重要的是,这种主动的编程方式会让学习者体验到学习的成就感的乐趣,有成就感,自然就强化对编程的自信了。这种内心的自信和强大在建模中会发挥意想不到的力量,所为信念的力量。+ X" p& k$ R2 w. Y! o
1 r2 r% ]- N% i' X& [
$ {9 L0 }" B5 J6 Y9 \





欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5