数学建模社区-数学中国

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

作者: 杨利霞    时间: 2019-4-10 15:18
标题: Matlab数学建模学习报告(一)
Matlab数学建模学习报告(一)
: b. N+ z, U* C/ b0 t* I, m一、学习目标。

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

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

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


8 `( ?6 r( m1 O二、实例演练。
: ~# X8 M) _: _0 U, Q( i; Z1 i9 b% G# ]" N, i: j
   1、谈谈你对Matlab与数学建模竞赛的了解。( A3 V2 t9 Z. t
. S! N+ ]' U; ]1 e
        Matlab在数学建模中使用广泛:MATLAB 是公认的最优秀的数学模型求解工具,在数学建模竞赛中超过 95% 的参赛队使用 MATLAB 作为求解工具,在国家奖队伍中,MATLAB 的使用率几乎 100%。虽然比较知名的数模软件不只 MATLAB。( u- K; U8 {8 e2 K: g7 D  n7 O

; r. {% n2 @  u! `$ [4 [        人们喜欢使用Matlab去数学建模的原因:8 Q) k$ s* x$ w- w/ p

0 X4 w6 s" D2 [# W(1)MATLAB 的数学函数全,包含人类社会的绝大多数数学知识。
- _' d2 e  }, p( r' }- e5 b7 z0 }0 y6 K% u3 }
(2)MATLAB 足够灵活,可以按照问题的需要,自主开发程序,解决问题。; f5 X, [* f# p/ t

) e% Z) Z5 F2 M3 X2 W* I(3)MATLAB易上手,本身很简单,不存在壁垒。掌握正确的 MATLAB 使用方法和实用的小技巧,在半小时内就可以很快地变成 MATLAB 高手了。
* z: Q9 a, I: Y3 k$ `2 I* H( r# j7 T* V5 n5 a# U. I
        正确且高效的 MATLAB 编程理念就是以问题为中心的主动编程。我们传统学习编程的方法是学习变量类型、语法结构、算法以及编程的其他知识,因为学习时候是没有目标的,也不知道学的知识什么时候能用到,收效甚微。而以问题为中心的主动编程,则是先找到问题的解决步骤,然后在 MATLAB 中一步一步地去实现。在每步实现的过程中,遇到问题,查找知识(互联网时代查询知识还是很容易的),定位方法,再根据方法,查询 MATLAB 中的对应函数,学习函数用法,回到程序,解决问题。在这个过程中,知识的获取都是为了解决问题的,也就是说每次学习的目标都是非常明确的,学完之后的应用就会强化对知识的理解和掌握,这样即学即用的学习方式是效率最高,也是最有效的方式。最重要的是,这种主动的编程方式会让学习者体验到学习的成就感的乐趣,有成就感,自然就强化对编程的自信了。这种内心的自信和强大在建模中会发挥意想不到的力量,所为信念的力量。
% T( I5 j+ b4 V& `# u: d0 H0 i# P
         数学建模竞赛中的 MATLAB 水平要求:' P- c7 q- O$ u2 ~! b6 h
, j( A% V0 ~( N& g7 E
要想在全国大学生数学建模竞赛中拿到国奖, MATLAB 技能是必备的。 具体的技能水平应达到:
' k+ Y4 d$ _7 _' }* G8 r8 G4 x) N5 z" p: {4 W
1)了解 MATLAB 的基本用法,包括几个常用的命令,如何获取帮助,脚本结构,程序的分节与注释,矩阵的基本操作,快捷绘图方式;
0 V$ P2 K) p+ e1 y* e
" \& U5 I% W8 i9 T9 J2)熟悉 MATLAB 的程序结构,编程模式,能自由地创建和引用函数(包括匿名函数);6 }) J+ B) x: x

" |% n1 L5 _$ ?. Z3)熟悉常见模型的求解算法和套路,包括连续模型,规划模型,数据建模类的模型;
# Y; u( {# ^# A
. }' X5 T0 j: d4)能够用 MALTAB 程序将机理建模的过程模拟出来,就是能够建立和求解没有套路的数学模型。
1 K5 q3 e1 O3 ~' V: Y% m" Z( z- _" r7 T  {: p* ?. \- R0 y
要想达到如上要求, 不能按照传统的学习方式一步一步地学习, 而要结合上述提到的学习理念制定科学的训练计划。% ]4 a% B0 K8 n; f1 z
1 _. `6 U3 s9 ~, C) ^
  2、已知股票的交易数据:日期、开盘价、最高价、最低价、收盘价、成交量和换手率,试用某种方法来评价这只股票的价值和风险。如何用MATLAB去求解该问题?(交易数据:点击此处获取数据)! L& o* W7 M# b5 L
( Q. k5 a/ L+ F) B
解题步骤:
% X' Y7 o# b2 j9 a4 x
0 f  c7 V' Y  f$ y" n/ `0 J2 m第一阶段:从外部读取数据" V/ z$ U% f' K( f7 C" ~+ {

' |& C6 ^! |/ C: e) X1 l$ ?Step1.1:把数据文件sz000004.xls拖曳进‘当前文件夹区’,选中数据文件sz000004.xls,右键,将弹出右键列表,很快可发现有个“导入数据”菜单,如图 1 所示。/ u( W6 m1 R8 k% ?2 M( Y
! I" `" d" ^+ X! u+ [
1 g/ o  o3 T/ c5 S8 g9 ^, Y
; I6 ]% D, F# I1 b2 e
                                                                  图1. 启动导入数据引擎示意图) W- n/ L3 [' K1 t$ B- k
& S4 B1 {: ~1 E! Q* D
Step1.2:单击“导入数据”这个按钮,则很快发现起到一个导入数据引擎,如图 4 所示。$ @8 H; J* g8 a/ j. Y+ P% y# ]4 v/ U
4 A. |: R0 v' o( N  |" e

* m) J) X+ N+ @# @8 N# C& h- h. C3 K9 T
                                                                    图2. 导入数据界面* X) I6 B3 M. V9 U
) u1 ?. P9 g" B: x2 ~* m" U
Step1.3:观察图 2,在右上角有个“导入所选内容”按钮,则可直接单击之。马上我们就会发现在 MATLAB 的工作区(当前内存中的变量)就会显示这些导入的数据,并以列向量的方式表示,因为默认的数据类型就是“列向量”,当然您可以可以选择其他的数据类型,大家不妨做几个实验,观察一下选择不同的数据类型后会结果会有什么不同。至此,第一步获取数据的工作的完成。
# o2 _5 H* R! l% R2 X) `! S# ^# W- x( a% x
9 i8 o: \1 }0 R% @# Q& w) h2 q
( c, ~2 v( U7 P) `/ L
第二阶段:数据探索和建模9 \" H5 G0 s  b. s/ B$ L. }7 o. ]

: o8 P7 v% ~0 {) N4 Z/ ^# s现在重新回到问题,对于该问题,我们的目标是能够评估股票的价值和风险,但现在我们还不知道该如何去评估,MATLAB 是工具,不能代替我们决策用何种方法来评估,但是可以辅助我们得到合适的方法,这就是数据探索部分的工作。下面我们就来尝试如何在 MATLAB 中进行数据的探索和建模。( j+ s' t+ X* F' i" ~2 N, y

* K, W: D8 l0 _' Q/ E+ FStep2.1:查看数据的统计信息,了解我们的数据。具体操作方式是双击工具区(直接双击这三个字),此时会得到所有变量的详细统计信息。通过查看这些基本的统计信息,有助于快速在第一层面认识我们所正在研究的数据。当然,只要大体浏览即可,除非这些统计信息对某个问题都有很重要的意义。数据的统计信息是认识数据的基础,但不够直观,更直观也更容易发现数据规律的方式就是数据可视化,也就是以图的形式呈现数据的信息。下面我们将尝试用 MATLAB 对这些数据进行可视化。
! Y% N* F9 `7 ]1 V* B' _
! _- d1 \. z0 F0 x/ `6 ?3 }由于变量比较多,所以还有必要对这些变量进行初步的梳理。对于这个问题,我们一般关心收盘价随时间的变化趋势,这样我们就可以初步选定日期(DateNum)和收盘价(Pclose)作为重点研究对象。也就是说下一步,要对这这两个变量进行可视化。
. g$ ?3 v4 w* _3 V5 z, ?
% Y  |# e7 `' }" F/ |对于一个新手,我们还不知道如何绘图。但不要紧,新版 MATLAB 提供了更强大的绘图功能——“绘图”面板,这里提供了非常丰富的图形原型,如图 3 所示。) q( f8 v* L% ^3 \+ l* P2 V

