数学建模社区-数学中国

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

作者: 杨利霞    时间: 2019-4-10 15:18
标题: Matlab数学建模学习报告(一)
Matlab数学建模学习报告(一)
7 Q( T5 F) M  \3 A" a& ~一、学习目标。

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

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

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


: e2 u0 o- u+ `! h二、实例演练。
# `: _2 [  l; {. A  P# c
* m( _5 e2 p9 p% u& m/ _+ _$ f   1、谈谈你对Matlab与数学建模竞赛的了解。
) i# R3 [! Q+ u
# }" N4 I/ g1 F  s0 d        Matlab在数学建模中使用广泛:MATLAB 是公认的最优秀的数学模型求解工具,在数学建模竞赛中超过 95% 的参赛队使用 MATLAB 作为求解工具,在国家奖队伍中,MATLAB 的使用率几乎 100%。虽然比较知名的数模软件不只 MATLAB。* t7 H) u1 Y9 C
% }9 x& t0 Y- J. t- q
        人们喜欢使用Matlab去数学建模的原因:& c/ ^8 L* {, W7 q# P8 [

! r" a- W0 n) ]5 y; @1 G6 d2 j(1)MATLAB 的数学函数全,包含人类社会的绝大多数数学知识。3 J6 m  I+ S1 n! w

% q$ W6 Z2 ^  _0 E7 X(2)MATLAB 足够灵活,可以按照问题的需要,自主开发程序,解决问题。# s7 t2 J; b5 w. S/ }9 W3 @

$ V/ _4 V( G  W/ g! M9 @(3)MATLAB易上手,本身很简单,不存在壁垒。掌握正确的 MATLAB 使用方法和实用的小技巧,在半小时内就可以很快地变成 MATLAB 高手了。( W  k. ~2 A! n7 G) j: Y! N

/ w1 X8 @+ n0 N' L$ B        正确且高效的 MATLAB 编程理念就是以问题为中心的主动编程。我们传统学习编程的方法是学习变量类型、语法结构、算法以及编程的其他知识,因为学习时候是没有目标的,也不知道学的知识什么时候能用到,收效甚微。而以问题为中心的主动编程,则是先找到问题的解决步骤,然后在 MATLAB 中一步一步地去实现。在每步实现的过程中,遇到问题,查找知识(互联网时代查询知识还是很容易的),定位方法,再根据方法,查询 MATLAB 中的对应函数,学习函数用法,回到程序,解决问题。在这个过程中,知识的获取都是为了解决问题的,也就是说每次学习的目标都是非常明确的,学完之后的应用就会强化对知识的理解和掌握,这样即学即用的学习方式是效率最高,也是最有效的方式。最重要的是,这种主动的编程方式会让学习者体验到学习的成就感的乐趣,有成就感,自然就强化对编程的自信了。这种内心的自信和强大在建模中会发挥意想不到的力量,所为信念的力量。
* ]: Q. l% X. h+ z
/ R, g+ d5 r" Q6 G( b         数学建模竞赛中的 MATLAB 水平要求:+ v8 {9 p: ~5 @$ V& s! h9 Z! h
; n. J. H$ Z$ j  B  C* Q
要想在全国大学生数学建模竞赛中拿到国奖, MATLAB 技能是必备的。 具体的技能水平应达到:9 d, V+ R0 p. w. k- y: Q6 J2 E

' G1 u7 @3 j) Y2 L# \1)了解 MATLAB 的基本用法,包括几个常用的命令,如何获取帮助,脚本结构,程序的分节与注释,矩阵的基本操作,快捷绘图方式;) ^& M0 J  o- t- V  t; J2 s

' S1 ?' v4 P( a7 _2 d* i- r1 Z% ~  i2)熟悉 MATLAB 的程序结构,编程模式,能自由地创建和引用函数(包括匿名函数);
' Y5 _( A+ y2 h
& x, S6 }, E9 i8 N7 [0 x1 H- O3)熟悉常见模型的求解算法和套路,包括连续模型,规划模型,数据建模类的模型;  G9 v& L/ z& P* j4 V6 e' O
0 u0 }- P* d" r  {
4)能够用 MALTAB 程序将机理建模的过程模拟出来,就是能够建立和求解没有套路的数学模型。
  }8 o* v) j3 M& f+ n/ F$ Z# F
% r% [+ C( r' {$ D/ T( {要想达到如上要求, 不能按照传统的学习方式一步一步地学习, 而要结合上述提到的学习理念制定科学的训练计划。
& }$ ^* }8 I  T& j. ~8 S; e1 |! ~5 H, p% V: S" n  h8 R! P
  2、已知股票的交易数据:日期、开盘价、最高价、最低价、收盘价、成交量和换手率,试用某种方法来评价这只股票的价值和风险。如何用MATLAB去求解该问题?(交易数据:点击此处获取数据)
