Matlab与数学建模

一、学习目标。

（1）了解Matlab与数学建模竞赛的关系。

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

（3）掌握Matlab数学建模的回归算法。

二、实例演练。

1、谈谈你对Matlab与数学建模竞赛的了解。

Matlab在数学建模中使用广泛：MATLAB 是公认的最优秀的数学模型求解工具，在数学建模竞赛中超过 95% 的参赛队使用 MATLAB 作为求解工具，在国家奖队伍中，MATLAB 的使用率几乎 100%。虽然比较知名的数模软件不只 MATLAB。

人们喜欢使用Matlab去数学建模的原因：

（1）MATLAB 的数学函数全，包含人类社会的绝大多数数学知识。

（2）MATLAB 足够灵活，可以按照问题的需要，自主开发程序，解决问题。

（3）MATLAB易上手，本身很简单，不存在壁垒。掌握正确的 MATLAB 使用方法和实用的小技巧，在半小时内就可以很快地变成 MATLAB 高手了。

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

数学建模竞赛中的 MATLAB 水平要求：

要想在全国大学生数学建模竞赛中拿到国奖， MATLAB 技能是必备的。 具体的技能水平应达到：

1）了解 MATLAB 的基本用法，包括几个常用的命令，如何获取帮助，脚本结构，程序的分节与注释，矩阵的基本操作，快捷绘图方式；

2）熟悉 MATLAB 的程序结构，编程模式，能自由地创建和引用函数（包括匿名函数）；

3）熟悉常见模型的求解算法和套路，包括连续模型，规划模型，数据建模类的模型；

4）能够用 MALTAB 程序将机理建模的过程模拟出来，就是能够建立和求解没有套路的数学模型。

要想达到如上要求， 不能按照传统的学习方式一步一步地学习， 而要结合上述提到的学习理念制定科学的训练计划。

2、已知股票的交易数据：日期、开盘价、最高价、最低价、收盘价、成交量和换手率，试用某种方法来评价这只股票的价值和风险。如何用MATLAB去求解该问题？（交易数据：点击此处获取数据）

解题步骤：

第一阶段：从外部读取数据

Step1.1：把数据文件sz000004.xls拖曳进‘当前文件夹区’，选中数据文件sz000004.xls，右键，将弹出右键列表，很快可发现有个“导入数据”菜单，如图 1 所示。

图1. 启动导入数据引擎示意图

Step1.2：单击“导入数据”这个按钮，则很快发现起到一个导入数据引擎，如图 4 所示。

图2. 导入数据界面

Step1.3：观察图 2，在右上角有个“导入所选内容”按钮，则可直接单击之。马上我们就会发现在 MATLAB 的工作区（当前内存中的变量）就会显示这些导入的数据，并以列向量的方式表示，因为默认的数据类型就是“列向量”，当然您可以可以选择其他的数据类型，大家不妨做几个实验，观察一下选择不同的数据类型后会结果会有什么不同。至此，第一步获取数据的工作的完成。

 

第二阶段：数据探索和建模

现在重新回到问题，对于该问题，我们的目标是能够评估股票的价值和风险，但现在我们还不知道该如何去评估，MATLAB 是工具，不能代替我们决策用何种方法来评估，但是可以辅助我们得到合适的方法，这就是数据探索部分的工作。下面我们就来尝试如何在 MATLAB 中进行数据的探索和建模。

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

由于变量比较多，所以还有必要对这些变量进行初步的梳理。对于这个问题，我们一般关心收盘价随时间的变化趋势，这样我们就可以初步选定日期（DateNum）和收盘价（Pclose）作为重点研究对象。也就是说下一步，要对这这两个变量进行可视化。

对于一个新手，我们还不知道如何绘图。但不要紧，新版 MATLAB 提供了更强大的绘图功能——“绘图”面板，这里提供了非常丰富的图形原型，如图 3 所示。

图3 MATLAB绘图面板中的图例

要注意，需要在工作区选中变量后绘图面板中的这些图标才会激活。接下来就可以选中一个中意的图标进行绘图，一般都直接先选第一个（plot）看一下效果，然后再浏览整个面板，看看有没有更合适的。下面我们进行绘图操作。

Step2.2：选中变量 DataNum 和 Pclose，在绘图面板中单机 plot 图标，马上可以得到这两个变量的可视化结果，如图 4 所示，同时还可以在命令窗口区看到绘制此图的命令：

