数学建模竞赛中的两个技巧

admin管理员组

文章数量:1660200

文章目录

• • • 1. 数据拟合与最优化方法
    • • 1.1. 数据拟合
      • • 1.1.1. 线性回归 fitlm
        • 1.1.2. 多项式线性回归 regress
        • 1.1.3. 多项式非线性回归 polyfit
        • 1.1.4. 自定义非线性回归 fit
        • 1.1.5. 自定义非线性最小二乘法回归 lsqcurvefit
        • 1.1.6. 自定义非线性牛顿法回归 nlinfit
        • 1.1.7. 拟合工具箱 Curve Fitting
      • 1.2. 最优化方法
      • • 1.2.1. 无约束优化问题
        • 1.2.2. 约束优化问题
    • 2. 神经网络与特征选择
    • • 2.1. 神经网络
      • • 2.1.1. 分类神经网络
        • 2.1.2. 回归神经网络
      • 2.2. 特征选择

1. 数据拟合与最优化方法

1.1. 数据拟合

数据拟合是为了连续化离散的数据，以方便微积分工具在数据上的使用。

例如有数据x代表温度，y代表产量：

x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

如图：

怎么去计算温度 x = 5.5 x=5.5 x=5.5时，产量 y y y的大小呢？这就需要用到数据拟合了，通过拟合来将离散的数据点串成连续的函数，就可以计算出函数上任意一点的值。当然拟合的作用不止于此。

下面我们用一组二维数据和一组三维数据来展示MatLab中多种回归拟合模型的使用。

二维数据
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

三维数据
x1 = [1 2 3 4 5 6 7 8 9];
x2 = [5 4 3 2 0 9 10 15 19];
y = [9 7 6 3 -1 2 5 7 20];

1.1.1. 线性回归 fitlm

https://ww2.mathworks/help/stats/fitlm.html?lang=en

film的输入：

1. X，即自变量矩阵；
2. Y，即因变量值。

film的输出：

1. 回归表达式(Linear regression model)，即线性回归表达式；
2. 常数项(Intercept)，即回归表达式中常数项的值；
3. 系数估计值(Estimate)，即回归表达式中的系数值；
4. 标准误差(Standard Error, SE)，即系数的标准误差；
5. t统计量(tStat)，即每个系数的 t 统计量，tStat = Estimate/SE；
6. t 统计量的 p 值(pValue)，即假设检验的 t 统计量的 p 值；
7. 其余参数查看matlab官网 https://ww2.mathworks/help/stats/fitlm.html?lang=en。

% 拟合x1与y
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

mdl = fitlm(x1,y);
disp(mdl);

输出：

Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                   Estimate      SE        tStat     pValue 
                   ________    _______    _______    _______

    (Intercept)     3.0278      4.3607    0.69434    0.50985
    x1             0.68333     0.77491    0.88182    0.40713

表达式为： y = 3.0278 + 0.6833 x 1 y=3.0278+0.6833x_1 y=3.0278+0.6833x1​

% 拟合x1、x2与y
x1 = [1 2 3 4 5 6 7 8 9];
x2 = [5 4 3 2 0 9 10 15 19];
y = [9 7 6 3 -1 2 5 7 20];

mdl = fitlm([x1', x2'],y); % x1'是x1向量转置的意思
disp(mdl);

输出：

Linear regression model:
    y ~ 1 + x1 + x2

Estimated Coefficients:
                   Estimate      SE        tStat      pValue 
                   ________    _______    _______    ________

    (Intercept)     5.0714      2.7551     1.8407     0.11526
    x1             -1.4922      0.7824    -1.9072     0.10512
    x2              1.1866     0.33762     3.5147    0.012599


Number of observations: 9, Error degrees of freedom: 6
Root Mean Squared Error: 3.71
R-squared: 0.706,  Adjusted R-Squared 0.608
F-statistic vs. constant model: 7.2, p-value = 0.0255