( S* S: p3 j. d$ f3 s* Z2 X( J( S" a  Z; @. O3 h
解题步骤:5 J' c8 Q, B: I7 c
6 |: v+ y5 b/ ]% K: m9 `0 P
第一阶段:从外部读取数据2 R# f& r4 Y; b. s8 H0 y$ ]6 w
& ^6 h( d% A8 F7 o0 w2 E6 e: j
Step1.1:把数据文件sz000004.xls拖曳进‘当前文件夹区’,选中数据文件sz000004.xls,右键,将弹出右键列表,很快可发现有个“导入数据”菜单,如图 1 所示。
4 R2 N* i4 I' C7 v, x
1 ?/ Z/ C3 ^$ Z0 z5 D: T- S
) A' z+ N" ~6 f* m2 q) g" S6 R' C3 M& A  ~% {+ ?/ H7 a
                                                                  图1. 启动导入数据引擎示意图
  s2 f% d' ^: Q% {. s( J% q
; X. |. g6 _% A6 jStep1.2:单击“导入数据”这个按钮,则很快发现起到一个导入数据引擎,如图 4 所示。7 r0 t# p6 L+ Z& ]' S/ _, f

4 I9 l1 o. a; D/ ^% o7 V8 H$ ]3 _, S) ?; y

0 e$ v4 h# v4 w, h* E- s$ g                                                                    图2. 导入数据界面
1 t6 |( h( D  E# |7 L8 r2 u/ _' D; O4 U. n% n2 A
Step1.3:观察图 2,在右上角有个“导入所选内容”按钮,则可直接单击之。马上我们就会发现在 MATLAB 的工作区(当前内存中的变量)就会显示这些导入的数据,并以列向量的方式表示,因为默认的数据类型就是“列向量”,当然您可以可以选择其他的数据类型,大家不妨做几个实验,观察一下选择不同的数据类型后会结果会有什么不同。至此,第一步获取数据的工作的完成。. J# y6 u: n) j1 K9 i& w# l) i2 g

* {+ U* e* Q- o  T1 r, ]  g- c, }8 u
* I  f) j) g8 m% I4 n
$ c. c7 G& P* {* C/ m第二阶段:数据探索和建模* b; `; N5 H( D6 ]+ ]
3 x# @+ A# ]# }( S
现在重新回到问题,对于该问题,我们的目标是能够评估股票的价值和风险,但现在我们还不知道该如何去评估,MATLAB 是工具,不能代替我们决策用何种方法来评估,但是可以辅助我们得到合适的方法,这就是数据探索部分的工作。下面我们就来尝试如何在 MATLAB 中进行数据的探索和建模。
, [0 i3 ]0 R& p; j( T: @- R- n3 f1 j$ o8 H+ x8 b# g
Step2.1:查看数据的统计信息,了解我们的数据。具体操作方式是双击工具区(直接双击这三个字),此时会得到所有变量的详细统计信息。通过查看这些基本的统计信息,有助于快速在第一层面认识我们所正在研究的数据。当然,只要大体浏览即可,除非这些统计信息对某个问题都有很重要的意义。数据的统计信息是认识数据的基础,但不够直观,更直观也更容易发现数据规律的方式就是数据可视化,也就是以图的形式呈现数据的信息。下面我们将尝试用 MATLAB 对这些数据进行可视化。
5 [" @2 K, K& q- L. y* L) U9 [# m' P5 a) p6 I
由于变量比较多,所以还有必要对这些变量进行初步的梳理。对于这个问题,我们一般关心收盘价随时间的变化趋势,这样我们就可以初步选定日期(DateNum)和收盘价(Pclose)作为重点研究对象。也就是说下一步,要对这这两个变量进行可视化。
. B5 M1 X4 u7 c2 F* s9 Y( N$ \7 _. v/ y
对于一个新手,我们还不知道如何绘图。但不要紧,新版 MATLAB 提供了更强大的绘图功能——“绘图”面板,这里提供了非常丰富的图形原型,如图 3 所示。
* I* C! ]5 U+ e' x' Q. p6 Y$ y* [' S+ |! m/ E! ~# d4 E5 }5 u

" O2 r- Y1 i; `% ~. i- J
% P7 s: u  i2 ]4 n7 T5 ~; u                                                                                 图3 MATLAB绘图面板中的图例% n6 Q, X. V& f4 N  W

4 w% M$ o7 O1 ~要注意,需要在工作区选中变量后绘图面板中的这些图标才会激活。接下来就可以选中一个中意的图标进行绘图,一般都直接先选第一个(plot)看一下效果,然后再浏览整个面板,看看有没有更合适的。下面我们进行绘图操作。- v+ Q; Y% Z: M4 V" j' p

" o$ \5 H! A: T6 RStep2.2:选中变量 DataNum 和 Pclose,在绘图面板中单机 plot 图标,马上可以得到这两个变量的可视化结果,如图 4 所示,同时还可以在命令窗口区看到绘制此图的命令:
2 L# C1 H* V, n8 h" w0 D8 H' ?: s- c2 J
>> plot(DateNum,Pclose)
$ T3 M5 i4 O) |. e
- C: M, w3 A' f+ b3 ]! w4 ~  [2 ^% K
( P* p5 K. G# W; o" `
                                                                                       图4 通过 plot 图标绘制的原图
* w- \' L9 G# h/ |" h& V% @6 G# [
  r4 ?$ w7 W( V, ?. e这样我们就知道了,下次再绘制这样的图直接用 plot 命令就可以了。一般情况下,用这种方式绘图的图往往不能满足我们的要求,比如我们希望更改:/ t, N: o: x! N: i3 V8 r

# J2 l. Z- ]) W$ |1 W(1)曲线的颜色、线宽、形状;
4 |8 P" x5 y# Q7 a& Y6 L# z" n! m  k* C5 Y' a+ v
(2)坐标轴的线宽、坐标,增加坐标轴描述;
0 i. ~; ]/ O2 X& Y) \7 D! r! S# O- A
(3)在同个坐标轴中绘制多条曲线。
; o8 p. L" G2 l% _8 m5 q, R* ?3 S
' S8 i* ~+ L! d此时我们就需要了解更多关于命令 plot 的用法,这时就可以通过 MATLAB 强大的帮助系统来帮助我们实现期望的结果。最直接获取帮助的两个命令是 doc 和 help,对于新手来说,推荐使用 doc,因为 doc 直接打开的是帮助系统中的某个命令的用法说明,不仅全,而且有应用实例,这样就可以“照猫画虎”,直接参考实例,从而将实例快速转化成自己需要的代码。3 J  n  d( C  z) j0 {4 e. [: j

  b4 K" D; z/ N( @* s2 }" V: @接下来我们就要考虑如何评估股票的价值和风险呢?
) p3 y3 c& g6 d4 k4 u- y1 g( T2 B' {" W9 o5 E8 P
         对于一只好的股票,我们希望股票的增幅越大越好,体现在数学上,就是曲线的斜率越大越好。
1 ~' T- z  a3 Z. _- `% w% ]3 k2 U! x3 p
         对于风险,则可用最大回撤率来描述更合适,什么是最大回撤率?: S' c( J9 P, H5 d# b% q% y5 z
3 @$ e) K+ e  E4 u$ I
         最大回撤率的公式可以这样表达:) y7 f* F$ k3 M1 r- W" N/ v8 d. z1 w$ ^
+ F* z+ n% p2 ~; a2 }1 Y
D为某一天的净值,i为某一天,j为i后的某一天,Di为第i天的产品净值,Dj则是Di后面某一天的净值; l/ w2 n9 I8 N3 j7 h0 X! E$ N% Z

* I. n& r7 f( ]drawdown=max(Di-Dj)/Di,drawdown就是最大回撤率。其实就是对每一个净值进行回撤率求值,然后找出最大的。可以使用程序实现。最大回撤率越大,说明该股票的风险越高。所以最大回撤率越小,股票越好。1 F& ~& |8 i' c9 C
: G7 J( ~" j: \# W5 D
           斜率和最大回撤率不妨一个一个来解决。我们先来看如何计算曲线的斜率。对于这个问题,比较简单,由于从数据的可视化结果来看,数据近似成线性,所以不妨用多项式拟合的方法来拟合该改组数据的方程,这样我们就可以得到斜率。: y: d6 q! }! U* y
1 t; N0 k/ T* `% H- z+ h8 A% B1 H) {
Step2.3:通过polyfit()多项式拟合的命令,并计算股票的价值,具体代码为:
9 J: X8 R* i- {! M: \2 u( X
7 g" z4 P5 y1 l$ ?7 M>> p = polyfit(DateNum,Pclose,1); % 多项式拟合
6 i) F/ M4 L7 W7 l3 O4 A' M8 w: q  T+ m, d3 {# h9 X! y
>> value = p(1) % 将斜率赋值给value,作为股票的价值
: a% e- y0 i/ H/ I3 s9 d: ?- G9 m/ o/ n0 H) d
value =* \0 z* s* D# v) z& ]7 S& h
( \& B, U2 q' h+ Y7 v* `" i: \
    0.12129 B9 O6 W; C& x! Y

/ Q# O6 {2 l! S代码分析:%后面的内容是注释。polyfit()有三个参数,前两个大家都能明白是什么意思,那第三个参数是什么意思呢?它表示多项式的阶数,也就是最高次数。比如:在本例中,第三个参数为1,说明其为一次项,即一次函数。第三个参数为你要拟合的阶数,一阶直线拟合,二阶抛物线拟合,并非阶次越高越好,看拟合情况而定。polyfit()返回阶数为 n 的多项式 p(x) 的系数,p 中的系数按降幂排列。在本例中的P(1)指的是最高项的系数,即斜率。
' \/ g) u! d8 _6 C+ \4 @
# [+ |3 p7 T. R9 D) j' ]Step2.4:用相似的方法,可以很快得到计算最大回撤的代码:8 n$ i' V& O% q2 a, ]