>> plot(DateNum,Pclose)

图4 通过 plot 图标绘制的原图

这样我们就知道了，下次再绘制这样的图直接用 plot 命令就可以了。一般情况下，用这种方式绘图的图往往不能满足我们的要求，比如我们希望更改：

（1）曲线的颜色、线宽、形状；

（2）坐标轴的线宽、坐标，增加坐标轴描述；

（3）在同个坐标轴中绘制多条曲线。

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

接下来我们就要考虑如何评估股票的价值和风险呢？

对于一只好的股票，我们希望股票的增幅越大越好，体现在数学上，就是曲线的斜率越大越好。

对于风险，则可用最大回撤率来描述更合适，什么是最大回撤率？

最大回撤率的公式可以这样表达:

D为某一天的净值，i为某一天，j为i后的某一天，Di为第i天的产品净值，Dj则是Di后面某一天的净值

drawdown=max(Di-Dj)/Di，drawdown就是最大回撤率。其实就是对每一个净值进行回撤率求值，然后找出最大的。可以使用程序实现。最大回撤率越大，说明该股票的风险越高。所以最大回撤率越小，股票越好。

斜率和最大回撤率不妨一个一个来解决。我们先来看如何计算曲线的斜率。对于这个问题，比较简单，由于从数据的可视化结果来看，数据近似成线性，所以不妨用多项式拟合的方法来拟合该改组数据的方程，这样我们就可以得到斜率。

Step2.3：通过polyfit()多项式拟合的命令，并计算股票的价值，具体代码为：

>> p = polyfit(DateNum,Pclose,1); % 多项式拟合

>> value = p(1) % 将斜率赋值给value，作为股票的价值

value =

0.1212

代码分析：%后面的内容是注释。polyfit()有三个参数，前两个大家都能明白是什么意思，那第三个参数是什么意思呢？它表示多项式的阶数，也就是最高次数。比如：在本例中，第三个参数为1，说明其为一次项，即一次函数。第三个参数为你要拟合的阶数，一阶直线拟合，二阶抛物线拟合，并非阶次越高越好，看拟合情况而定。polyfit()返回阶数为 n 的多项式 p(x) 的系数，p 中的系数按降幂排列。在本例中的P(1)指的是最高项的系数，即斜率。

Step2.4：用相似的方法，可以很快得到计算最大回撤的代码：

>> MaxDD = maxdrawdown(Pclose); % 计算最大回撤

>> risk = MaxDD % 将最大回撤赋值给risk，作为股票的风险

risk =

0.1155

代码分析：最大回撤率当然计算的是每天收盘时的股价。最大回撤率越大，说明该股票的风险越高。所以最大回撤率越小，股票越好。

到此处，我们已经找到了评估股票价值和风险的方法，并能用 MALTAB 来实现了。但是，我们都是在命令行中实现的，并不能很方便地修改代码。而 MATLAB 最经典的一种用法就是脚本，因为脚本不仅能够完整地呈现整个问题的解决方法，同时更便于维护、完善、执行，优点很多。所以当我们的探索和开发工作比较成熟后，通常都会将这些有用的程序归纳整理起来，形成脚本。现在我们就来看如何快速开发解决该问题的脚本。

Step2.5：像 Step1.1 一样，重新选中数据文件，右键并单击“导入数据”菜单，待启动导入数据引擎后，选择“生成脚本”，然后就会得到导入数据的脚本，并保存该脚本。

脚本源代码中有些地方要注意：

%%在matlab代码中的作用是将代码分块，上下两个%%之间的部分作为一块，在运行代码的时候可以分块运行，查看每一块代码的运行情况。常用于调试程序。%%相当于jupyter notebook中的cell。

%后的内容是注释。

每句代码后面的分号作用为不在命令窗口显示执行结果。

脚本源代码：

%% 预测股票的价值与风险

%% 导入数据
clc, clear, close all
% clc：清除命令窗口的内容，对工作环境中的全部变量无任何影响 
% clear：清除工作空间的所有变量 
% close all:关闭所有的Figure窗口