+ a2 L$ G4 s% B& y& T
% P$ e2 l( s6 u! P" P% W; e# U4 M3 C2 `6 ^. }9 |- m2 x
                                                                                 图3 MATLAB绘图面板中的图例
+ D3 R% N9 ?6 ^: F% o# N* M
. M) t0 ~) I$ o1 i- q! T6 C8 X要注意,需要在工作区选中变量后绘图面板中的这些图标才会激活。接下来就可以选中一个中意的图标进行绘图,一般都直接先选第一个(plot)看一下效果,然后再浏览整个面板,看看有没有更合适的。下面我们进行绘图操作。8 P) O# x) g7 k8 m* u8 O+ R

1 ]6 D  v1 b5 cStep2.2:选中变量 DataNum 和 Pclose,在绘图面板中单机 plot 图标,马上可以得到这两个变量的可视化结果,如图 4 所示,同时还可以在命令窗口区看到绘制此图的命令:
$ s% U0 l1 a4 Z0 q5 z( u, _& N
9 p' c" l  l: A+ J& @8 J>> plot(DateNum,Pclose)2 R5 ~. k! q: ^1 }3 i
# M: N9 i8 x# [$ _
7 R- G3 M1 ?1 z6 u8 g
/ Z: ?! {* F: j8 K
                                                                                       图4 通过 plot 图标绘制的原图7 H) ?' I5 E) a. f

: q. |+ O: Z. N# s- y- a这样我们就知道了,下次再绘制这样的图直接用 plot 命令就可以了。一般情况下,用这种方式绘图的图往往不能满足我们的要求,比如我们希望更改:; q  o5 _; i! S2 w8 s4 r1 w

& s  e8 j  }; |4 t8 u7 r* z(1)曲线的颜色、线宽、形状;
% E/ e. M: _4 m3 ~1 z, O6 Y  Y. z$ W% g$ R) {( X9 h, d
(2)坐标轴的线宽、坐标,增加坐标轴描述;
0 L- T, r& L# a! ?( b% Q1 R+ N1 P1 R- J+ s3 o4 y+ y
(3)在同个坐标轴中绘制多条曲线。" q1 v0 J4 B5 E$ h

, z# k7 N* Z& z8 c$ W此时我们就需要了解更多关于命令 plot 的用法,这时就可以通过 MATLAB 强大的帮助系统来帮助我们实现期望的结果。最直接获取帮助的两个命令是 doc 和 help,对于新手来说,推荐使用 doc,因为 doc 直接打开的是帮助系统中的某个命令的用法说明,不仅全,而且有应用实例,这样就可以“照猫画虎”,直接参考实例,从而将实例快速转化成自己需要的代码。- C* q' m' q6 \0 K

9 B) N9 g/ [8 G: l$ q' q1 Y  C接下来我们就要考虑如何评估股票的价值和风险呢?* Z8 m' s  u0 k9 l( p
: x# V+ ~& g( U, ]' {( [
         对于一只好的股票,我们希望股票的增幅越大越好,体现在数学上,就是曲线的斜率越大越好。
$ A0 _/ u2 n7 ]: q4 g6 k7 y' S! q/ q% z- i: ^" l& P+ w
         对于风险,则可用最大回撤率来描述更合适,什么是最大回撤率?2 A$ y! o# k0 E; d' _7 I! g. `

. G, G  z2 Z6 `# [  b. E         最大回撤率的公式可以这样表达:- ]8 d: y, ^# V) r. Z: q7 @
+ U* Y  V+ q* N1 w" V
D为某一天的净值,i为某一天,j为i后的某一天,Di为第i天的产品净值,Dj则是Di后面某一天的净值0 O4 Q9 B7 n- |

+ U# Y4 V* Z" y7 q* Hdrawdown=max(Di-Dj)/Di,drawdown就是最大回撤率。其实就是对每一个净值进行回撤率求值,然后找出最大的。可以使用程序实现。最大回撤率越大,说明该股票的风险越高。所以最大回撤率越小,股票越好。# q$ S+ E8 }9 Y& Z0 y% D9 A  \+ v& H/ a
0 {5 Y7 ~3 u' d: }2 j. a$ p
           斜率和最大回撤率不妨一个一个来解决。我们先来看如何计算曲线的斜率。对于这个问题,比较简单,由于从数据的可视化结果来看,数据近似成线性,所以不妨用多项式拟合的方法来拟合该改组数据的方程,这样我们就可以得到斜率。
; M" \! n7 w1 X) x: Y5 |1 y
% R; q. T+ T) V2 }; F2 N/ bStep2.3:通过polyfit()多项式拟合的命令,并计算股票的价值,具体代码为:1 m8 |3 z# j- s" F' X/ E+ o, {: ^
6 I! \  ~* g! V0 U3 Q3 u# ?1 u
>> p = polyfit(DateNum,Pclose,1); % 多项式拟合
* g  r1 b% Q6 W0 F' g4 l8 `9 Z. s( I
>> value = p(1) % 将斜率赋值给value,作为股票的价值
5 v! U; o: i5 Q4 K! p' S8 ?8 M5 C0 `) w8 r- I
value =
3 T/ s# ?* z: _! p, z( H( P
7 Q& S% R% k; ~2 @. ^, G6 k5 |    0.12128 a/ M  X% }& n3 ^& V
0 K& |' t# m1 s# I
代码分析:%后面的内容是注释。polyfit()有三个参数,前两个大家都能明白是什么意思,那第三个参数是什么意思呢?它表示多项式的阶数,也就是最高次数。比如:在本例中,第三个参数为1,说明其为一次项,即一次函数。第三个参数为你要拟合的阶数,一阶直线拟合,二阶抛物线拟合,并非阶次越高越好,看拟合情况而定。polyfit()返回阶数为 n 的多项式 p(x) 的系数,p 中的系数按降幂排列。在本例中的P(1)指的是最高项的系数,即斜率。; j: l/ j( W6 N( Z5 g( p! d- M- t& U

( D, R( N9 r* K+ g* KStep2.4:用相似的方法,可以很快得到计算最大回撤的代码:  `. f; \& `" i7 N: Q
' d( {; a& W, v% r* x8 }
>> MaxDD = maxdrawdown(Pclose); % 计算最大回撤& j+ O" [2 \( q3 R

0 y7 ~+ H9 C5 ~+ `) ^' l>> risk = MaxDD  % 将最大回撤赋值给risk,作为股票的风险9 Q/ T$ r4 h! u' K, L7 {
+ n5 h" n# P3 S+ G( d9 B% W8 }
risk =/ \: k. m7 g( I: F! ^* i
' G1 ^( `* R. O6 Q; Q
    0.1155) H$ x3 o# b% e9 g
& t/ N. H2 v4 T* H) ]- t& r: W
代码分析:最大回撤率当然计算的是每天收盘时的股价。最大回撤率越大,说明该股票的风险越高。所以最大回撤率越小,股票越好。1 S; K' o# [& a. V" ?) h, x! @
. m  C7 M) t/ R
到此处,我们已经找到了评估股票价值和风险的方法,并能用 MALTAB 来实现了。但是,我们都是在命令行中实现的,并不能很方便地修改代码。而 MATLAB 最经典的一种用法就是脚本,因为脚本不仅能够完整地呈现整个问题的解决方法,同时更便于维护、完善、执行,优点很多。所以当我们的探索和开发工作比较成熟后,通常都会将这些有用的程序归纳整理起来,形成脚本。现在我们就来看如何快速开发解决该问题的脚本。
' d2 J9 J& B" h4 q2 m
1 Z2 A$ a. ]! ~4 g" uStep2.5:像 Step1.1 一样,重新选中数据文件,右键并单击“导入数据”菜单,待启动导入数据引擎后,选择“生成脚本”,然后就会得到导入数据的脚本,并保存该脚本。; a1 N7 O' y- }+ Q8 s3 R

% G: h# l4 O8 w脚本源代码中有些地方要注意:
4 A5 X, d: V' K
" d2 x5 P) ~3 U0 S2 S0 i8 T5 V       %%在matlab代码中的作用是将代码分块,上下两个%%之间的部分作为一块,在运行代码的时候可以分块运行,查看每一块代码的运行情况。常用于调试程序。%%相当于jupyter notebook中的cell。
  }$ R4 o) j% D/ p) ]5 d