( Q9 @* k* u: D( z/ _>> MaxDD = maxdrawdown(Pclose); % 计算最大回撤& b9 |8 A0 m. @' }1 ~7 u
- j7 ~' q: {' P% ~8 r
>> risk = MaxDD  % 将最大回撤赋值给risk,作为股票的风险: s/ [6 z' |9 O3 U- C1 D9 d
9 g0 X- t* k5 r* O
risk =
8 s' l$ H" w: U3 ^. w( D( ]9 @! q6 I0 V2 w7 H
    0.1155
/ C6 j+ k9 {; @, u% b1 }
3 z2 i5 O+ g/ ~5 k/ J/ i代码分析:最大回撤率当然计算的是每天收盘时的股价。最大回撤率越大,说明该股票的风险越高。所以最大回撤率越小,股票越好。& F8 R7 D/ ^' G8 F* Z

5 Y. m9 u; z! ]: t: M6 a1 u到此处,我们已经找到了评估股票价值和风险的方法,并能用 MALTAB 来实现了。但是,我们都是在命令行中实现的,并不能很方便地修改代码。而 MATLAB 最经典的一种用法就是脚本,因为脚本不仅能够完整地呈现整个问题的解决方法,同时更便于维护、完善、执行,优点很多。所以当我们的探索和开发工作比较成熟后,通常都会将这些有用的程序归纳整理起来,形成脚本。现在我们就来看如何快速开发解决该问题的脚本。; c! Z+ P# V( A
( H; u% k! y* \( R8 E, s5 F# F
Step2.5:像 Step1.1 一样,重新选中数据文件,右键并单击“导入数据”菜单,待启动导入数据引擎后,选择“生成脚本”,然后就会得到导入数据的脚本,并保存该脚本。
$ n9 Z) g* v0 B1 Z: g3 x! m
; Z0 @4 @. \7 a5 a" f. N( V5 s" h  c脚本源代码中有些地方要注意:$ K5 l7 }+ w$ B

( ^: J- r1 z. `8 h4 f       %%在matlab代码中的作用是将代码分块,上下两个%%之间的部分作为一块,在运行代码的时候可以分块运行,查看每一块代码的运行情况。常用于调试程序。%%相当于jupyter notebook中的cell。
% h" [) x1 z" f8 e8 P6 r% K
3 o; ~. R- K' W" [' \       %后的内容是注释。
+ R6 Y3 [' g2 i" X+ g4 E6 o- h+ G9 s# u; H; e; X( n2 {8 ~
        每句代码后面的分号作用为不在命令窗口显示执行结果。
* a; c; b7 N0 [8 J! n& ^1 V* n& @& D
脚本源代码:
7 |# A9 v. Y; w) _2 _& I
% S' Y% R3 `  L9 I%% 预测股票的价值与风险
  T- x- Y0 A# x, Z0 i4 s' q& Z4 L, i( l
  b. w! p+ K1 Y%% 导入数据
* g' T5 r) u9 e  xclc, clear, close all
8 z2 E$ @$ {* ]- f+ P3 g/ V% clc:清除命令窗口的内容,对工作环境中的全部变量无任何影响
. E* b  x) E) ^% clear:清除工作空间的所有变量
$ _/ l: h0 c/ A3 E9 N# K- n9 f% close all:关闭所有的Figure窗口6 P7 j9 x7 P8 b; Z+ W/ \0 O8 I

, v$ L1 _! F+ g( B7 n% 导入数据) Q9 H1 `* s1 i
[~, ~, raw] = xlsread('sz000004.xlsx', 'Sheet1', 'A2:H7');) H+ r4 \' M* Q$ Y
% [num,txt,raw],~表示省略该部分的返回值) W* P- M7 s! q7 ~! g" A
% xlsread('filename','sheet', 'range'),第二个参数指数据在sheet1还是其他sheet部分,range表示单元格范围: _- e1 `, y) t4 W
2 k2 T9 z: U3 m- D" u; H! H' v3 ?& a
% 创建输出变量
$ z: K; i# c# ]: }" ?$ Ndata = reshape([raw{:}],size(raw));" |, U% K! c$ S! V; j# ^9 |  ~- ~4 ]
% [raw{:}]指raw里的所有数据,size(raw):6 x 8 ,该语句把6x8的cell类型数据转换为6x8 double类型数据
" @& B6 y. a2 L+ h" c) |& u/ W
% 将导入的数组分配列变量名称0 T; p& D1 d% d$ o) R$ ^
Date = data(:, 1); % 第一个参数表示从第一行到最后一行,第二个参数表示第一列
- O9 f- W5 B* w6 y; G* kDateNum = data(:, 2);) v9 ?; y8 G* `& `8 W
Popen = data(:, 3);* I* w% B: S! X# d- b
Phigh = data(:, 4);+ d. N7 s, r! W$ y1 u# d9 V$ |8 @* w+ X
Plow = data(:, 5);( W/ g. {0 A5 ]6 M
Pclose = data(:, 6);  
: K3 j7 N2 }: }Volum = data(:, 7); % Volume 表示股票成交量的意思,成交量=成交股数*成交价格 再加权求和
' f* Z4 X" J" r2 Q) ]- sTurn = data(:, 8); % turn表示股票周转率,股票周转率越高,意味着该股股性越活泼,也就是投资人所谓的热门股
$ R9 Y4 R! z! D, p5 d: B
0 x: b0 G4 k8 P+ W% 清除临时变量data和raw4 q( M8 m; Y# v* H. L  J6 {
clearvars data raw;: M1 t/ e' B( A

  F6 M3 z' f3 g9 l" ?%% 数据探索
8 W: C% ^  g. D$ \; U6 Q) r+ ?1 y
1 Q  q: ^' E2 _% ~9 |figure % 创建一个新的图像窗口
; O& [5 ^6 f  n, fplot(DateNum, Pclose, 'k'); % 'k',曲线是黑色的,打印后不失真% `) F3 t% b6 }( [* `7 [3 k1 A, Y
datetick('x','mm-dd'); % 更改日期显示类型。参数x表示x轴,mm-dd表示月份和日。yyyy-mm-dd,如2018-10-27+ H- B( _, Z! S3 O# N5 i
xlabel('日期') % x轴
6 d. m) s; ]/ W7 N  j' Iylabel('收盘价') % y轴- g& _3 `9 B9 r* L5 d/ {: H, z! {
figure
& g* D% h# m4 G& y# {) Ebar(Pclose) % 作为对照图形3 b5 B7 R! g1 w# a# F% b
: T$ U4 d2 |  E& T1 O$ {4 ^
%% 股票价值的评估
+ _; w& {! y' g* q: W
3 u% x: C9 x; V4 Pp = polyfit(DateNum, Pclose, 1); % 多项式拟合4 S! F& l! l/ F+ J' r
% polyfit()返回阶数为 n 的多项式 p(x) 的系数,p 中的系数按降幂排列
% f* j( K% a' o9 GP1 = polyval(p,DateNum); % 得到多项式模型的结果, u% r& Y) }% K
figure
9 X) m8 x! {5 S4 Yplot(DateNum,P1,DateNum,Pclose,'*g'); % 模型与原始数据的对照, '*g'表示绿色的*
+ _: }6 U8 x- ?# z* pvalue = p(1) % 将斜率赋值给value,作为股票的价值。p(1)最高项的次数
& ?4 `* W. U: |# Y; N* I* Y/ l  q% z3 _6 d, ]5 W' y
%% 股票风险的评估0 y. P. d; y  G9 j) t
MaxDD = maxdrawdown(Pclose); % 计算最大回撤
/ Z% E% h( v6 s/ L7 f. _  c1 lrisk = MaxDD  % 将最大回撤赋值给risk,作为股票的风险$ O, T/ ~8 G% r1 \# n, |* x
  3、回归算法演练。
; |( c3 l$ H; z. x
, t- i/ W) h5 [/ U' c8 C(1)一元线性回归" B; f6 X6 x' g& T' c$ {! G

  c& ?: a* C5 P7 {- d0 l- R" f# ?[ 例1 ] 近 10 年来,某市社会商品零售总额与职工工资总额(单位:亿元)的数据见表1,请建立社会商品零售总额与职工工资总额数据的回归模型。
" n5 a0 F5 b  w9 o9 Q3 _% ]8 o
% O/ S0 d9 p, k, Y
3 e# N: O  O7 N7 j, Z  d' ~' Y9 a  J9 c! z
该问题是典型的一元回归问题,但先要确定是线性还是非线性,然后就可以利用对应的回归方法建立他们之间的回归模型了,具体实现的 MATLAB 代码如下:
9 S1 j" M# w# g! H3 G5 A1 X7 I9 W4 ]
(1)输入数据
+ a) H7 F8 T7 q6 j% Q% Z9 e" V& z- m1 `& H2 \8 Z( O" l; Q: t
%% 输入数据4 _; F1 {; b) r2 G6 ]  I/ s% P$ A
clc, clear, close all1 z3 n3 g8 \. d# D) k$ z# ?
% 职工工资总额0 i' J, C6 I3 B2 R
x = [23.8,27.6,31.6,32.4,33.7,34.90,43.2,52.8,63.8,73.4];( R4 h! |) l6 t% u! {- |
% 商品零售总额
0 S# c$ ]. t3 c) ^y = [41.4,51.8,61.7,67.9,68.7,77.5,95.9,137.4,155.0,175.0];$ }8 N. E5 Q7 Y4 N: @8 w
(2)采用最小二乘回归* t8 J0 {4 s# |
- Y* v. ~% `& f4 A& l) f7 p0 Y% Z" t
%% 采用最小二乘法回归- \9 A, r1 |) t
% 作散点图
/ p8 v( m- T* s5 x* ffigure
6 i" k+ E3 y" E$ t& Yplot(x,y,'r*') % 散点图,散点为红色: ?: x9 A% ?# [
xlabel('x(职工工资总额)','fontsize',12)' _$ r; L) P1 N8 C8 X/ |: t4 z. e6 a
ylabel('y(商品零售总额)','fontsize',12)
, R. O" c9 p) ]set(gca, 'linewidth',2) % 坐标轴线宽为2
8 y, I4 y. U* _" }# Z4 S* |( a# O, v: R$ y
% 采用最小二乘法拟合) Q: o, ^6 s  F# [( i3 P5 }: ]2 v4 y) y
Lxx = sum((x-mean(x)).^2); %在列表运算中,^与.^不同
* K- Z& T+ U8 mLxy = sum((x-mean(x)).*(y-mean(y)));# q, N1 k7 h- K/ z
b1 = Lxy/Lxx;
% g% H9 C) B' t5 P- {% pb0 = mean(y) - b1 * mean(x);
) a  A/ p. _% v3 s0 x4 vy1 = b1 * x + b0;2 S4 u: q7 B( b& M; R9 s
; Z1 S! ^* X! a8 b5 z: K  a
hold on % hold on是当前轴及图像保持而不被刷新,准备接受此后将绘制的图形,多图共存
5 J, H6 f' s& A) i3 R4 v* G2 l! Uplot(x,y1, 'linewidth',2);4 r4 L4 T5 \( I( i
运行本节程序,会得到如图5所示的回归图形。在用最小二乘回归之前,先绘制了数据的散点图,这样就可以从图形上判断这些数据是否近似成线性关系。当发现它们的确近似在一条线上后,再用线性回归的方法进行回归,这样也更符合我们分析数据的一般思路。
, T+ B# T+ c, i. @, I9 k8 b
0 ]8 R, @% `( V- c8 T2 b! o6 z2 b
' \* s1 f$ d) V+ l2 W5 I/ V4 e3 @5 [" J+ f% B
                                                                                                    图5
: _7 ?6 v6 ~" T( J$ ]# e' [3 r. u$ o. K3 \/ R  v6 s
(3)采用 LinearModel.fit 函数进行线性回归; I" X$ r6 S% D: B  J; ?. ~/ g: T
, i3 g9 h1 w  M5 x
%% 采用 LinearModel.fit 函数进行线性回归
* m! _  X6 x( p' j- u' R5 O; Km2 = LinearModel.fit(x, y)  G' W+ n" X; `" v. L7 X
运行结果如下:
* N+ y7 v2 E5 E) a& q9 W' E
1 I, o6 c4 z5 P. jm2 =) J) _% m2 o6 P( w+ a5 U: Y. a
* z; N: J1 }5 C7 E; W) I7 o. \
Linear regression model:
+ [0 g  y6 f% J3 A
. X/ u8 U$ C0 P" u) \    y ~ 1 + x1" u7 A3 j* p' K7 w" m: E7 f4 C
Estimated Coefficients:
$ i: v2 v+ }( j1 o% y# F& G: ?2 X. H% l
               Estimate      SE       tStat       pValue " R& Z+ s8 n$ K6 o. f; }
& t: w$ g2 X! e- F8 E
    (Intercept)    -23.549      5.1028    -4.615     0.0017215
* ?  J0 S: d" J( s, i4 f
8 I% P( M6 c/ l    x1           2.7991     0.11456    24.435    8.4014e-09  k% |$ ?& h$ L

, A; L: m2 Y  A$ JR-squared: 0.987,  Adjusted R-Squared 0.9859 C% o# B$ h) ?6 q1 H1 _

7 |+ b/ F. V+ M( s$ [. r+ V4 M( YF-statistic vs. constant model: 597, p-value = 8.4e-09. b7 \/ }% ~' @, A- P8 x

* L8 z: ?  Y; @' s& v/ Q1 c, a1 R如下图,我们只需记住-23.594是一次函数的中x的系数,2.7991是一次函数中的常数项即可,其它的不用理会。6 g& Y0 I" U6 Q$ z; n. Y1 b' S( N& w

" s% m8 j: B) _4 B% H
5 U9 y0 K6 g5 c4 q
/ A8 d" T' {$ k" U! @4)采用 regress 函数进行回归
. i9 D8 k( e5 B2 d  u
7 e3 ~- a( \/ I; R- D) U& J%% 采用 regress 函数进行回归
- ]- O9 F8 X# n2 vY = y'
( ~4 z9 e: K! ^8 O, z' nX = [ones(size(x,2),1),x']2 J/ D2 j; o* y1 y
[b,bint,r,rint,s] = regress(Y,X)* L& E5 v, u9 l' `
运行结果如下:# G& I4 j9 r# R/ N! l0 W
* b) d1 N+ _$ E4 g) g
b =, k& X6 ]9 `4 B! J4 N' D
6 l. n% e8 N7 |
  -23.5493  I8 b# m8 D8 e* W5 p" ^+ N

5 B; B- ^( O* D2 l' E! x    2.7991, k9 j8 N/ ^) q* G* _$ F3 ^7 D
9 y3 F  w4 U* `1 U' c+ s, |9 V
我们只需记住-23.594是一次函数的中x的系数,2.7991是一次函数中的常数项即可,其它的不用理会。+ M; l! t" P1 E+ y2 j1 u5 O( R

% k; ]; T) o0 ~9 k7 g(2)一元非线性回归* J/ {; L( p7 ?" Z
9 Z3 z/ G( ]3 {& b
[ 例2 ] 为了解百货商店销售额 x 与流通率(这是反映商业活动的一个质量指标,指每元商品流转额所分摊的流通费用)y 之间的关系,收集了九个商店的有关数据(见表2)。请建立它们关系的数学模型。
/ q0 g. |6 ~- A( Q. p7 X
$ f( ^. x% |" i$ q, s9 _, p# j0 I6 Q# [5 J3 t) e
; K; J2 P6 S" Z. P/ j
/ A. m  ]. @! C+ `
% s; c/ K3 E: M- J# ~1 I8 x6 t0 r
        为了得到 x 与 y 之间的关系,先绘制出它们之间的散点图,如图 2 所示的“雪花”点图。由该图可以判断它们之间的关系近似为对数关系或指数关系,为此可以利用这两种函数形式进行非线性拟合,具体实现步骤及每个步骤的结果如下:
9 b0 @# e, T# G
2 T9 F3 o4 B6 f; \(1)输入数据
/ S  {& [; N6 X, E# }& p0 A; R& l0 ~3 d
%% 输入数据. ?9 k0 ?+ g( k2 j- l
clc, clear all, close all
& C3 M0 U( P! s7 R8 d; {1 t3 cx = [1.5, 4.5, 7.5,10.5,13.5,16.5,19.5,22.5,25.5];
# O2 }" L" K1 T+ k; O: F8 i5 by = [7.0,4.8,3.6,3.1,2.7,2.5,2.4,2.3,2.2];
: {* n, `! V; d3 K! {plot(x, y, '*', 'linewidth', 1) % 这里的linewidth指的是散点大小9 p! v9 F  r, [; v6 \; k
set(gca,'linewidth',2) % 设置坐标轴的线宽为2/ `7 L" l) a! n
xlabel('销售额x/万元','fontsize',12)
$ M% J9 o4 ]4 Jylabel('流通率y/%','fontsize',12)
9 J' ]$ j( z$ ~(2)对数形式非线性回归1 K/ |  Q( ]$ Y6 Z) u# K- B, m& s
% G/ S/ r8 c4 m6 `$ t
%% 对数形式非线性回归/ z* o# U5 D/ N' t
m1 = @(b,x) b(1) + b(2)*log(x);
! _$ q9 }- f3 Dnonlinfit1 = fitnlm(x,y,m1,[0.01;0.01])
! l- q* T0 N0 W" e! f4 M' Yb = nonlinfit1.Coefficients.Estimate;
4 v1 A9 I) P2 }! a4 i1 f: w% BY1 = b(1,1) + b(2,1)*log(x);
4 A6 i: I8 P* m) phold on , T/ N) j4 J2 V0 j) Y( N# V
plot(x, Y1, '--k', 'linewidth',2)
8 p* M) T% M1 V6 a3 Q运行结果如下:8 }% v# b; i) T) A* W" @0 h
6 y% U: r6 `2 C. j
nonlinfit1 =( j' D! n3 H; g0 |; Y7 ^9 Z
: e. N2 Y/ `1 g/ p) v  |# ~
Nonlinear regression model:3 M) H( z0 D3 s: t1 A5 m

0 x2 H' s% ?6 M+ c0 G$ m# k    y ~ b1 + b2*log(x)0 _) ?. m7 T1 z7 b

3 l# S' R+ t8 @3 oEstimated Coefficients:5 [7 W; {( }) `) k# a1 T. T

1 N/ B1 x, }* L+ ]4 z          Estimate      SE        tStat       pValue
5 l3 H* ^0 n5 [( A9 u
7 h+ o% v; A4 \0 \# b    b1    7.3979      0.26667     27.742    2.0303e-086 y2 c7 p; U" d* H

0 y* q8 J, N5 D  V    b2    -1.713      0.10724    -15.974    9.1465e-07
, H& J  X" g5 y; C; S! F9 j1 }/ i% W; z$ s  L9 J
R-Squared: 0.973,  Adjusted R-Squared 0.969
1 ]: M- e1 M* s% t& J; M& M- h# J+ W3 b  @, r: W: a# R: S
F-statistic vs. constant model: 255, p-value = 9.15e-07
8 B* c2 l( }8 p. I8 `$ K7 q/ P4 T3 B; D+ [  h3 _
(3)指数形式非线性回归# g& H9 F4 I' q- Q- e4 C8 x

' i: o: _# j2 a0 k7 N%% 指数形式非线性回归
7 a* r2 K9 Z$ b: M4 V) ?m2 = 'y ~ b1*x^b2';! |" _. M, Z% |4 i
nonlinfit2 = fitnlm(x,y,m2, [1;1])
6 N* E; d1 h: u1 W" ub1 = nonlinfit2.Coefficients.Estimate(1,1);3 N6 F; R; q! B. f. }! D' J; e
b2 = nonlinfit2.Coefficients.Estimate(2,1): V. |, D, r* ~- w3 k
Y2 = b1*x.^b2;" i& r; M' G5 v! a# G8 w
hold on;
( T. b9 M/ N9 y4 A3 {plot(x,Y2,'r','linewidth',2)
# y) I3 o( y( A. X( S1 z! Olegend('原始数据','a+b*lnx','a*x^b') % 图例# u* Q# ~8 G: G( x/ p
运行结果如下:
& @; {2 o0 T8 @+ ^1 P/ @6 N) o; T+ \) Q+ c( V6 T$ }8 H
nonlinfit2 =
  K$ |! v5 b4 E) p5 n' W- N; b; a
" ~$ y* f/ i8 P; W7 e1 \' gNonlinear regression model:
; e( U1 k5 _. l! l, M3 b
( i7 p0 Z* p/ o! W$ F, ~    y ~ b1*x^b2
: e$ N: ?! y4 k/ o; C
: l& n5 B3 c; W- r- CEstimated Coefficients:) {( c8 L; I5 N
! |1 G* a' {* L0 W
          Estimate       SE        tStat       pValue
% N$ F% d  o! N. u4 U
" l( i7 [1 g8 K  Y" E0 _    b1      8.4112     0.19176     43.862    8.3606e-10
0 f: n3 G+ b! i0 H
2 I- o4 m$ v" D( J4 ]4 |, E    b2    -0.41893    0.012382    -33.834    5.1061e-09& ?7 `, _9 `: H2 X: p

5 h! t+ [2 V$ S& d( O: mR-Squared: 0.993,  Adjusted R-Squared 0.992% s5 Z% B5 U$ d" I

+ ^2 z7 ^4 a& H* ]( v. L0 K9 }F-statistic vs. zero model: 3.05e+03, p-value = 5.1e-115 P; o& t# l. N/ `+ y

; B1 j6 M- n% V5 ?' S在该案例中,选择两种函数形式进行非线性回归,从回归结果来看,对数形式的决定系数为 0.973 ,而指数形式的为 0.993 ,优于前者,所以可以认为指数形式的函数形式更符合 y 与 x 之间的关系,这样就可以确定他们之间的函数关系形式了。! N/ W; p6 d8 F+ H
0 E7 z* m1 F! E( _' x% n9 k( K
2.多元回归
( ~3 U# D' p1 t/ m" n) Y- u! X* ^2 _) |9 H
1.多元线性回归3 v' O6 E8 b3 c( E* H1 e& K5 a
6 J) O7 }4 F/ _
[ 例3 ] 某科学基金会希望估计从事某研究的学者的年薪 Y 与他们的研究成果(论文、著作等)的质量指标 X1、从事研究工作的时间 X2、能成功获得资助的指标 X3 之间的关系,为此按一定的实验设计方法调查了 24 位研究学者,得到如表3 所示的数据( i 为学者序号),试建立 Y 与 X1 , X2 , X3 之间关系的数学模型,并得出有关结论和作统计分析。+ r6 o7 Q: F2 R! _
2 z2 ]' R) [3 t/ r* W0 S% {
+ w" {$ o# ]8 n' g9 h* @- I
" l- y1 C0 I+ H" _
该问题是典型的多元回归问题,但能否应用多元线性回归,最好先通过数据可视化判断他们之间的变化趋势,如果近似满足线性关系,则可以执行利用多元线性回归方法对该问题进行回归。具体步骤如下:5 ?2 g. o* T2 X& i6 g# {
8 N. f/ q; g8 j
(1)作出因变量 Y 与各自变量的样本散点图
1 A+ t3 B4 n% E* j0 @. Q9 A; H6 c3 `5 I% s
作散点图的目的主要是观察因变量 Y 与各自变量间是否有比较好的线性关系,以便选择恰当的数学模型形式。图3 分别为年薪 Y 与成果质量指标 X1、研究工作时间 X2、获得资助的指标 X3 之间的散点图。从图中可以看出这些点大致分布在一条直线旁边,因此,有比较好的线性关系,可以采用线性回归。绘制图3的代码如下:
. m' A' P8 `7 }% G# o  z& C
8 o  Z! u' v/ t%% 作出因变量Y与各自变量的样本散点图1 Y1 T/ g, t2 F+ L3 F
% x1,x2,x3,Y的数据
* v2 X" J9 d3 U' M  gx1=[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];1 [, A6 N* B( o3 x
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];
9 [; ~  `  r! }/ L) L+ tx3=[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];* \* q5 a- d, H  h. |6 K$ ?
Y=[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];; |+ C2 c# G' o: Y
% 绘图,三幅图横向并排7 h; v/ [9 y( d+ ?! u
subplot(1,3,1),plot(x1,Y,'g*')# W! t$ m* b" K" L* A6 C4 f9 P
subplot(1,3,2),plot(x2,Y,'k+')1 i# T6 a7 ~1 X/ q) Q# u! C6 Q; f
subplot(1,3,3),plot(x3,Y,'ro')2 B7 k- S* B% H. g" \+ K: u. O9 t5 C
绘制的图形如下:
- h0 t9 @% M! d, B& D
2 w2 a% N/ C0 ]* J2 a7 W, Y& X1 x9 F
+ y0 g% X0 F1 H+ A4 Q' R
(2)进行多元线性回归! R% L! w& A9 a

! b7 `  p+ N0 o& W这里可以直接使用 regress 函数执行多元线性回归,注意以下代码模板,以后碰到多元线性问题直接套用代码,具体代码如下:
6 W9 \% e8 v7 A4 _, y5 T+ D6 @2 W+ P' S& j' f, ~; L% l. D& w
%% 进行多元线性回归
3 F: F' {6 R8 o- G5 U, \n = 24; m = 3; % 每个变量均有24个数据,共有3个变量" h( v! J( \* ?: c/ ?9 T4 r
X = [ones(n,1),x1',x2',x3'];
7 ~2 B! V. w7 C+ S- H: t4 I[b,bint,r,rint,s]=regress(Y',X,0.05) % 0.05为预定显著水平,判断因变量y与自变量之间是否具有显著的线性相关关系需要用到。  Y; ]/ ~- m' l( i# y# f# h
运行结果如下:
. F4 h/ u) ?4 l- v) l" a. v2 B
! ]2 @! r+ h% D) @$ Z3 Pb =, V  F6 W7 K; [- i

2 c% K& m7 d7 ]: f9 S- v   18.01579 ^  H0 x: {1 W
    1.0817
. O/ n/ x! b, R: B3 M7 w    0.3212
2 B7 L5 Q" Z, L7 p. }( q8 u    1.2835
4 L; I7 `. V. P0 _* G! D" [
9 m: Z3 R2 v1 {6 T* K5 p
% H2 {. ~0 e% X3 ]  ?; b) k6 r! Fbint =6 n* h( q( I4 C. E/ O) B

* t$ r  m+ A1 U8 z, e) t   13.9052   22.1262* `8 Q: A1 p! n6 D. A2 Y/ z
    0.3900    1.7733. R. ?* C: F: b0 i
    0.2440    0.3984) q0 h2 P- `- e* q
    0.6691    1.8979
5 q+ L8 G9 w% v( J" C" L& I: H- S0 a& h" l5 q% [

. C, J) a( w  Tr =
+ h) d9 K1 r* I% O+ z/ t8 F) ]% {" b( E$ u9 m- |
    0.6781( j3 g. ^% i: i4 V
    1.9129
2 W! h$ j0 ^% U. _- a- @! w$ H! E   -0.1119
7 s* H$ D. O' t+ Z# ]    3.3114& i( l; p" ?+ }4 S- p' o" p2 I) ]
   -0.7424
& k. N6 c! f; |& S6 ^    1.2459
# k0 ~% o1 [1 h   -2.1022
, w. z$ @" ^6 L) s    1.9650
$ H9 n9 Q% K* ~1 S   -0.3193
3 [0 s8 f, H3 K' D* H* D5 D    1.3466
# {$ c  e: v: U) A% K3 m    0.8691
7 f% D$ ~5 a. G   -3.2637/ l/ z! M7 s" m  C: ^* D
   -0.5115
; B7 r! X8 T# u0 `9 a0 P   -1.1733
( Q! `; h: c1 }   -1.49105 G4 z+ d. {8 F+ q5 ~
   -0.2972
/ g* {; ]. F/ u0 Q, T    0.17026 d+ e8 j7 F$ A1 [
    0.5799
# o4 G- E9 T2 [   -3.2856: t5 _8 s. w+ ~/ \/ V4 S% |; d
    1.1368
. b+ b: w/ U1 r9 h7 @9 I) D   -0.8864
, U2 W; p: J/ i% o; ]3 C   -1.4646
) y3 P% u# E' D9 z, S, G    0.8032! s7 ^0 `4 j7 _  g! x
    1.63012 A  Q) I7 v' H8 G) K3 ~; |