% 导入数据
[~, ~, raw] = xlsread('sz000004.xlsx', 'Sheet1', 'A2:H7');
% [num,txt,raw],~表示省略该部分的返回值
% xlsread('filename'，'sheet'， 'range')，第二个参数指数据在sheet1还是其他sheet部分，range表示单元格范围

% 创建输出变量
data = reshape([raw{:}],size(raw));
% [raw{:}]指raw里的所有数据，size(raw)：6 x 8 ，该语句把6x8的cell类型数据转换为6x8 double类型数据

% 将导入的数组分配列变量名称
Date = data(:, 1); % 第一个参数表示从第一行到最后一行，第二个参数表示第一列
DateNum = data(:, 2);
Popen = data(:, 3);
Phigh = data(:, 4);
Plow = data(:, 5);
Pclose = data(:, 6);  
Volum = data(:, 7); % Volume 表示股票成交量的意思，成交量=成交股数*成交价格 再加权求和
Turn = data(:, 8); % turn表示股票周转率，股票周转率越高，意味着该股股性越活泼，也就是投资人所谓的热门股

% 清除临时变量data和raw
clearvars data raw;

%% 数据探索

figure % 创建一个新的图像窗口
plot(DateNum, Pclose, 'k'); % 'k',曲线是黑色的，打印后不失真
datetick('x','mm-dd'); % 更改日期显示类型。参数x表示x轴，mm-dd表示月份和日。yyyy-mm-dd,如2018-10-27
xlabel('日期') % x轴
ylabel('收盘价') % y轴
figure
bar(Pclose) % 作为对照图形

%% 股票价值的评估

p = polyfit(DateNum, Pclose, 1); % 多项式拟合
% polyfit()返回阶数为 n 的多项式 p(x) 的系数，p 中的系数按降幂排列
P1 = polyval(p,DateNum); % 得到多项式模型的结果
figure
plot(DateNum,P1,DateNum,Pclose,'*g'); % 模型与原始数据的对照, '*g'表示绿色的*
value = p(1) % 将斜率赋值给value，作为股票的价值。p(1)最高项的次数

%% 股票风险的评估
MaxDD = maxdrawdown(Pclose); % 计算最大回撤
risk = MaxDD  % 将最大回撤赋值给risk，作为股票的风险

3、回归算法演练。

（1）一元线性回归

[ 例1 ] 近 10 年来，某市社会商品零售总额与职工工资总额（单位：亿元）的数据见表1，请建立社会商品零售总额与职工工资总额数据的回归模型。

该问题是典型的一元回归问题，但先要确定是线性还是非线性，然后就可以利用对应的回归方法建立他们之间的回归模型了，具体实现的 MATLAB 代码如下：

（1）输入数据

%% 输入数据
clc, clear, close all
% 职工工资总额
x = [23.8,27.6,31.6,32.4,33.7,34.90,43.2,52.8,63.8,73.4];
% 商品零售总额
y = [41.4,51.8,61.7,67.9,68.7,77.5,95.9,137.4,155.0,175.0];

（2）采用最小二乘回归

%% 采用最小二乘法回归
% 作散点图
figure
plot(x,y,'r*') % 散点图，散点为红色
xlabel('x(职工工资总额)','fontsize',12)
ylabel('y(商品零售总额)','fontsize',12)
set(gca, 'linewidth',2) % 坐标轴线宽为2

% 采用最小二乘法拟合
Lxx = sum((x-mean(x)).^2); %在列表运算中，^与.^不同
Lxy = sum((x-mean(x)).*(y-mean(y)));
b1 = Lxy/Lxx;
b0 = mean(y) - b1 * mean(x);
y1 = b1 * x + b0;

hold on % hold on是当前轴及图像保持而不被刷新，准备接受此后将绘制的图形，多图共存
plot(x,y1, 'linewidth',2);

运行本节程序，会得到如图5所示的回归图形。在用最小二乘回归之前，先绘制了数据的散点图，这样就可以从图形上判断这些数据是否近似成线性关系。当发现它们的确近似在一条线上后，再用线性回归的方法进行回归，这样也更符合我们分析数据的一般思路。

图5

（3）采用 LinearModel.fit 函数进行线性回归

%% 采用 LinearModel.fit 函数进行线性回归
m2 = LinearModel.fit(x, y)

运行结果如下：

m2 =

Linear regression model:

y ~ 1 + x1
Estimated Coefficients:

Estimate SE tStat pValue