6 U: A# |% h* b5 l' v+ y  ~, s' `       %后的内容是注释。( q0 a5 S& X# ^
, ^" B- B8 R: k7 e& [1 ^
        每句代码后面的分号作用为不在命令窗口显示执行结果。  D0 P/ @% t. j; d
) Q$ q. M4 m: G8 z" K
脚本源代码:2 \) \" o. c$ ?) p$ }

0 L2 R) @# f1 ~%% 预测股票的价值与风险- i* f4 X4 r1 u% R# q' Z# C# Z  g

  c4 Y' X! N, n% w% D! S. r  ?3 M" r%% 导入数据& l) L+ L) h6 r" s! Y
clc, clear, close all1 }5 p7 e+ {% R0 r
% clc:清除命令窗口的内容,对工作环境中的全部变量无任何影响 & G9 p% t: T) }8 Z' z! S
% clear:清除工作空间的所有变量
- D% X/ I) M8 d! p6 y% close all:关闭所有的Figure窗口
: B9 k5 F6 ~6 A& t) J: R" Q& S4 s
1 @  M- a1 V, T: k8 n% 导入数据8 i8 G, L8 L- H) v  e( ~
[~, ~, raw] = xlsread('sz000004.xlsx', 'Sheet1', 'A2:H7');
# s6 w% y3 M+ P* F  g% L% [num,txt,raw],~表示省略该部分的返回值0 C+ }8 v7 Z% Y2 |6 Q1 ^
% xlsread('filename','sheet', 'range'),第二个参数指数据在sheet1还是其他sheet部分,range表示单元格范围
  a* |1 T6 x2 M4 m" v2 ]' D3 p' p0 _% }7 y
% 创建输出变量7 I; z0 y2 `- V, D
data = reshape([raw{:}],size(raw));' L$ V+ c8 N3 l
% [raw{:}]指raw里的所有数据,size(raw):6 x 8 ,该语句把6x8的cell类型数据转换为6x8 double类型数据6 U+ g" Y$ j. g* p( n7 D8 v* Q
4 a* N. O5 _& c9 k! h7 r1 Q% d
% 将导入的数组分配列变量名称' [; C7 B' D( a0 |. C/ U
Date = data(:, 1); % 第一个参数表示从第一行到最后一行,第二个参数表示第一列
& l  k1 \4 j1 G3 X& W$ t: r( vDateNum = data(:, 2);4 G* F8 F! @+ a5 L6 N0 R
Popen = data(:, 3);0 m" g/ ?$ @" p- U- m
Phigh = data(:, 4);
1 M/ t% ]4 k, g# ?, `& i6 RPlow = data(:, 5);
$ f2 K, h" h' V, Q$ \1 K) RPclose = data(:, 6);  
& r% b- x+ i! T9 }Volum = data(:, 7); % Volume 表示股票成交量的意思,成交量=成交股数*成交价格 再加权求和
6 K5 k) |* \0 q6 b' KTurn = data(:, 8); % turn表示股票周转率,股票周转率越高,意味着该股股性越活泼,也就是投资人所谓的热门股, K( S6 ?6 ]: q# X3 _* d# \

! n! f- h  ~7 K& J9 N0 }% 清除临时变量data和raw
9 n4 X9 l2 h* }, _clearvars data raw;4 C" C4 F# [* s, I+ U