: C! i2 M+ a- ^9 h6 k9 o: I  G

' m! D) M9 S0 ]/ Brint =
1 `" I1 ?( B- x4 x) s3 G
) z8 a8 q0 [0 g   -2.7017    4.0580+ X# T2 `7 Y  Y: X
   -1.6203    5.4461
2 Y- R# }/ g/ L+ y   -3.6190    3.39517 H' `5 O2 n8 J9 @7 [& a) g
    0.0498    6.5729
$ X" h( y2 j/ j4 B. v4 b: b" u   -4.0560    2.5712
0 X, g( v6 t# y" \+ c( S! s2 f   -2.1800    4.6717. Z5 e/ K, \. t  \9 k0 R
   -5.4947    1.2902
0 x! A  F( |. o' |   -1.3231    5.2531
0 h0 E1 [9 x+ V( c   -3.5894    2.9507
% v4 [% M" k! o' ^# C   -1.7678    4.4609
% V5 T8 p. x% g5 @2 z   -2.7146    4.4529
5 z) Q/ t$ Z" _" l   -6.4090   -0.1183: e; X: o  `( f- D3 ?; w
   -3.6088    2.5859
* F: U  D. ^2 d, _/ s% b   -4.7040    2.3575
. |- S, w2 [+ t& b# M! ]: j1 b" W   -4.8249    1.8429
# ]8 n$ Q! Q5 o) {! E) c% d7 b   -3.7129    3.1185; R& N/ a. S- h, L6 R4 r: D- H
   -3.0504    3.39075 o2 e/ ~& P& ]
   -2.8855    4.0453
8 h7 e! z% p" C, @. x' \) }0 Z   -6.2644   -0.3067  U; _: B) R7 @, b! z0 S" V
   -2.1893    4.4630  u: A, L' X& v2 t, E
   -4.4002    2.6273
3 P+ M' ]5 `' }# }+ R   -4.8991    1.9699
* e* M2 c# U+ t; s   -2.4872    4.0937
* N& @4 K8 _% e! X) ^- z: ~( K& S   -1.8351    5.0954
9 V1 \1 c3 X  x, w3 Q
7 v1 Y4 t1 ~6 Z6 v1 ~" j) N
1 T$ s* f7 }' X" Js =
0 P3 H  _1 `, u3 k0 ~6 _+ G. e* \. X/ S4 }0 `) S
    0.9106   67.9195    0.0000    3.0719
& q$ h: t% Z! {6 J9 |0 L6 j. x看到如此长的运行结果,我们不要害怕,因为里面很多数据是没用的,我们只需提取有用的数据。9 w  ^3 i' `8 x, @
# }% o$ i0 D. M
在运行结果中,很多数据我们不需理会,我们真正需要用到的数据如下:9 y1 x1 }2 ]1 K9 _5 L* f
8 j* C# x8 t+ ^
b =! s9 {( Y8 ]6 Z* A; f
5 E% _3 z/ x& u+ C0 y: ]
   18.0157
3 O1 K( ?! A! W) f    1.08175 C1 w  P) |- i6 k$ u- o: t
    0.32129 M% f6 e2 k3 K" d4 k
    1.28359 p0 I$ g. v2 H! @9 c5 S