(Intercept) -23.549 5.1028 -4.615 0.0017215

x1 2.7991 0.11456 24.435 8.4014e-09

R-squared: 0.987, Adjusted R-Squared 0.985

F-statistic vs. constant model: 597, p-value = 8.4e-09

如下图，我们只需记住-23.594是一次函数的中x的系数，2.7991是一次函数中的常数项即可，其它的不用理会。

4）采用 regress 函数进行回归

%% 采用 regress 函数进行回归
Y = y'
X = [ones(size(x,2),1),x']
[b,bint,r,rint,s] = regress(Y,X)

运行结果如下：

b =

-23.5493

2.7991

我们只需记住-23.594是一次函数的中x的系数，2.7991是一次函数中的常数项即可，其它的不用理会。

（2）一元非线性回归

[ 例2 ] 为了解百货商店销售额 x 与流通率（这是反映商业活动的一个质量指标，指每元商品流转额所分摊的流通费用）y 之间的关系，收集了九个商店的有关数据（见表2）。请建立它们关系的数学模型。

为了得到 x 与 y 之间的关系，先绘制出它们之间的散点图，如图 2 所示的“雪花”点图。由该图可以判断它们之间的关系近似为对数关系或指数关系，为此可以利用这两种函数形式进行非线性拟合，具体实现步骤及每个步骤的结果如下：

（1）输入数据

%% 输入数据
clc, clear all, close all
x = [1.5, 4.5, 7.5,10.5,13.5,16.5,19.5,22.5,25.5];
y = [7.0,4.8,3.6,3.1,2.7,2.5,2.4,2.3,2.2];
plot(x, y, '*', 'linewidth', 1) % 这里的linewidth指的是散点大小
set(gca,'linewidth',2) % 设置坐标轴的线宽为2
xlabel('销售额x/万元','fontsize',12)
ylabel('流通率y/%','fontsize',12)

（2）对数形式非线性回归

%% 对数形式非线性回归
m1 = @(b,x) b(1) + b(2)*log(x);
nonlinfit1 = fitnlm(x,y,m1,[0.01;0.01])
b = nonlinfit1.Coefficients.Estimate;
Y1 = b(1,1) + b(2,1)*log(x);
hold on 
plot(x, Y1, '--k', 'linewidth',2)

运行结果如下：

nonlinfit1 =

Nonlinear regression model:

y ~ b1 + b2*log(x)

Estimated Coefficients:

Estimate SE tStat pValue

b1 7.3979 0.26667 27.742 2.0303e-08

b2 -1.713 0.10724 -15.974 9.1465e-07

R-Squared: 0.973, Adjusted R-Squared 0.969

F-statistic vs. constant model: 255, p-value = 9.15e-07

（3）指数形式非线性回归

%% 指数形式非线性回归
m2 = 'y ~ b1*x^b2';
nonlinfit2 = fitnlm(x,y,m2, [1;1])
b1 = nonlinfit2.Coefficients.Estimate(1,1);
b2 = nonlinfit2.Coefficients.Estimate(2,1)
Y2 = b1*x.^b2;
hold on;
plot(x,Y2,'r','linewidth',2)
legend('原始数据','a+b*lnx','a*x^b') % 图例

运行结果如下：

nonlinfit2 =

Nonlinear regression model:

y ~ b1*x^b2

Estimated Coefficients:

Estimate SE tStat pValue

b1 8.4112 0.19176 43.862 8.3606e-10

b2 -0.41893 0.012382 -33.834 5.1061e-09

R-Squared: 0.993, Adjusted R-Squared 0.992

F-statistic vs. zero model: 3.05e+03, p-value = 5.1e-11

在该案例中，选择两种函数形式进行非线性回归，从回归结果来看，对数形式的决定系数为 0.973 ，而指数形式的为 0.993 ，优于前者，所以可以认为指数形式的函数形式更符合 y 与 x 之间的关系，这样就可以确定他们之间的函数关系形式了。

2.多元回归

1.多元线性回归

[ 例3 ] 某科学基金会希望估计从事某研究的学者的年薪 Y 与他们的研究成果（论文、著作等）的质量指标 X1、从事研究工作的时间 X2、能成功获得资助的指标 X3 之间的关系，为此按一定的实验设计方法调查了 24 位研究学者，得到如表3 所示的数据（ i 为学者序号），试建立 Y 与 X1 , X2 , X3 之间关系的数学模型，并得出有关结论和作统计分析。

该问题是典型的多元回归问题，但能否应用多元线性回归，最好先通过数据可视化判断他们之间的变化趋势，如果近似满足线性关系，则可以执行利用多元线性回归方法对该问题进行回归。具体步骤如下：

（1）作出因变量 Y 与各自变量的样本散点图

作散点图的目的主要是观察因变量 Y 与各自变量间是否有比较好的线性关系，以便选择恰当的数学模型形式。图3 分别为年薪 Y 与成果质量指标 X1、研究工作时间 X2、获得资助的指标 X3 之间的散点图。从图中可以看出这些点大致分布在一条直线旁边，因此，有比较好的线性关系，可以采用线性回归。绘制图3的代码如下：

%% 作出因变量Y与各自变量的样本散点图
% x1,x2,x3,Y的数据
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];
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];
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];
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];
% 绘图，三幅图横向并排
subplot(1,3,1),plot(x1,Y,'g*')
subplot(1,3,2),plot(x2,Y,'k+')
subplot(1,3,3),plot(x3,Y,'ro')