( e* O4 L! A6 h  C) M4 m+ A%% 数据探索. {# i; [& S/ k2 \$ ~
8 j1 y9 z8 s1 F5 N& M
figure % 创建一个新的图像窗口' p3 X4 b7 G, \" h
plot(DateNum, Pclose, 'k'); % 'k',曲线是黑色的,打印后不失真
& K# R5 z4 O0 F" M: L3 }datetick('x','mm-dd'); % 更改日期显示类型。参数x表示x轴,mm-dd表示月份和日。yyyy-mm-dd,如2018-10-27
; q$ N% j$ _  W' M! N% U! {, [% vxlabel('日期') % x轴9 [* C$ E: T  W# K3 u
ylabel('收盘价') % y轴
8 w- \5 J, @! D/ Z% V3 afigure( q& p$ \# e: s0 i% B9 }- ?
bar(Pclose) % 作为对照图形7 @4 q0 H! M$ j; a9 F
' R8 l9 B* x+ Q" T- u  |' }- n
%% 股票价值的评估
2 p% Z2 J6 O2 A! Y4 ?& d7 C
6 i" ~' {; N- D; D- f: w" Zp = polyfit(DateNum, Pclose, 1); % 多项式拟合
/ w* D9 ^7 Y( ?/ |, i/ q% polyfit()返回阶数为 n 的多项式 p(x) 的系数,p 中的系数按降幂排列+ g1 n" x. T+ L' e
P1 = polyval(p,DateNum); % 得到多项式模型的结果
, @7 E' J1 j. ]/ n- P1 Efigure/ F$ m. A2 p# {) X! i+ O6 J6 ~7 e
plot(DateNum,P1,DateNum,Pclose,'*g'); % 模型与原始数据的对照, '*g'表示绿色的*- V4 D0 I# L, U* m2 ~' g
value = p(1) % 将斜率赋值给value,作为股票的价值。p(1)最高项的次数6 Q# T% t$ ^, s6 E9 y9 z
. F- V& Y4 `! q5 N. d* d
%% 股票风险的评估( U6 H* n( a2 T+ ?$ j
MaxDD = maxdrawdown(Pclose); % 计算最大回撤
5 E; j# r7 ^0 P7 z5 urisk = MaxDD  % 将最大回撤赋值给risk,作为股票的风险8 s, D. M4 Q7 `5 R
  3、回归算法演练。8 w8 F$ w5 E+ |. L
8 V# |: x, n) h7 _8 x% s$ B
(1)一元线性回归
4 s8 D( ]$ ?2 J0 w. a7 L! l2 A( R4 a: }& @8 O) T
[ 例1 ] 近 10 年来,某市社会商品零售总额与职工工资总额(单位:亿元)的数据见表1,请建立社会商品零售总额与职工工资总额数据的回归模型。
  J7 c5 J/ ]3 w. H) F6 I! c1 s  t3 C; v7 M  S# E/ d. F7 m