2 d7 z6 L5 X& {+ ?7 i
s =
: A! t) v6 G7 L" o2 U) I/ o6 U& Q  Z3 d6 j( u2 X; O3 P% v
    0.9106   67.9195    0.0000    3.0719: D, x: Q; L; ^5 Q
回归系数 b = (β0,β1,β2,β3) = (18.0157, 1.0817, 0.3212, 1.2835),回归系数的置信区间,以及统计变量 stats(它包含四个检验统计量:相关系数的平方R^2,假设检验统计量 F,与 F 对应的概率 p,s^2 的值)。观察表4的数据,会发现它来源于运行结果中的b和s:" k$ i& P# _) x  n3 E

: M5 a# \0 @$ f8 U7 p/ Y8 I
( L, A0 `% t  S( M8 D* {) O2 W3 X& G6 Q+ d( v
根据β0,β1,β2,β3,我们初步得出回归方程为:
7 j! h1 \, }! F6 O; h4 }" _3 ?: A. d1 ^7 r

3 \7 n+ e6 b3 ]4 C$ N; O. @1 D0 S" Y& y. W& F6 y: F: T
如何判断该回归方程是否符合该模型呢?有以下3种方法:  m# i) Q$ l" o- `9 T  L& J

1 ~* E# z1 j4 ^, g4 t1)相关系数 R 的评价:本例 R 的绝对值为 0.9542 ,表明线性相关性较强。* R' E: H/ n* ~) u- N* |& v

. j0 D, o# ^9 O( ?( k% g$ Z2)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。
% U' {) _. b1 c8 Z- O% }. m+ w0 V, M" Y& j. A. t3 j! ~) N
3)p 值检验:若 p < α(α 为预定显著水平),则说明因变量 y 与自变量 x1,x2,...,xm之间显著地有线性相关关系。本例输出结果,p<0.0001,显然满足 p<α=0.05。
8 t& Y7 u3 i" x3 E. T4 \$ d0 A# U4 ~$ X
以上三种统计推断方法推断的结果是一致的,说明因变量 y 与自变量之间显著地有线性相关关系,所得线性回归模型可用。s^2 当然越小越好,这主要在模型改进时作为参考。
" j# ^3 R$ W3 X% L
3 a# D' K: k% Y3. 逐步回归! L" r" {' h' x( W3 r. [$ k
! g' d0 m6 |' S1 L* Z- g
[ 例4 ] (Hald,1960)Hald 数据是关于水泥生产的数据。某种水泥在凝固时放出的热量 Y(单位:卡/克)与水泥中 4 种化学成品所占的百分比有关:' b+ q/ X8 `+ z! \# N  l
/ }, p" K" l9 x1 R) y0 [  |0 y2 O7 M, H

5 S; N$ z& c7 e9 x. Y' F* }) ?9 w  m" K! ]% `
在生产中测得 12 组数据,见表5,试建立 Y 关于这些因子的“最优”回归方程。
) n+ M9 }; Z+ e+ N" {$ h$ e4 Q  g0 ~9 Y' e$ A) z
# I% z" @/ x$ ^* S# A7 x
& O; _. E% A8 f" q. ?9 e6 G
对于例 4 中的问题,可以使用多元线性回归、多元多项式回归,但也可以考虑使用逐步回归。从逐步回归的原理来看,逐步回归是以上两种回归方法的结合,可以自动使得方程的因子设置最合理。对于该问题,逐步回归的代码如下:& \3 I0 m5 G# U) }& Q) ?
1 i, |, i$ \) M0 P
%% 逐步回归
/ g; S( e' `' I: v, n3 @( A3 L' e; ]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];   %自变量数据
* r0 r. [$ ^, K# V0 {. ]) m3 WY=[78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3];  %因变量数据
% c- H+ r# p; ^# pstepwise(X,Y,[1,2,3,4],0.05,0.10)% in=[1,2,3,4]表示X1、X2、X3、X4均保留在模型中1 q/ w$ ]) C& Y3 s. j2 {
程序执行后得到下列逐步回归的窗口,如图 4 所示。1 W& w5 ]& [  n/ A2 b

1 k, z( p" t( b* d$ F0 X
0 Y$ ]( i) i8 z7 e2 @3 z
6 Z1 i/ ^( }0 J3 _4 [                                                                                                             图41 J7 S: [0 b) w/ ^+ D3 F6 E. d

" N) ]' z' S( V; t$ |$ D8 ?在图 4 中,用蓝色行显示变量 X1、X2、X3、X4 均保留在模型中,窗口的右侧按钮上方提示:将变量X4剔除回归方程(Move X4 out),单击 Next Step 按钮,即进行下一步运算,将第 4 列数据对应的变量 X4 剔除回归方程。单击 Next Step 按钮后,剔除的变量 X3 所对应的行用红色表示,同时又得到提示:将变量 X3 剔除回归方程(Move X3 out),单击 Next Step 按钮,这样一直重复操作,直到 “Next Step” 按钮变灰,表明逐步回归结束,此时得到的模型即为逐步回归最终的结果。最终结果如下:+ [) |. y9 k1 j7 i( Q( u# L

  \) t2 O) X# c/ q* b/ S) f
' B, I/ e# x* {2 N4 s( _3 v4 j0 d$ D! P, H; X0 D
4. 逻辑回归
6 E% f* S& i) H7 G& w" G6 z9 c6 a& R  l# B( J3 o  I% E
[ 例5 ] 企业到金融商业机构贷款,金融商业机构需要对企业进行评估。评估结果为 0 , 1 两种形式,0 表示企业两年后破产,将拒绝贷款,而 1 表示企业 2 年后具备还款能力,可以贷款。在表 6 中,已知前 20 家企业的三项评价指标值和评估结果,试建立模型对其他 5 家企业(企业 21-25)进行评估。
" X8 a  R7 b: [+ K0 g
, l" t1 T5 u1 }
" V( t4 P% B" M; r8 @" V0 y
0 Q* j/ {! }- b4 E- `! a' p- I7 D9 T对于该问题,很明显可以用 Logistic 模型来回归,具体求解程序如下:! p+ W  q; D4 B2 n: [+ j$ T* J

3 D* g$ [' B+ E% E  L: i+ s/ }程序中需要用到的数据文件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, G5 D% g$ e7 A0 ^7 p- h$ n% [
! ]0 ?% Z' I9 D2 c* \! t  s/ a
% logistic回归
1 C& s$ C* s. G6 }/ e; L8 {/ B7 A+ X4 L
+ N/ P* y3 _  N1 B# r%% 导入数据
, r# U. C, @& J- s( B3 lclc,clear,close all
# A  k: x/ }# cX0 = xlsread('logistic_ex1.xlsx','A2:C21'); % 前20家企业的三项评价指标值,即回归模型的输入  c8 ]1 `& O# u
Y0 = xlsread('logistic_ex1.xlsx','D221'); % 前20家企业的评估结果,即回归模型的输出6 X$ q- f$ W2 Q: F
X1 = xlsread('logistic_ex1.xlsx','A2:C26'); % 预测数据输入
+ A5 a9 T% l; N: A2 n, S6 y) ]1 ~7 N: X# W
%% 逻辑函数1 L, v1 {8 c- k9 B& U
GM = fitglm(X0,Y0,'Distribution','binomial');* j# Z; K$ c" g  m( t; R! h& R
Y1 = predict(GM,X1);
( R6 j; Y2 x# A7 v. l0 D8 j" n7 U7 \( S0 j0 E' D
%% 模型的评估
. p2 ]0 D- M3 fN0 = 1:size(Y0,1); % N0 = [1,2,3,4,……,20]
3 A1 b  B, z  x' P; A; Q7 p4 YN1 = 1:size(Y1,1); % N1 = [1,2,3,4,……,25]
* N6 R/ `! f: U8 Bplot(N0',Y0,'-kd'); % N0'指的是对N0'进行转置,N0'和Y0的形式相同,该行代码绘制的是前20家企业的评估结果
& [0 p1 W3 d- z0 u  e% plot()中的参数'-kd'的解析:-代表直线,k代表黑色,d代表菱形符号
, u, R- W% g7 ^" G" bhold on;! j: C0 K% J; h! ^/ A
scatter(N1',Y1,'b'); % N1'指的是对N1'进行转置,N1'和Y1的形式相同
$ t, A. T4 I; |* S/ Zxlabel('企业编号');; ]# _3 ]; U8 P' ?/ X% ]
ylabel('输出值');
7 G5 F" c: F0 {2 n2 h/ o得到的回归结果与原始数据的比较如图5所示。0 \3 O; W$ f' i7 G6 I
3 e& W6 k( N/ r1 m" g2 h

/ S( t2 J5 {% E& z& s6 X1 T& W  L% W/ r
- j& O8 z2 [, g1 W                                                                   图54 u* S5 R8 s' A( I+ e/ [

; F- ~, U7 w# Q' {" l三、总结与感悟。
, H  ]5 e9 S4 {5 U9 V
* V; m0 n1 N' V  o8 c        总结:通过这次学习,我了解到Matlab在数学建模竞赛中使用广泛;在评估股票价值与风险的小实例中,我掌握了用Matlab去建模的基本方法和步骤;在回归算法的学习过程中,我掌握了一元线性回归、一元非线性回归、多元线性回归、逐步回归、逻辑回归的算法。
* A( H# N# s& _2 X; r7 p1 G
' p5 @% o* R0 @  x, x0 m% t: _( @        感悟:正确且高效的 MATLAB 编程理念就是以问题为中心的主动编程。我们传统学习编程的方法是学习变量类型、语法结构、算法以及编程的其他知识,因为学习时候是没有目标的,也不知道学的知识什么时候能用到,收效甚微。而以问题为中心的主动编程,则是先找到问题的解决步骤,然后在 MATLAB 中一步一步地去实现。在每步实现的过程中,遇到问题,查找知识(互联网时代查询知识还是很容易的),定位方法,再根据方法,查询 MATLAB 中的对应函数,学习函数用法,回到程序,解决问题。在这个过程中,知识的获取都是为了解决问题的,也就是说每次学习的目标都是非常明确的,学完之后的应用就会强化对知识的理解和掌握,这样即学即用的学习方式是效率最高,也是最有效的方式。最重要的是,这种主动的编程方式会让学习者体验到学习的成就感的乐趣,有成就感,自然就强化对编程的自信了。这种内心的自信和强大在建模中会发挥意想不到的力量,所为信念的力量。/ U; p( J' g; v9 }7 f% ^  F

9 c, P4 o+ m1 Y2 s3 J9 q0 T! f3 u( J' I& [( f





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