绘制的图形如下：

（2）进行多元线性回归

这里可以直接使用 regress 函数执行多元线性回归，注意以下代码模板，以后碰到多元线性问题直接套用代码，具体代码如下：

%% 进行多元线性回归
n = 24; m = 3; % 每个变量均有24个数据，共有3个变量
X = [ones(n,1),x1',x2',x3'];
[b,bint,r,rint,s]=regress(Y',X,0.05) % 0.05为预定显著水平，判断因变量y与自变量之间是否具有显著的线性相关关系需要用到。

运行结果如下：

b =

   18.0157
    1.0817
    0.3212
    1.2835


bint =

   13.9052   22.1262
    0.3900    1.7733
    0.2440    0.3984
    0.6691    1.8979


r =

    0.6781
    1.9129
   -0.1119
    3.3114
   -0.7424
    1.2459
   -2.1022
    1.9650
   -0.3193
    1.3466
    0.8691
   -3.2637
   -0.5115
   -1.1733
   -1.4910
   -0.2972
    0.1702
    0.5799
   -3.2856
    1.1368
   -0.8864
   -1.4646
    0.8032
    1.6301


rint =

   -2.7017    4.0580
   -1.6203    5.4461
   -3.6190    3.3951
    0.0498    6.5729
   -4.0560    2.5712
   -2.1800    4.6717
   -5.4947    1.2902
   -1.3231    5.2531
   -3.5894    2.9507
   -1.7678    4.4609
   -2.7146    4.4529
   -6.4090   -0.1183
   -3.6088    2.5859
   -4.7040    2.3575
   -4.8249    1.8429
   -3.7129    3.1185
   -3.0504    3.3907
   -2.8855    4.0453
   -6.2644   -0.3067
   -2.1893    4.4630
   -4.4002    2.6273
   -4.8991    1.9699
   -2.4872    4.0937
   -1.8351    5.0954


s =

    0.9106   67.9195    0.0000    3.0719

看到如此长的运行结果，我们不要害怕，因为里面很多数据是没用的，我们只需提取有用的数据。

在运行结果中，很多数据我们不需理会，我们真正需要用到的数据如下：

b =

   18.0157
    1.0817
    0.3212
    1.2835

s =

    0.9106   67.9195    0.0000    3.0719

回归系数 b = (β0,β1,β2,β3) = (18.0157, 1.0817, 0.3212, 1.2835)，回归系数的置信区间，以及统计变量 stats（它包含四个检验统计量：相关系数的平方R^2，假设检验统计量 Ｆ，与 F 对应的概率 p，s^2 的值）。观察表4的数据，会发现它来源于运行结果中的b和s:

根据β0,β1,β2,β3,我们初步得出回归方程为：

如何判断该回归方程是否符合该模型呢？有以下3种方法：

1）相关系数 Ｒ 的评价：本例 Ｒ 的绝对值为 0.9542 ，表明线性相关性较强。

2）F 检验法：当 F > F1-α(m,n-m-1) ，即认为因变量 y 与自变量 x1,x2,...,xm 之间有显著的线性相关关系；否则认为因变量 y 与自变量 x1,x2,...,xm 之间线性相关关系不显著。本例 Ｆ＝67.919 > F1-0.05( 3,20 ) = 3.10。

3）p 值检验：若 p < α（α 为预定显著水平），则说明因变量 y 与自变量 x1,x2,...,xm之间显著地有线性相关关系。本例输出结果，p<0.0001,显然满足 p<α=0.05。

以上三种统计推断方法推断的结果是一致的，说明因变量 y 与自变量之间显著地有线性相关关系，所得线性回归模型可用。s^2 当然越小越好，这主要在模型改进时作为参考。

3. 逐步回归

[ 例4 ] （Hald,1960）Hald 数据是关于水泥生产的数据。某种水泥在凝固时放出的热量 Y（单位：卡/克）与水泥中 4 种化学成品所占的百分比有关：

在生产中测得 12 组数据，见表5，试建立 Y 关于这些因子的“最优”回归方程。

对于例 4 中的问题，可以使用多元线性回归、多元多项式回归，但也可以考虑使用逐步回归。从逐步回归的原理来看，逐步回归是以上两种回归方法的结合，可以自动使得方程的因子设置最合理。对于该问题，逐步回归的代码如下：

%% 逐步回归
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];   %自变量数据
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];  %因变量数据
stepwise(X,Y,[1,2,3,4],0.05,0.10)% in=[1,2,3,4]表示X1、X2、X3、X4均保留在模型中