/ b; ?/ @, q" [& K. P" M% Q
5 Z& ]% z. d: D' u  v
该问题是典型的一元回归问题,但先要确定是线性还是非线性,然后就可以利用对应的回归方法建立他们之间的回归模型了,具体实现的 MATLAB 代码如下:& y0 h4 a& \6 b3 [4 v- G: i9 N3 m

, M( L+ Z; x0 O$ K7 e  T(1)输入数据6 L% f4 r  a" H4 L  w* s

0 a4 G' \9 f$ C; s4 Z%% 输入数据
# O7 ?3 {0 _  dclc, clear, close all
+ ^6 A# [5 _/ c8 v% 职工工资总额
, Q/ h5 g, m9 }x = [23.8,27.6,31.6,32.4,33.7,34.90,43.2,52.8,63.8,73.4];( c: _- \/ Z) ?! Q5 Q0 j$ t% O; m
% 商品零售总额
% `5 M: [# i7 N1 `- _- ]y = [41.4,51.8,61.7,67.9,68.7,77.5,95.9,137.4,155.0,175.0];
! H4 ?, ]3 T# D(2)采用最小二乘回归
% @/ L2 ^" A. ?( R% D$ l: X/ Q: e3 a# N" |
%% 采用最小二乘法回归
1 z9 P5 a5 p/ k. V8 b% 作散点图
+ @* C  K; R( e+ _2 sfigure
$ U! k8 C5 n* e) Pplot(x,y,'r*') % 散点图,散点为红色& Q' y5 i. ^3 ~1 N2 D, x% S
xlabel('x(职工工资总额)','fontsize',12)' z. l: d% w! \2 t
ylabel('y(商品零售总额)','fontsize',12)! u6 E% v3 ]3 Z( i$ A3 Q
set(gca, 'linewidth',2) % 坐标轴线宽为25 g8 K" m1 J2 y$ E1 l2 M

' m1 a: o$ z3 I) l+ @0 u% 采用最小二乘法拟合
8 o) T# c1 {" J$ H, Q" TLxx = sum((x-mean(x)).^2); %在列表运算中,^与.^不同
' X. D* `; l4 ?7 W- L9 _: v) i$ RLxy = sum((x-mean(x)).*(y-mean(y)));3 ?8 c$ q$ V/ ^/ ?4 m" a( D  J
b1 = Lxy/Lxx;
! V$ m& c$ X1 |9 Ab0 = mean(y) - b1 * mean(x);
7 U" Y4 q6 ^1 f) L* h, iy1 = b1 * x + b0;
* m  a' f/ Z: C* i
; j8 h  n! j4 Y. u( |& E. dhold on % hold on是当前轴及图像保持而不被刷新,准备接受此后将绘制的图形,多图共存
0 j( H) B8 {  R, ?plot(x,y1, 'linewidth',2);. g9 B+ v) \3 i
运行本节程序,会得到如图5所示的回归图形。在用最小二乘回归之前,先绘制了数据的散点图,这样就可以从图形上判断这些数据是否近似成线性关系。当发现它们的确近似在一条线上后,再用线性回归的方法进行回归,这样也更符合我们分析数据的一般思路。
( Q& ]; X! ~; y  B3 A$ |7 m# b, W) O" o# F9 h8 Z& H8 {
4 W- T! b4 {: R2 B
' n7 }) ]6 M1 b6 p- v9 c1 ~, z: t
                                                                                                    图56 Y" [% T' c2 K2 @) V
/ A9 d& F2 a/ y" I" `# e
(3)采用 LinearModel.fit 函数进行线性回归
8 y% O3 j8 A/ w
/ s7 v* j  v+ C9 ?%% 采用 LinearModel.fit 函数进行线性回归; y6 T8 F# M! W4 [# ^  r1 U5 j
m2 = LinearModel.fit(x, y)
3 S) `! M2 i4 I7 l运行结果如下:
* Z( b% [! \# c% p2 J# H7 w! W) C$ q5 D2 P" N  r5 l/ C; R
m2 =
: G" l" C" o" m$ u
5 s1 W( ?7 z$ J- q; qLinear regression model:
" _  B# f; j+ G0 r5 T, f* B$ B. v) h9 D
    y ~ 1 + x1
3 l8 N+ E8 K; XEstimated Coefficients:" p. {! R7 }" N1 a) x' `5 D; c$ I7 `

" I, h" u  j+ t" |0 |: s$ X               Estimate      SE       tStat       pValue
8 z& O# d: _1 Z' F! Y& m
3 q: R, N" j$ V! ~    (Intercept)    -23.549      5.1028    -4.615     0.0017215
1 m  J6 Y8 l& K4 T5 ^& a0 c8 N4 z
    x1           2.7991     0.11456    24.435    8.4014e-09  O2 D; f0 c; I# Z8 @% ?& t- X% I

) z' J! ^9 _8 A7 v2 S0 PR-squared: 0.987,  Adjusted R-Squared 0.985" Y. W3 n' v2 b! M8 \
6 t% k; u: ~$ Z' l) a2 c% |
F-statistic vs. constant model: 597, p-value = 8.4e-09
4 O9 X) N( T- ]4 s, I# o5 G8 h7 y# m1 v
如下图,我们只需记住-23.594是一次函数的中x的系数,2.7991是一次函数中的常数项即可,其它的不用理会。$ r& v3 T/ p1 l4 [4 ]1 z
) @2 D# i5 e( B% m# v8 H

$ ?9 z0 W. x) B8 W8 d+ X* J9 V9 C/ _1 N  h8 c5 _% W/ b: _5 k
4)采用 regress 函数进行回归% Z( K) M* r9 B$ x/ p
+ h: }2 S4 n( g" d3 Z
%% 采用 regress 函数进行回归3 h; Y' `0 }$ _" E6 C: @% U
Y = y'7 h7 ~# }6 G- g2 A8 v+ E7 N8 z
X = [ones(size(x,2),1),x']
& _" O: J4 h% Q0 l5 H: i[b,bint,r,rint,s] = regress(Y,X)9 a$ @3 ^# T0 r- H" o6 g# ~; s2 z8 a! N
运行结果如下:0 t% m! O' |4 A8 \. Q/ A+ z

" C0 j% X8 l* B' F- Cb =7 k& v+ @7 y4 ^) j5 M8 t3 ]! t

( @1 b- o& [# D5 v% ]6 S  -23.5493( a. N2 V, s# k3 j2 K4 R# K

, x) g6 g" y" ~% Y! f    2.7991( o3 [% m' s$ F

. s% a) f& g* n我们只需记住-23.594是一次函数的中x的系数,2.7991是一次函数中的常数项即可,其它的不用理会。
6 [4 t; `4 L! ]8 r2 H% O7 o, r
2 b- y! m% o5 R# J: v* i. W/ c7 ](2)一元非线性回归
" {% w) i$ L$ Y% L1 E# o& y* I! p
8 W4 Y  Z* |* [. d# s+ b[ 例2 ] 为了解百货商店销售额 x 与流通率(这是反映商业活动的一个质量指标,指每元商品流转额所分摊的流通费用)y 之间的关系,收集了九个商店的有关数据(见表2)。请建立它们关系的数学模型。  u) Z& G/ G  F/ k" J9 U' s" x
9 t6 t+ B" s" i( X! e
; h8 ^+ c: p) J: A" D
/ G: o& d3 X5 [4 N2 |" f2 {
# N1 G0 \& p+ r" y; _! W9 a

/ r% S' h9 u* @" {9 A" N5 q& p# M        为了得到 x 与 y 之间的关系,先绘制出它们之间的散点图,如图 2 所示的“雪花”点图。由该图可以判断它们之间的关系近似为对数关系或指数关系,为此可以利用这两种函数形式进行非线性拟合,具体实现步骤及每个步骤的结果如下:
! Z% n+ a2 }, H- G3 O# S0 d% u
$ m& `4 l  ]0 J; u. Y(1)输入数据
2 K  v& y  |) q$ G* d: E5 X* A$ `% }1 }* }
%% 输入数据% ]5 [# Z9 z9 j& L. p5 _, i) z
clc, clear all, close all
- l5 ~' w% U/ v+ Y& I8 w1 Kx = [1.5, 4.5, 7.5,10.5,13.5,16.5,19.5,22.5,25.5];& Z& V3 D& Q6 x: X$ @1 k3 ^
y = [7.0,4.8,3.6,3.1,2.7,2.5,2.4,2.3,2.2];
" [; R3 z5 k: a, u8 W2 [plot(x, y, '*', 'linewidth', 1) % 这里的linewidth指的是散点大小* d- }) U: n' M
set(gca,'linewidth',2) % 设置坐标轴的线宽为2& W' ]* T$ a6 _( A+ o# ^7 G# Z$ B
xlabel('销售额x/万元','fontsize',12)
$ c4 ^5 r+ A. h% J7 W9 d$ Cylabel('流通率y/%','fontsize',12)- A  s1 x% Z8 _' i( u) Q
(2)对数形式非线性回归9 y% K' d5 Z) X
% e0 C3 y8 Z# y/ ~+ \
%% 对数形式非线性回归  B# ?% j" M3 z1 k1 W4 e" {
m1 = @(b,x) b(1) + b(2)*log(x);
/ _' N* b' r. ?5 Knonlinfit1 = fitnlm(x,y,m1,[0.01;0.01])% y" P' x* c: A
b = nonlinfit1.Coefficients.Estimate;! O$ ^6 B7 w8 N
Y1 = b(1,1) + b(2,1)*log(x);
' R) E8 [9 E& f( i: @2 ~, N+ ehold on
1 _6 y7 W+ y2 \  rplot(x, Y1, '--k', 'linewidth',2)$ p% ~8 M/ r7 T- w2 @8 [7 B) W( h
运行结果如下:7 ]: l* t+ b' K+ T: X' t' H

  F& [" c8 M% G/ \nonlinfit1 =
" {  U3 R# {( t2 R9 ]
2 X/ M: X) v/ u0 b0 U& D, B! XNonlinear regression model:
: x7 D' K' }7 K
+ V: w' p2 x1 @3 y# P0 C8 Q$ y    y ~ b1 + b2*log(x)
- H) Z; C1 a# s3 k' w; d$ m" ^) ~# V: q9 z7 L6 w: |, k% [) E
Estimated Coefficients:
+ J7 q6 \) Z4 |- Z& p4 E6 w" ~" R5 N, z! ^
          Estimate      SE        tStat       pValue
6 k- u8 J2 R8 b" s/ j
) N" _, G$ ?. F  Y    b1    7.3979      0.26667     27.742    2.0303e-08
+ ^9 s( L9 n- K( {' Z: M
. l/ T4 z2 y9 A5 O    b2    -1.713      0.10724    -15.974    9.1465e-07
& L/ c$ M% W, z! t5 K# y% w
5 a* \1 I) I$ ?& P; f, F# aR-Squared: 0.973,  Adjusted R-Squared 0.969
% }: l/ f# x1 L, |1 N# Z8 \- c' v
5 t# ~% ?# M8 l* `( e" d2 X* ?" Q, q3 ?8 YF-statistic vs. constant model: 255, p-value = 9.15e-07
2 R4 G4 p% e7 u( c  ^" C+ c5 ^( N3 T# I, P+ j4 h. S" S
(3)指数形式非线性回归6 z6 a1 F! f! l
, ~) a' x' `0 l! [
%% 指数形式非线性回归
  a6 D- I% p0 bm2 = 'y ~ b1*x^b2';
: f6 B" R/ z1 T7 Y* D7 ononlinfit2 = fitnlm(x,y,m2, [1;1])% U9 n$ t# A7 |  o
b1 = nonlinfit2.Coefficients.Estimate(1,1);
- C. g1 a* c5 P5 ]  ~5 Cb2 = nonlinfit2.Coefficients.Estimate(2,1)
, o, K- g0 B, U$ C* X7 j/ d: y9 |5 xY2 = b1*x.^b2;
; {  c; Z9 g9 y; o6 g3 s4 V5 Whold on;
1 `5 g+ N" m8 e& cplot(x,Y2,'r','linewidth',2)! W! x7 X+ Q. S
legend('原始数据','a+b*lnx','a*x^b') % 图例
6 R/ G- p% A( M( X2 `7 `" j* Q: ]运行结果如下:) s* z. H$ K; S

" `5 u" X7 ]0 K% enonlinfit2 =, v4 O+ l4 E, r$ w

) @" j; l: X! x  p' PNonlinear regression model:
& j( V8 t( H% @1 [# O0 k8 M/ u! j
    y ~ b1*x^b2& ^- h" j8 p7 O/ r( Q1 j: O) q

- B& u4 c2 }0 K& W7 I+ b$ {Estimated Coefficients:
) S3 {- u8 u7 L/ c
% M: u" v0 g# n$ |0 ^          Estimate       SE        tStat       pValue ' B* l$ h8 q2 H, }

5 u! W, p/ k9 b* E    b1      8.4112     0.19176     43.862    8.3606e-10
" K; ]  m( r- A6 @, R9 N( {7 H! p1 a) G' d) s
    b2    -0.41893    0.012382    -33.834    5.1061e-09) R" }% [& C9 k/ w0 y

' v1 o& k! Z5 E' C, ?7 i1 X9 bR-Squared: 0.993,  Adjusted R-Squared 0.992: H+ o9 l4 a9 M" I1 |9 D4 z0 N
6 W4 Y/ r% G) T% H' ~3 @
F-statistic vs. zero model: 3.05e+03, p-value = 5.1e-11
' f/ }( [4 ]: P; G) r8 I5 [+ ?- B. I4 @, r9 [
在该案例中,选择两种函数形式进行非线性回归,从回归结果来看,对数形式的决定系数为 0.973 ,而指数形式的为 0.993 ,优于前者,所以可以认为指数形式的函数形式更符合 y 与 x 之间的关系,这样就可以确定他们之间的函数关系形式了。
' i! B& l3 l5 A, J. |, C) M$ U5 D: i2 ?3 p3 i% X, C2 D$ q
2.多元回归8 K# G- L2 Q; c4 e" S0 ?2 u

. ?4 K% C. y3 z5 {5 r1.多元线性回归
& F, M2 J$ W- h: w
3 N0 G# D* x4 m; `( f# {, [6 M[ 例3 ] 某科学基金会希望估计从事某研究的学者的年薪 Y 与他们的研究成果(论文、著作等)的质量指标 X1、从事研究工作的时间 X2、能成功获得资助的指标 X3 之间的关系,为此按一定的实验设计方法调查了 24 位研究学者,得到如表3 所示的数据( i 为学者序号),试建立 Y 与 X1 , X2 , X3 之间关系的数学模型,并得出有关结论和作统计分析。$ _" r6 K) F! E# z! H2 B4 s

' h, T# T6 \" t3 n- Y
! i9 m, `) v9 h5 ~' |" _
5 w1 e. M; v* U( o1 m该问题是典型的多元回归问题,但能否应用多元线性回归,最好先通过数据可视化判断他们之间的变化趋势,如果近似满足线性关系,则可以执行利用多元线性回归方法对该问题进行回归。具体步骤如下:! [, g! Y, A' M4 v% @0 d7 c
* x2 ~( T, C9 D$ y
(1)作出因变量 Y 与各自变量的样本散点图6 Z% b: Z) v# l5 ^
& M" ~2 _- m3 d0 X; a
作散点图的目的主要是观察因变量 Y 与各自变量间是否有比较好的线性关系,以便选择恰当的数学模型形式。图3 分别为年薪 Y 与成果质量指标 X1、研究工作时间 X2、获得资助的指标 X3 之间的散点图。从图中可以看出这些点大致分布在一条直线旁边,因此,有比较好的线性关系,可以采用线性回归。绘制图3的代码如下:, w. l6 }5 v# z) U" z
' B7 E/ }  v8 M3 i: A, O
%% 作出因变量Y与各自变量的样本散点图7 ]1 c$ W& h& e5 S* ^
% x1,x2,x3,Y的数据/ N" @4 v$ i# |- u, S" r; }
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];* I& p$ P& Z1 o1 K/ M. l0 C
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];5 m9 j$ m/ M. a! u1 d& ^
x3=[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];- N' z3 {* w, R  [( f
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];, n% V7 N$ _9 J8 f# X1 E' C( Z
% 绘图,三幅图横向并排: Y7 j1 @- y' V" i& S& D# f  L
subplot(1,3,1),plot(x1,Y,'g*')
( W8 H- o6 |, \: osubplot(1,3,2),plot(x2,Y,'k+')
7 s7 w# U( M: O% j# Fsubplot(1,3,3),plot(x3,Y,'ro')  T7 z# f! @+ u; y" ]5 P
绘制的图形如下:
) f' O5 u6 p8 }$ c, C4 H+ m: y
7 Q, }- {3 w3 r8 \# n: T/ x4 F! w6 A# a4 U
/ c5 q6 g6 |& ~# u
(2)进行多元线性回归7 M) T8 ^, `1 i1 P7 k
: J5 x, _) u3 p' x
这里可以直接使用 regress 函数执行多元线性回归,注意以下代码模板,以后碰到多元线性问题直接套用代码,具体代码如下:
9 V: `- X6 Z. a0 x  `/ Q. b1 {  F! l- h( n9 @
%% 进行多元线性回归7 G, ]* }) i) o* Z9 v, T% f+ X
n = 24; m = 3; % 每个变量均有24个数据,共有3个变量
/ V, d6 [- B9 A* x, c& mX = [ones(n,1),x1',x2',x3'];' L( r1 w/ n" ^
[b,bint,r,rint,s]=regress(Y',X,0.05) % 0.05为预定显著水平,判断因变量y与自变量之间是否具有显著的线性相关关系需要用到。
: u' Z3 G! n6 D. {4 S: q运行结果如下:
5 ~0 s& p6 U2 f$ I7 R- ?6 F
; o/ R  c2 T! I  sb =3 v: V9 M# a7 Z( a

/ `0 X( z/ u$ c4 `, j9 X   18.0157
1 I; D, ^* m# U* F- P1 X    1.0817
4 s/ {, @3 e9 c! c  r2 A8 r0 N2 Z    0.3212/ r7 Y8 n  M6 I
    1.2835
' u) y" p* |5 q
# j1 q: B/ M9 D  |1 U7 a3 h1 b; E1 z" o
bint =- ^7 Y9 O$ w2 B7 m/ R9 D% i
. I1 c9 P& v, n  n: o
   13.9052   22.1262
& q0 Q' M9 Q5 u. }+ {# c    0.3900    1.7733
2 U8 v8 @2 `5 v7 E    0.2440    0.3984
1 K: g, {/ F5 o  @    0.6691    1.89791 x5 [, |* g% z" B
# _' E  x' @% a/ m5 h

, e, _; k$ p6 B& z3 \( A( P$ X8 @r =" ]5 C% N, `! }( P
& W1 f7 t* H7 Q* g- p& u& p, ~
    0.6781
1 r5 @- o) I6 u    1.91294 l8 ^- f; `: d$ O
   -0.1119
/ t) j/ K2 N/ P: |9 T( {) c    3.3114
1 s  P. Q  M/ S& s   -0.7424" C0 _3 i( Q$ ]& j
    1.2459+ a; b7 R8 a4 [0 f/ d
   -2.1022
! d3 k. a; W, @0 @    1.9650
& {, [7 a) ~/ Y   -0.3193; R2 a) |* r! R8 ^& Q9 [' k* r
    1.3466
5 ]% Z  O% d5 c8 h    0.8691
6 P; a0 i2 z# c   -3.2637- Q5 v( y! ?' W7 }* }# g
   -0.5115
( W. U( o, n6 ?+ p/ q% L   -1.1733
. n" Q3 t4 u5 s   -1.4910- G! e. P/ F+ ^* w% E
   -0.2972
" O0 Y- B8 k) S" v2 i5 W6 s/ f    0.1702* Q( S9 ]$ ~! N, Q$ d( T
    0.5799
/ L5 s0 Q% `/ C  K; H   -3.2856& ]6 o3 i& e/ V8 C0 w$ {
    1.1368
+ ]  @8 z% \( X9 H   -0.8864
# m; V  s: _0 A; x: b* k   -1.4646
* k$ ?9 r4 W: d9 }! e    0.8032
5 O3 @4 |6 N# t4 v    1.63015 j7 k2 }& }. Z! o/ n/ N& m3 x: Q
: E( v4 b- r  [" t0 ]; F0 N( W5 u9 [

6 q3 P7 d5 Z! frint =* H7 r# A0 K2 O
9 a" I1 `; b! Q% |$ u2 G+ I/ f
   -2.7017    4.0580
3 R" o; Q/ c6 {, D: t. t! I1 }4 i   -1.6203    5.4461
# ?% ^& }$ F! z% T' L! I# u7 ]5 |   -3.6190    3.3951
& D9 h' }5 L  i$ K2 W. Q    0.0498    6.57295 M  ~% n9 g) U: h" w8 k
   -4.0560    2.5712* q2 \7 p# s8 n3 J( P$ Q* Z! O
   -2.1800    4.6717  _5 m& e  ?' d5 t: F
   -5.4947    1.2902
6 P' _5 M( X, t; \3 h   -1.3231    5.2531
% @' G, [# d4 Q. X   -3.5894    2.9507
/ O6 r/ S0 F! E3 W   -1.7678    4.4609' t- o8 ]$ r- |- m
   -2.7146    4.4529
  N7 N, Q, r  m$ ^3 [# i3 t   -6.4090   -0.11833 b2 x3 w6 G! R% F4 A- Q
   -3.6088    2.5859
6 l" `9 v0 O5 i; t+ ]& q6 M% ?5 F- }   -4.7040    2.3575! L3 k: J; e; t2 e
   -4.8249    1.84293 m+ p! N. o( i  G' y, A
   -3.7129    3.11859 _7 w, ?) W1 x, `: a
   -3.0504    3.3907
; N# {; S8 a& G   -2.8855    4.0453
) J2 X+ ]" m! T$ R; q, d) o8 a   -6.2644   -0.3067; X; C9 e9 d+ j, T5 h/ P8 x
   -2.1893    4.4630
+ m4 I$ q4 V0 J: r% E) X! b   -4.4002    2.6273
, a( L) O6 o6 K7 ]3 |   -4.8991    1.9699
$ h  X; z7 J: y7 Z* v- \   -2.4872    4.0937! c) q$ @. e2 _( W! t) }
   -1.8351    5.0954
! `$ K# [2 D- j% `( J
6 \- N$ o; [" M  Z8 d
$ c, |" p6 p* D: o5 d, ss =" ?! }+ T: j' D

' t5 {- x7 h/ v7 ~2 {    0.9106   67.9195    0.0000    3.0719
9 _. I/ B0 e4 I. T: \看到如此长的运行结果,我们不要害怕,因为里面很多数据是没用的,我们只需提取有用的数据。
' t5 Y/ d" q9 T2 Z  I* J3 K# [7 K1 H& w( y/ h  ]! {
在运行结果中,很多数据我们不需理会,我们真正需要用到的数据如下:
/ p: [5 c. z% [; G$ M- q6 @% [3 P" W$ t9 u% F4 `% X! E
b =
. Z4 `! U$ ~1 f
& I# \8 x6 b& I8 Y7 G( D* \   18.0157
$ i* S* Y4 [6 _8 b* t    1.08171 ^. m  x/ C  \4 [% Y  |
    0.3212
& b8 ]0 m, x. Z1 S    1.28353 w) z8 n, \) q" H! T5 O  C
6 o/ {0 B, V4 D
s =
) E* t# Y( B) H: A- @# h' T; a) b* g- n% E) ~$ X
    0.9106   67.9195    0.0000    3.0719, N8 F, }' \1 d! h) f2 Y$ f: V
回归系数 b = (β0,β1,β2,β3) = (18.0157, 1.0817, 0.3212, 1.2835),回归系数的置信区间,以及统计变量 stats(它包含四个检验统计量:相关系数的平方R^2,假设检验统计量 F,与 F 对应的概率 p,s^2 的值)。观察表4的数据,会发现它来源于运行结果中的b和s:2 p8 ]+ s6 }9 x

3 v. ?5 }! m4 R% b/ W7 w' V6 Z* Y$ P5 h; T4 P: \. A8 n# z6 K2 J

8 m. t- [: ]4 }, a根据β0,β1,β2,β3,我们初步得出回归方程为:
- A' f8 O4 k7 j! |5 e
% G0 T, }( S1 n# ^6 U0 K- [5 T; Y. b3 Q" `) f% J1 y; R2 d
# U+ c# A5 v2 l
如何判断该回归方程是否符合该模型呢?有以下3种方法:
+ x  J5 \' f* W( S' F2 I& v; Y. H, Z; S8 m
1)相关系数 R 的评价:本例 R 的绝对值为 0.9542 ,表明线性相关性较强。1 M9 W0 b: ~: h0 v% {
9 p. @$ W: y. C8 {6 W
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。* c$ l' m* u. C! Y) z" p5 f

! R, s& b* x; Z3)p 值检验:若 p < α(α 为预定显著水平),则说明因变量 y 与自变量 x1,x2,...,xm之间显著地有线性相关关系。本例输出结果,p<0.0001,显然满足 p<α=0.05。  Y) W$ B) v! U, b5 d: S# L
& Z: U7 k: F2 I/ [) C
以上三种统计推断方法推断的结果是一致的,说明因变量 y 与自变量之间显著地有线性相关关系,所得线性回归模型可用。s^2 当然越小越好,这主要在模型改进时作为参考。0 ?9 V5 F: W% N6 S! p2 e

4 ?9 R% ^4 I9 f9 ~; U3. 逐步回归- U% y; c2 S% D  S% K. u/ n
( B6 `  l  R6 P2 T
[ 例4 ] (Hald,1960)Hald 数据是关于水泥生产的数据。某种水泥在凝固时放出的热量 Y(单位:卡/克)与水泥中 4 种化学成品所占的百分比有关:' L- _& M. M( h! T- n

: b5 Y. n* y: v/ K" f3 N" r( b. z4 t6 N; Z* Z
( d2 \9 h" r! J
在生产中测得 12 组数据,见表5,试建立 Y 关于这些因子的“最优”回归方程。- O& {$ s) v7 u8 k' @+ v

% p9 M) I3 L3 Q5 d: X8 F  z9 L% h

- _; ?. X1 G3 b  X9 m9 C) W3 ?3 a对于例 4 中的问题,可以使用多元线性回归、多元多项式回归,但也可以考虑使用逐步回归。从逐步回归的原理来看,逐步回归是以上两种回归方法的结合,可以自动使得方程的因子设置最合理。对于该问题,逐步回归的代码如下:
" W4 u; i  u* m3 q8 X" T% Y4 m
! W4 Y2 d$ |, S2 d, r%% 逐步回归/ `' Q' B; B3 {3 g, M7 ^
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];   %自变量数据* U! e7 p) X& _0 z9 S0 l
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];  %因变量数据7 t: W4 m8 {- |* i5 C" e) F- Y0 V
stepwise(X,Y,[1,2,3,4],0.05,0.10)% in=[1,2,3,4]表示X1、X2、X3、X4均保留在模型中6 Z. r9 F7 a/ H( s9 R
程序执行后得到下列逐步回归的窗口,如图 4 所示。
* P' W' A/ |1 O% _& t" D
1 V9 _2 d- i6 R& A. s, s5 [0 W& C9 U' f# K

: x# W( ~. u; g( @                                                                                                             图4+ i7 J  l& Q# j" d0 a

0 r, \& w5 S& }( a0 `1 t在图 4 中,用蓝色行显示变量 X1、X2、X3、X4 均保留在模型中,窗口的右侧按钮上方提示:将变量X4剔除回归方程(Move X4 out),单击 Next Step 按钮,即进行下一步运算,将第 4 列数据对应的变量 X4 剔除回归方程。单击 Next Step 按钮后,剔除的变量 X3 所对应的行用红色表示,同时又得到提示:将变量 X3 剔除回归方程(Move X3 out),单击 Next Step 按钮,这样一直重复操作,直到 “Next Step” 按钮变灰,表明逐步回归结束,此时得到的模型即为逐步回归最终的结果。最终结果如下:
( `! C) Y3 r, J; F; L/ z
* B5 J0 i! e. S4 G, ]4 O- O' J% o5 k4 v7 L9 H- P
9 X: `' V7 J6 }: Y
4. 逻辑回归1 i& L$ E) c2 l9 W( O" ?
" V; Q, t  ^3 \, D8 _
[ 例5 ] 企业到金融商业机构贷款,金融商业机构需要对企业进行评估。评估结果为 0 , 1 两种形式,0 表示企业两年后破产,将拒绝贷款,而 1 表示企业 2 年后具备还款能力,可以贷款。在表 6 中,已知前 20 家企业的三项评价指标值和评估结果,试建立模型对其他 5 家企业(企业 21-25)进行评估。
+ j& ?7 E% ]2 _: U9 y+ B# n2 G) I3 [( w- G  W& ]
# C3 o, i2 w8 D7 _4 s

6 S  `( B3 H# g对于该问题,很明显可以用 Logistic 模型来回归,具体求解程序如下:. A7 ?% W: A$ R1 _2 i( s
  b) H- r  d6 X  N# l* O* {
程序中需要用到的数据文件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
  B$ J9 t/ {0 ]* @; l! u  I+ N- t" D7 O+ b  _3 {; G; O0 K2 z
% logistic回归) O! L# I/ P" k- W& x) u
: {8 c' V" Y; z/ h" p; Y
%% 导入数据
' `4 C5 Q* I* qclc,clear,close all
# m+ g: W2 _" R# E6 B! }: D8 F7 eX0 = xlsread('logistic_ex1.xlsx','A2:C21'); % 前20家企业的三项评价指标值,即回归模型的输入" x) I; }7 \/ |" r
Y0 = xlsread('logistic_ex1.xlsx','D221'); % 前20家企业的评估结果,即回归模型的输出
  Z9 [* J9 H; A2 o  GX1 = xlsread('logistic_ex1.xlsx','A2:C26'); % 预测数据输入
8 j1 k% @! c9 Y7 R8 w- T
9 L# E2 w% r  ^5 P5 W+ S" ^- X%% 逻辑函数% ?! O) X8 A& V. N- u; W
GM = fitglm(X0,Y0,'Distribution','binomial');
, S7 j& I/ e* b; u. FY1 = predict(GM,X1);6 n& [# E/ \6 Q" m( P- @. @, g
* H' [: k- |' W2 c; Q2 p
%% 模型的评估  |$ ~4 j; ]3 K) ^
N0 = 1:size(Y0,1); % N0 = [1,2,3,4,……,20]
" H2 _) d! ^* vN1 = 1:size(Y1,1); % N1 = [1,2,3,4,……,25]  e7 [4 t9 y( J
plot(N0',Y0,'-kd'); % N0'指的是对N0'进行转置,N0'和Y0的形式相同,该行代码绘制的是前20家企业的评估结果
! r! `3 @2 n% ]4 j) S1 m% plot()中的参数'-kd'的解析:-代表直线,k代表黑色,d代表菱形符号/ S- {9 ?) Z5 g% W/ K- `3 l5 c
hold on;) n+ U- B5 t9 E
scatter(N1',Y1,'b'); % N1'指的是对N1'进行转置,N1'和Y1的形式相同
0 ~* q3 b1 \) ^4 o1 b+ N1 |# c2 q7 Dxlabel('企业编号');2 [5 b6 y' Z! ~$ j
ylabel('输出值');
- b4 }; G8 H, E% u. N4 p( x' @+ ^得到的回归结果与原始数据的比较如图5所示。) r3 [- K- E8 y% u# e8 p! W
+ ~% Y/ n; ^- @; L& m# y

' x0 a2 R; A. I
- G  G' ?9 H+ v2 ^- q                                                                   图5
$ I1 S' o1 R  I/ K( m( K
3 C  a+ @* a* E0 I三、总结与感悟。
8 Y1 Y. a: M# E6 f. `- P7 a+ i8 V, w6 C/ O) z5 P
        总结:通过这次学习,我了解到Matlab在数学建模竞赛中使用广泛;在评估股票价值与风险的小实例中,我掌握了用Matlab去建模的基本方法和步骤;在回归算法的学习过程中,我掌握了一元线性回归、一元非线性回归、多元线性回归、逐步回归、逻辑回归的算法。& ?- V4 U# _4 U

# |4 c+ b9 ^: g6 J# `$ J2 F        感悟:正确且高效的 MATLAB 编程理念就是以问题为中心的主动编程。我们传统学习编程的方法是学习变量类型、语法结构、算法以及编程的其他知识,因为学习时候是没有目标的,也不知道学的知识什么时候能用到,收效甚微。而以问题为中心的主动编程,则是先找到问题的解决步骤,然后在 MATLAB 中一步一步地去实现。在每步实现的过程中,遇到问题,查找知识(互联网时代查询知识还是很容易的),定位方法,再根据方法,查询 MATLAB 中的对应函数,学习函数用法,回到程序,解决问题。在这个过程中,知识的获取都是为了解决问题的,也就是说每次学习的目标都是非常明确的,学完之后的应用就会强化对知识的理解和掌握,这样即学即用的学习方式是效率最高,也是最有效的方式。最重要的是,这种主动的编程方式会让学习者体验到学习的成就感的乐趣,有成就感,自然就强化对编程的自信了。这种内心的自信和强大在建模中会发挥意想不到的力量,所为信念的力量。% f7 s6 ]; c+ f

4 R, [  c, ]- B& X+ }( U) F$ C5 S/ C6 X





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