表达式为： y = 5.0714 − 1.4922 x 1 + 1.1866 x 2 y=5.0714-1.4922x_1+1.1866x_2 y=5.0714−1.4922x1​+1.1866x2​

1.1.2. 多项式线性回归 regress

https://ww2.mathworks/help/stats/regress.html?lang=en

regress的输入：

1. Y，即因变量值;
2. X_matrix，即自变量表达式矩阵，可自定义。

regress的输出：

1. 系数估计值(b)，即多元线性回归的系数估计值；
2. 95% 置信区间(bint)，即系数估计值的 95% 置信区间的矩阵；
3. 其余参数查看matlab官网 https://ww2.mathworks/help/stats/regress.html?lang=en。

```python
% 拟合x1与y
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

X = [ones(size(x1')), x1'];
[b, bint] = regress(y', X);
disp(b);

输出：

b=
    3.0278
    0.6833

表达式为： y = 3.0278 + 0.6833 x 1 y=3.0278+0.6833x_1 y=3.0278+0.6833x1​

% 拟合x1、x2与y
x1 = [1 2 3 4 5 6 7 8 9];
x2 = [5 4 3 2 0 9 10 15 19];
y = [9 7 6 3 -1 2 5 7 20];

X = [ones(size(x1')), x1'.^3, x2']; % 注意这里给自变量x1加上了3次方
[b, bint] = regress(y', X);
disp(b);

输出：

b=
    0.9233
   -0.0089
    1.0100

表达式为： y = 0.9233 − 0.0089 x 1 3 + 1.0100 x 2 y=0.9233-0.0089x_1^3+1.0100x_2 y=0.9233−0.0089x13​+1.0100x2​

1.1.3. 多项式非线性回归 polyfit

https://ww2.mathworks/help/matlab/ref/polyfit.html?lang=en

polyfit的输入：

1. X，即自变量值；
2. Y，即因变量值；
3. n，即多项式次数；

polyfit的输出：

1. 系数估计值§，即次数为 n 的多项式 p(x) 的系数；
2. 误差估计值(s)，s是一个结构体，记录了误差估计值信息；
3. 矩 μ \mu μ，分别记录了一阶矩 μ ( 1 ) \mu(1) μ(1)和二阶矩 μ ( 2 ) \mu(2) μ(2)；
4. 其余参数查看matlab官网 https://ww2.mathworks/help/matlab/ref/polyfit.html?lang=en。

% 拟合x1与y
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

[p, S, mu] = polyfit(x1, y, 1);
disp(p);

输出：

p=
1.8714    6.4444

表达式为： y = 1.8714 ∗ x 1 + 6.4444 ∗ x 0 y=1.8714*x^1+6.4444*x^0 y=1.8714∗x1+6.4444∗x0

% 拟合x1与y
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

[p, S, mu] = polyfit(x2, y, 2); % 这里变成了2
disp(p);

输出：

p=
6.1445    1.8714    0.9827

表达式为： y = 6.1445 ∗ x 1 2 + 1.8714 ∗ x 1 1 + 0.9827 ∗ x 1 0 y=6.1445*x_1^2+1.8714*x_1^1+0.9827*x_1^0 y=6.1445∗x12​+1.8714∗x11​+0.9827∗x10​

1.1.4. 自定义非线性回归 fit

https://ww2.mathworks/help/curvefit/fit.html

fit的输入：

1. X，即自变量值；
2. Y，即因变量值；
3. fitType，即拟合函数；

fit的输出：

1. 拟合模型(fitobject)
2. 其余参数查看matlab官网 https://ww2.mathworks/help/curvefit/fit.html。

% 拟合x1与y
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

fitobject = fit(x1', y', 'poly2'); % ploy2表示x1的最高项次数为2
disp(fitobject);

输出：

Linear model Poly2:
     fitobject(x) = p1*x^2 + p2*x + p3
     Coefficients (with 95% confidence bounds):
       p1 =      0.8193  (0.4354, 1.203)
       p2 =      -7.509  (-11.45, -3.573)
       p3 =       18.05  (9.476, 26.62)

表达式为： y = 0.8193 ∗ x 1 2 − 7.509 ∗ x 1 1 + 18.05 y=0.8193*x_1^2-7.509*x_1^1+18.05 y=0.8193∗x12​−7.509∗x11​+18.05

% 拟合x1、x2与y
x1 = [1 2 3 4 5 6 7 8 9];
x2 = [5 4 3 2 0 9 10 15 19];
y = [9 7 6 3 -1 2 5 7 20];

fitobject = fit([x1', x2'], y', 'poly23'); % ploy23表示x1的最高项次数为2，x2的最高项次数为3
disp(fitobject);

输出：

     Linear model Poly23:
     fitobject(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2 + p21*x^2*y 
                    + p12*x*y^2 + p03*y^3
     Coefficients:
       p00 =      -286.4
       p10 =       121.3
       p01 =       62.76
       p20 =      -12.85
       p11 =      -21.57
       p02 =     -0.4537
       p21 =       1.523
       p12 =      0.7509
       p03 =     -0.2724

表达式为： y = − 286.4 + 121.3 ∗ x 1 + 62.76 ∗ x 2 + − 12.85 ∗ x 2 + − 21.57 ∗ x 1 ∗ x 2 + − 0.4537 ∗ x 2 2 + 1.523 ∗ x 1 ∗ x 2 + 0.7509 ∗ x 1 ∗ x 2 2 + − 0.2724 ∗ x 2 3 y=-286.4 + 121.3*x_1 + 62.76*x_2 + -12.85*x^2 + -21.57*x_1*x_2 + -0.4537*x_2^2 + 1.523*x_1*x_2 + 0.7509*x_1*x_2^2 + -0.2724*x_2^3 y=−286.4+121.3∗x1​+62.76∗x2​+−12.85∗x2+−21.57∗x1​∗x2​+−0.4537∗x22​+1.523∗x1​∗x2​+0.7509∗x1​∗x22​+−0.2724∗x23​

1.1.5. 自定义非线性最小二乘法回归 lsqcurvefit

https://ww2.mathworks/help/optim/ug/lsqcurvefit.html

lsqcurvefit的输入：

1. fun，即拟合函数；
2. x0，即初始化X值。
3. X，即自变量值；
4. Y，即因变量值；
5. lb，自变量的下限
6. ub，自变量的上限

lsqcurvefit的输出：

1. 自变量系数(x)
2. 其余参数查看matlab官网 https://ww2.mathworks/help/optim/ug/lsqcurvefit.html。

% 拟合x1与y
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

beta0 = [1 1 1];
fun = @(beta0, x1)( beta0(1).*x1.^2 + beta0(2).*x1 + beta0(3));
lb = [0];
ub = [10];
x_result = lsqcurvefit(fun, beta0, x1, y, lb, ub);
disp(x_result);

输出：

x=
    0.8193   -7.5093   18.0476

表达式为： y = 0.8193 ∗ x 1 2 − 7.509 ∗ x 1 1 + 18.0476 y=0.8193*x_1^2-7.509*x_1^1+18.0476 y=0.8193∗x12​−7.509∗x11​+18.0476

% 拟合x1、x2与y
x1 = [1 2 3 4 5 6 7 8 9];
x2 = [5 4 3 2 0 9 10 15 19];
y = [9 7 6 3 -1 2 5 7 20];

beta0 = [1 1 1];
X = [x1' x2'];
fun = @(beta0, X)( beta0(1).*X(:,1).^2 + beta0(2).*X(:,2).^2 + beta0(3));
lb = [0, -2];
ub = [10, 20];
x_result = lsqcurvefit(fun, beta0, X, y', lb, ub);
disp(x_result);

输出：

x=
    0.0000    0.0370    3.0720

表达式为： y = 0.0000 ∗ x 1 2 − 0.0370 ∗ x 2 2 + 3.0720 y=0.0000*x_1^2-0.0370 *x_2^2+3.0720 y=0.0000∗x12​−0.0370∗x22​+3.0720

1.1.6. 自定义非线性牛顿法回归 nlinfit

https://ww2.mathworks/help/stats/nlinfit.html?lang=en

nlinfit的输入：

1. X，即自变量值；
2. Y，即因变量值；
3. model，即拟合函数；
4. beta0，即初始化X值。

nlinfit的输出：

1. 自变量系数(beta)
2. 其余参数查看matlab官网 https://ww2.mathworks/help/stats/nlinfit.html?lang=en。

% 拟合x1与y
x1 = [1 2 3 4 5 6 7 8 9];
y = [9 7 6 3 -1 2 5 7 20];

beta0 = [1 1 1];
fun = @(beta0, x1)( beta0(1).*x1.^2 + beta0(2).*x1 + beta0(3));
beta = nlinfit(x1, y, fun, beta0);
disp(beta);

输出：

beta=
    0.8193   -7.5093   18.0476

表达式为： y = 0.8193 ∗ x 1 2 − 7.509 ∗ x 1 1 + 18.0476 y=0.8193*x_1^2-7.509*x_1^1+18.0476 y=0.8193∗x12​−7.509∗x11​+18.0476

% 拟合x1、x2与y
x1 = [1 2 3 4 5 6 7 8 9];
x2 = [5 4 3 2 0 9 10 15 19];
y = [9 7 6 3 -1 2 5 7 20];

beta0 = [1 1 1];
X = [x1' x2'];
fun = @(beta0, X)( beta0(1).*X(:,1).^2 + beta0(2).*X(:,2).^2 + beta0(3));
beta = nlinfit(X, y', fun, beta0);
disp(beta);

输出：

beta=
   -0.2973    0.0988    6.8492

表达式为： y = − 0.2973 ∗ x 1 2 − 0.0988 ∗ x 2 2 + 6.8492 y=-0.2973*x_1^2-0.0988*x_2^2+6.8492 y=−0.2973∗x12​−0.0988∗x22​+6.8492

1.1.7. 拟合工具箱 Curve Fitting

Curve Fitting是Matlab中一个方便快捷的函数拟合工具箱，优点就是方便，直观，可以自动出图，而且图像还可以任意角度旋转查看，但是无法拟合高于3维的数据。

1.2. 最优化方法

在有了函数表达式之后，如何求最优值通常是一个巨大的挑战，所以学习一些最优化知识是必要的。

1.2.1. 无约束优化问题

所谓无约束优化问题，就是指对一个函数求最优值，最优值可以出现在函数上任意一点，而不去限定查找最优值的范围。

无约束优化问题通常有四种优化方法：

1. Newton’s method(牛顿法)；
2. Levenberg-Marquardt’s method(LM)；
3. Broyden-Fletcher-Goldfarb-Shanno’s method(BFGS)；
4. Davidon-Fletcher-Powell’s method(DFP)

假设有如下无约束优化问题： f ( x 1 , x 2 ) = ( x 1 − 3 ) 4 + ( x 1 − 3 x 2 ) 2 f(x_1,x_2)=(x_1-3)^4+(x_1-3x_2)^2 f(x1​,x2​)=(x1​−3)4+(x1​−3x2​)2

下面我们来用matlab实现一下这四种优化算法：
无约束优化问题中牛顿法与拟牛顿法四种迭代方法的matlab实现

1.2.2. 约束优化问题

约束优化问题就是指给自变量 x x x的取值范围做限制，缩小优化范围，经典的约束优化算法有：

1. 内点法(interior-point)
2. 有效集法(active-set)
3. SQP算法(sqp)
4. 信赖域反射法(trust-region-reflective)

约束，又分线性约束与非线性约束，所谓线性约束，就是指约束条件中的自变量都是1次幂的，非线性约束即有高次幂的自变量。

假设有如下约束优化函数和约束：
min ⁡ f ( x 1 , x 2 ) = − x 1 x 2 \min f(x_1,x_2)=-x_1x_2 minf(x1​,x2​)=−x1​x2​ s . t . 1 − x 1 2 − x 2 2 ⩾ 0 s.t.\quad 1-x_1^2-x_2^2\geqslant 0 s.t.1−x12​−x22​⩾0

那么优化器方法可以定义为：

function [xsol, fval] = runfmincon
% 初始点
x0 = [-0.1 -0.1];

% 四种优化算法，选一种用，这里用的是内点法
% 'active-set', 'interior-point', 'sqp', or 'trust-region-reflective'.
options = optimset('Display', 'iter-detailed', 'Algorithm', 'interior-point', 'MaxIter', 8); % 优化器参数
% fmincon参数解释
% fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
% fun:要优化的函数
% x0:自变量的初始值
% A:非等式(<)线性约束的约束矩阵A
% b:非等式(<)线性约束的约束条件矩阵b
% Aeq:等式约束的线性约束的约束矩阵Aeq
% beq:等式线性约束的约束条件矩阵beq
% lb,ub:自变量的下限和上限
% nonlcon:非线性约束
[xsol, fval] = fmincon(@objfun, x0, [], [], [], [], [], [], @confun, options);
 
% 目标函数
 function f = objfun(x)
     f = - x(1) * x(2);
 end
 
% 非线性约束
 function [c, ceq] = confun(x)
     % Nonlinear inequality constraints
     c = [x(1)^2 + x(2)^2 - 1];
     % Nonlinear equality constraints
     ceq = [];
 end
end

优化结果：

[xsol,fval] = runfmincon;

% 迭代优化结果1
[x1,x2]=[-0.7071,-0.7071]
f=-0.5
% 迭代优化结果2
[x1,x2]=[0.7071,0.7071]
f=-0.5

2. 神经网络与特征选择

2.1. 神经网络

神经网络有两种功能，一是分类，二是回归，下面分别是两种神经网络的python代码。通常会使用pytorch框架来实现神经网络，可以参考博客：从细节过渡到实例 一天学会Pytorch

2.1.1. 分类神经网络

import torch
import torch.nn.functional as F
import matplotlib.pyplot as plt

# 构造数据
n_data = torch.ones(100, 2)
x0 = torch.normal(3*n_data, 1)
x1 = torch.normal(-3*n_data, 1)
# 标记为y0=0，y1=1两类标签
y0 = torch.zeros(100)
y1 = torch.ones(100)

# 通过.cat连接数据
x = torch.cat((x0, x1), 0).type(torch.float)
y = torch.cat((y0, y1), 0).type(torch.long)

# 构造一个简单的神经网络
class Net(torch.nn.Module):
    def __init__(self, n_feature, n_hidden, n_output):
        super(Net, self).__init__()
        self.inLayer = torch.nn.Linear(n_feature, n_hidden) # 输入层
        self.hiddenLayer = torch.nn.Linear(n_hidden, n_hidden) # 隐藏层
        self.outLayer = torch.nn.Linear(n_hidden, n_output) # 输出层

    # 前向计算函数，定义完成后会隐式地自动生成反向传播函数
    def forward(self, x):
        x = F.relu(self.hiddenLayer(self.inLayer(x)))
        x = self.outLayer(x)
        return x

net = Net(2, 10, 2) # 初始化一个网络，2个输入层节点，10个隐藏层节点，2个输出层节点
loss_func = torch.nn.CrossEntropyLoss() # 定义交叉熵损失函数
optimizer = torch.optim.SGD(net.parameters(), lr=0.2) # 配置网络优化器

# 将网络模型、损失函数和输入张量迁入GPU
if(torch.cuda.is_available()):
    net = net.cuda()
    loss_func = loss_func.cuda()
    x, y = x.cuda(), y.cuda()

# 训练模型
out = net(x)
for t in range(300):
    out = net(x) # 将数据输入网络，得到输出
    loss = loss_func(out, y) # 得到损失
    optimizer.zero_grad() # 梯度初始化
    loss.backward() # 反向传播
    optimizer.step() # 对网络进行优化

# 使用模型进行预测
prediction = torch.max(F.softmax(out, dim=0), 1)[1]
pred_y = prediction.data.cpu().numpy().squeeze()

# 可视化
plt.scatter(x.data.cpu().numpy()[:, 0], x.data.cpu().numpy()[:, 1], c=pred_y, s=100, lw=0, cmap='RdYlBu')
plt.show()

2.1.2. 回归神经网络

import torch
import torch.nn.functional as F
import matplotlib.pyplot as plt

# 构造数据
x = torch.unsqueeze(torch.linspace(-1, 1, 100), dim=1)
y = x.pow(2) + 0.2 * torch.rand(x.size())

# 构造一个简单的神经网络
class Net(torch.nn.Module):
    def __init__(self, n_feature, n_hidden, n_output):
        super(Net, self).__init__()
        self.inLayer = torch.nn.Linear(n_feature, n_hidden) # 输入层
        self.hiddenLayer = torch.nn.Linear(n_hidden, n_hidden) # 隐藏层
        self.outLayer = torch.nn.Linear(n_hidden, n_output) # 输出层

    # 前向计算函数，定义完成后会隐式地自动生成反向传播函数
    def forward(self, x):
        x = F.relu(self.hiddenLayer(self.inLayer(x)))
        x = self.outLayer(x)
        return x

net = Net(1, 10, 1) # 初始化一个网络，1个输入层节点，10个隐藏层节点，1个输出层节点
loss_func = torch.nn.MSELoss() # 定义均方误差损失函数
optimizer = torch.optim.SGD(net.parameters(), lr=0.2) # 配置网络优化器

# 将网络模型、损失函数和输入张量迁入GPU
if(torch.cuda.is_available()):
    net = net.cuda()
    loss_func = loss_func.cuda()
    x, y = x.cuda(), y.cuda()

# 训练模型
out = net(x)
for t in range(300):
    out = net(x) # 将数据输入网络，得到输出
    loss = loss_func(out, y) # 得到损失
    optimizer.zero_grad() # 梯度初始化
    loss.backward() # 反向传播
    optimizer.step() # 对网络进行优化

# 可视化
plt.scatter(x.data.cpu().numpy(), y.data.cpu().numpy())
plt.plot(x.data.cpu().numpy(), out.data.cpu().numpy(), 'r-', lw=5)
plt.show()

2.2. 特征选择

尽管神经网络是一个优秀的分类回归器，但是也不意为着它是万能的，在神经网络拟合的效果不好的情况下，还需要先进行特征筛选，也就是选出最适合的特征，再将它们输入神经网络进行训练。一个优秀的特征选择器是决策树，或者随机森林。

scikit-learn：决策树
scikit-learn：随机森林

在决策树中，有几个重要的概念就是“信息熵、条件熵、信息增益及信息增益比”

1. 信息熵
   信息熵表示了信息的不确定性，也就是一个特征对结果影响的不确定性；
2. 条件熵
   条件熵指的是在其他特征作用的情况下，某特征对结果影响的不确定性。
3. 信息增益
   信息增益 = 信息熵 - 条件熵。信息增益指的是在其他特征作用的情况下，某特征对结果的不确定性的减小程度。
4. 信息增益比
   信息增益比是一种对信息增益的归一化。

信息熵、条件熵、信息增益及信息增益比的应用举例
决策树中熵、条件熵、信息增益及信息增益比的python实现

本文标签： 建模两个数学技巧