程序执行后得到下列逐步回归的窗口，如图 4 所示。

图4

在图 4 中，用蓝色行显示变量 X1、X2、X3、X4 均保留在模型中，窗口的右侧按钮上方提示：将变量X4剔除回归方程（Move X4 out），单击 Next Step 按钮，即进行下一步运算，将第 4 列数据对应的变量 X4 剔除回归方程。单击 Next Step 按钮后，剔除的变量 X3 所对应的行用红色表示，同时又得到提示：将变量 X3 剔除回归方程（Move X3 out），单击 Next Step 按钮，这样一直重复操作，直到 “Next Step” 按钮变灰，表明逐步回归结束，此时得到的模型即为逐步回归最终的结果。最终结果如下：

4. 逻辑回归

[ 例5 ] 企业到金融商业机构贷款，金融商业机构需要对企业进行评估。评估结果为 0 , 1 两种形式，0 表示企业两年后破产，将拒绝贷款，而 1 表示企业 2 年后具备还款能力，可以贷款。在表 6 中，已知前 20 家企业的三项评价指标值和评估结果，试建立模型对其他 5 家企业（企业 21-25）进行评估。

对于该问题，很明显可以用 Logistic 模型来回归，具体求解程序如下：

程序中需要用到的数据文件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

% logistic回归

%% 导入数据
clc,clear,close all
X0 = xlsread('logistic_ex1.xlsx','A2:C21'); % 前20家企业的三项评价指标值,即回归模型的输入
Y0 = xlsread('logistic_ex1.xlsx','D2:D21'); % 前20家企业的评估结果，即回归模型的输出
X1 = xlsread('logistic_ex1.xlsx','A2:C26'); % 预测数据输入

%% 逻辑函数
GM = fitglm(X0,Y0,'Distribution','binomial');
Y1 = predict(GM,X1);

%% 模型的评估
N0 = 1:size(Y0,1); % N0 = [1,2,3,4,……，20]
N1 = 1:size(Y1,1); % N1 = [1,2,3,4,……，25]
plot(N0',Y0,'-kd'); % N0'指的是对N0'进行转置，N0'和Y0的形式相同，该行代码绘制的是前20家企业的评估结果
% plot()中的参数'-kd'的解析：-代表直线，k代表黑色，d代表菱形符号
hold on;
scatter(N1',Y1,'b'); % N1'指的是对N1'进行转置，N1'和Y1的形式相同
xlabel('企业编号');
ylabel('输出值');

得到的回归结果与原始数据的比较如图5所示。

图5

三、总结与感悟。

总结：通过这次学习，我了解到Matlab在数学建模竞赛中使用广泛；在评估股票价值与风险的小实例中，我掌握了用Matlab去建模的基本方法和步骤；在回归算法的学习过程中，我掌握了一元线性回归、一元非线性回归、多元线性回归、逐步回归、逻辑回归的算法。

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