个人在数学建模过程中的学习梳理的笔记，编程软件和语言以RStudio/R为主。
供学习参考。部分内容待补充。

Chapter1 预测模型

1.1 线性回归模型

总体步骤

注¹

step1【 建立模型】

• 使用内置函数 lm(Y~X,data=bp) 函数

step2【 F、t、拟合优度检验】

• 在建立模型后，我们需要对该模型进行检验，该检验包括：从整体效果上来看，该模型是否能够接受？从每个解释变量上来看，是否必要？同时，此模型的拟合效果又如何？
• 为讨论上述问题，需要使用 summary(lm.bp) 函数对模型进行以下检验：
  1. 模型整体效果评价（回归方程显著性检验）：F检验（ H 0 : β 1 = β 2 = . . . = 0 H_0:\beta_1=\beta_2=...=0 H0​:β1​=β2​=...=0）： p < α p<\alpha p<α 时检验通过；
  2. 解释变量个体必要评价（回归系数显著性检验）：t检验（ H 0 : β i = 0 ( i = 1 , . . . , p ) H_0:\beta_i=0(i=1,...,p) H0​:βi​=0(i=1,...,p)）： p < α p<\alpha p<α 时检验通过；
  3. 拟合优度检验： a d j u s t R α 2 adjust\space R^{2}_{\alpha} adjust Rα2​ 越接近1 整体拟合效果（模型所能解释的）越好

step3 【异常值检验和处理】

• 针对 Y Y Y 使用 学生化残差 rstudent(lm.bp) ，当所有 ∣ r s ∣ < 3 |rs|<3 ∣rs∣<3 时检验通过
• 针对 X X X 使用 库克距离 cooks.distance(lm.bp) ，当所有 ∣ D i ∣ < 1 |D_i|<1 ∣Di​∣<1 时检验通过

step4 【异方差与自相关检验】

• 首先要明确异方差的定义：残差项随自变量的变动情况

• 我们对模型有如下要求： E ( δ i ) = 0 E(\delta_i)=0 E(δi​)=0， c o v ( δ i , δ j ) = 0 cov(\delta_i,\delta_j)=0 cov(δi​,δj​)=0， D ( δ i ) = 0 D(\delta_i)=0 D(δi​)=0，也即：零均值，不相关，等方差

• 对于异方差检验，我们采取 cor.test() 函数的spearman检验，当 p > α p>\alpha p>α 时检验通过。具体形式代码如下：

  cor.test(abs(lm.bp$residual),bp$X,method='spearman')
  

• 对于自相关检验，我们采取 lmtest 宏包中的 dwtest() 函数，当 p > α p>\alpha p>α 时检验通过。具体形式代码如下：

  library(lmtest)
  dwtest(lm.bp)
  

另有，对于dw检验，若dw值在 2 2 2 附近则可初步认为残差不相关（经验值判定法）

step5 【模型优化】

• 模型的优化一般有以下几种情况：

  1. 在进行summary检验时，三类检验发现效果不好，可以考虑剔除某个变量（谨慎）；
  2. 模型未通过异方差和自相关检验（即残差不是同方差或残差自相关）；
  3. 某些变量可能通过了t检验，但由于其存在可能会导致模型复杂化，选择优化删去；或者是有些变量本身就不显著，此时不能轻易剔除，选择逐步回归；

• 对于 c a s e Ⅱ caseⅡ caseⅡ ，我们考虑 box-cox 变换：

  何谓 box-cox 变换？是一种将变量转化为近似正态分布的数学变换。
  y ( λ ) = { y i λ − 1 λ , i f λ ≠ 0 ; l n ( y i ) , i f λ = 0 ; y(\lambda)= \begin{cases} &\frac{y_i^\lambda-1}{\lambda}&,if\space \lambda\neq0;\\ &ln(y_i)& ,if\space \lambda=0; \end{cases} y(λ)={​λyiλ​−1​ln(yi​)​,if λ​=0;,if λ=0;​

  box-cox变换主要有两项工作，第一是做变换，第二是确定 λ \lambda λ 的值。要用最大似然估计才能确定 λ \lambda λ 的值。

  具体形式代码如下：

  library(MASS)
  box<-boxcox(lm.bp,lambda=seq(-2,2,0.1))#λ取值为区间[-2,2]上步长为0.01的值，box中保存了λ的值及其对应的对数似然函数值
  lambda<-box$x[which.max(box$y)]#将使对数似然函数值达到最大的λ赋值给lambda
  lm((Y^lambda-1)/lambda~X,data=bp)->lm.bp2
  

  部分内容引用自：https://www.pianshen.com/article/4849182009/

• 对于 c a s e Ⅲ caseⅢ caseⅢ ，我们考虑 AIC分析：

  AIC分析 ：AIC 值越小，说明模型效果越好。AIC 分析用于剔除与 y y y 相关性很弱的 x i x_i xi​ ，简化模型。

  lm.sbp<-step(lm.bp)
  

• 同时需要明确，不能仅仅根据t检验得到的显著性水平去剔除解释变量，还需要结合实际意义及分析目标来判断。更极端的说，甚至有可能某个模型每个解释变量本身对于线性模型不显著，但所有的解释变量加起来对因变量有很强的线性影响。所以我通常不说接受原假设，而是说无法拒绝原假设，不显著并不等同于没有影响。

step5 【模型预测】

• 在给定自变量的情况下，利用 predict() 函数可以根据已有的回归模型给出因变量的一个预测值，可以求出相应的置信区间和预测区间

preds<-data.frame(x=200)  #给定解释变量x的值
predict(lm.sbp,newdata=preds,interval="c",level=0.95)#置信区间
predict(lm.sbp,newdata=preds,interval="prediction",level=0.95)#预测区间

案例分析

数据源：土地含磷量.data

下为某地区土壤内含 植物可给态磷 y y y 与土壤内所含无机磷浓度 x 1 x_1 x1​ ，土壤内溶于 K 2 C O 3 K_2CO_3 K2​CO3​ 溶液并受溴化物水解的有机磷浓度 x 2 x_2 x2​ 以及土壤内溶于 K 2 C O 3 K_2CO_3 K2​CO3​ 溶液但不溶于溴化物的有机磷 x 3 x_3 x3​ 的观察数据。求 y y y 对 x 1 , x 2 , x 3 x_1,x_2 ,x_3 x1​,x2​,x3​ 的最优线性回归方程

     X1 X2  X3   Y
1   0.4 52 158  64
2   0.4 23 163  60
3   3.1 19  37  71
4   0.6 34 157  61
5   4.7 24  59  54
6   1.7 65 123  77
7   9.4 44  46  81
8  10.1 31 117  93
9  11.6 29 173  93
10 12.6 58 112  51
11 10.9 37 111  76
12 23.1 46 114  96
13 23.1 50 134  77
14 21.6 44  73  93
15 23.1 56 168  95
16  1.9 36 143  54
17 26.8 58 202 168
18 29.9 51 124  99

解答：

sp<-read.table(file.choose())
lm.sp<-lm(Y~X1+X2+X3,data=sp)

summary(lm.sp)

Call:
lm(formula = Y ~ X1 + X2 + X3, data = sp)Residuals:Min      1Q  Median      3Q     Max
-28.349 -11.383  -2.659  12.095  48.807 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 43.65007   18.05442   2.418  0.02984 *
X1           1.78534    0.53977   3.308  0.00518 **
X2          -0.08329    0.42037  -0.198  0.84579
X3           0.16102    0.11158   1.443  0.17098
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 19.97 on 14 degrees of freedom
Multiple R-squared:  0.5493,    Adjusted R-squared:  0.4527
F-statistic: 5.688 on 3 and 14 DF,  p-value: 0.009227

若显著性水平² α = 0.05 \alpha=0.05 α=0.05，由上述结果可知，回归方程为：
y = 43.65007 + 1.78534 ∗ x 1 − 0.08329 ∗ x 2 + 0.16102 ∗ x 3 y=43.65007+1.78534*x_1-0.08329*x_2+0.16102*x_3 y=43.65007+1.78534∗x1​−0.08329∗x2​+0.16102∗x3​
由F检验， p = 0.009227 < 0.05 p=0.009227<0.05 p=0.009227<0.05，说明模型整体效果较好，也即 x 1 , x 2 , x 3 x_1,x_2,x_3 x1​,x2​,x3​ 对 y y y 有显著的线性影响；

由拟合优度 A d j u s t e d R − s q u a r e d : 0.4527 Adjusted R-squared: 0.4527 AdjustedR−squared:0.4527 可知拟合效果较好，有 45.27 % 45.27\% 45.27% ；

再由t检验，发现 x 2 , x 3 x_2,x_3 x2​,x3​ 均不显著，下去使用逐步回归法得到最优回归模型：

lm.ssp<-step(lm.sp)
summary(lm.ssp)

Call:
lm(formula = Y ~ X1 + X3, data = sp)Residuals:Min      1Q  Median      3Q     Max
-29.713 -11.324  -2.953  11.286  48.679 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept)  41.4794    13.8834   2.988  0.00920 **
X1            1.7374     0.4669   3.721  0.00205 **
X3            0.1548     0.1036   1.494  0.15592
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 19.32 on 15 degrees of freedom
Multiple R-squared:  0.5481,    Adjusted R-squared:  0.4878
F-statistic: 9.095 on 2 and 15 DF,  p-value: 0.002589

此时有最优回归方程：
y = 41.4394 + 1.7374 x 1 + 0.1548 x 3 y=41.4394+1.7374x_1+0.1548x_3 y=41.4394+1.7374x1​+0.1548x3​

下面我们再进行异常值检验：

rstudent(lm.ssp)

          1           2           3           4           5           6           7           8           9
-0.14626939 -0.41640212  1.13680072 -0.32308014 -0.26560994  0.73311566  0.93916538  0.83713927  0.24572879 10          11          12          13          14          15          16          17          18
-1.67932661 -0.08278959 -0.17574935 -1.45643942  0.15050863 -0.70384407 -0.70661467  4.78264426 -0.80110979

学生化残差检验通过；

cooks.distance(lm.ssp)

           1            2            3            4            5            6            7            8
0.0017266538 0.0150647846 0.1676799013 0.0080769441 0.0057210303 0.0244556277 0.0846742907 0.0147139595 9           10           11           12           13           14           15           16
0.0031588557 0.0531045677 0.0001557366 0.0017106967 0.0955980736 0.0020314229 0.0341893200 0.0266185810 17           18
1.4447633641 0.0704056506

库克距离检验通过；

下再进行异方差和自相关检验，在进行异方差检验时，由于本问题为多元线性回归问题，残差 δ i \delta_i δi​ 可能随 任一个 x i x_i xi​ 变动：

cor.test(abs(lm.ssp$residual),sp$X1,method='spearman')

 Spearman's rank correlation rhodata:  abs(lm.ssp$residual) and sp$X1
S = 705.32, p-value = 0.2747
alternative hypothesis: true rho is not equal to 0
sample estimates:rho
0.2721168

说明误差关于解释变量 x 1 x_1 x1​ 是同方差的；

cor.test(abs(lm.ssp$residual),sp$X3,method='spearman')

 Spearman's rank correlation rhodata:  abs(lm.ssp$residual) and sp$X3
S = 980, p-value = 0.9672
alternative hypothesis: true rho is not equal to 0
sample estimates:rho
-0.01135191

说明误差关于解释变量 x 3 x_3 x3​ 是同方差的；

library(lmtest)
dwtest(lm.ssp)

 Durbin-Watson testdata:  lm.ssp
DW = 2.3043, p-value = 0.7162
alternative hypothesis: true autocorrelation is greater than 0

说明不同样本间不相关，自相关检验通过。

此时可以认为该模型通过检验，且最优线性回归方程为：
y = 41.4394 + 1.7374 ∗ x 1 + 0.1548 ∗ x 3 y=41.4394+1.7374*x_1+0.1548*x_3 y=41.4394+1.7374∗x1​+0.1548∗x3​

1.2 赛跑模型

1.3 时间序列模型

1.4 人口模型

1 指数增长模型

在假设年人口增长率为恒定常数的情况下，下式可用求出于经过 k k k 年的人口数量：
x k = x 0 ∗ ( 1 + r ) k x_k=x_0*(1+r)^k xk​=x0​∗(1+r)k
我们可由微分表达式解出相应的指数增长模型解：
{ d x d t ∗ x = r x ( 0 ) = x 0 ⇒ x ( t ) = x 0 e r t \begin{cases} \frac{dx}{dt*x}=r\\ x(0)=x_0\\ \end{cases} \Rightarrow x(t)=x_0e^{rt} {dt∗xdx​=rx(0)=x0​​⇒x(t)=x0​ert
需要对已知指数增长模型 x ( t ) = x 0 e r t x(t)=x_0e^{rt} x(t)=x0​ert 中的参数 r r r 进行参数估计:

参数 r r r 的估计

• 线性最小二乘法

  对指数增长模型 x ( t ) = x 0 e r t x(t)=x_0e^{rt} x(t)=x0​ert 两边取对数，我们有：
  l n ( x ( t ) ) = r t + l n x 0 ln(x(t))=rt+lnx_0 ln(x(t))=rt+lnx0​
  可以转化成线性函数的形式，容易看出， r r r 为这样一个线性函数的斜率：
  y = r t + a , ( y = ln ⁡ x ( t ) , a = ln ⁡ x 0 ) y=rt+a,(y=\ln{x(t)},a=\ln{x_0}) y=rt+a,(y=lnx(t),a=lnx0​)
  在直接进行线性最小二乘的情况下，我们只需要想办法找出一条直线使散点尽可能的分布在这条直线附近。直接去考虑lnx和t之间的线性函数，用 polyfit(x,y,degree)来求出斜率 r r r 的一个估计值，得到的结果分别是 r r r 和 a a a 的值，然后代入原模型求得解。

  %% 数据点
  x = [3.9 5.3 7.2 9.6 12.9 17.1 23.2    31.4 38.6 50.2 62.9 76.0 92.0 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4 308.7];
  x(end) = [];
  % t = 1790:10:2010;
  t = 0:length(x)-1;
  %%
  y = log(x); % 多项式拟合
  p = polyfit(t,y,1);
  %%
  % 绘图
  r = p(1); x0 = exp(p(2));
  f = @(t) x0*exp(r*t);
  tt = linspace(min(t), max(t), 101);% 绘图
  plot(t, x, 'o', tt, f(tt), '-')
  xlabel('年份-1790'); ylabel('人口数(百万)')
  legend('数据点','模型值', 'location', 'best')
  grid on
  title(sprintf('r=%f, x0=%f, sse : %f\n', r,x0,sum((x-f(t)).^2)))
  

• 数值微分法（p142表格中人口增长率来源）

  我们对每个不同阶段的人口数据作数值微分估计增长率，考虑最终取平均值得到的 r r r 更为可靠。

  容易知道，一般情况下，我们以 1 1 1 为步长， x ( t ) x(t) x(t) 在 t k t_k tk​ 处用数值微分（拉格朗日插值）中点法求得的一般斜率为
  x ′ ( t k ) = x ( t k + 1 ) − x ( t k − 1 ) 2 x'(t_k)=\frac{x(t_k+1)-x(t_k-1)}{2} x′(tk​)=2x(tk​+1)−x(tk​−1)​
  其中，对于端点有特殊形式：
  x ′ ( t 0 ) = − 3 x 0 + 4 x 1 − x 2 2 , x ′ ( t n ) = x n − 2 − 4 x n − 1 + 3 x n 2 x'(t_0)=\frac{-3x_0+4x_1-x_2}{2},x'(t_n)=\frac{x_{n-2}-4x_{n-1}+3x_n}{2} x′(t0​)=2−3x0​+4x1​−x2​​,x′(tn​)=2xn−2​−4xn−1​+3xn​​

  % 数据点
  x = [3.9 5.3 7.2 9.6 12.9 17.1 23.2    31.4 38.6 50.2 62.9 76.0 92.0 105.7 122.8 131.7 150.7 179.3 203.2 226.5 248.7 281.4 308.7];
  x(end) = [];
  t = 0:length(x)-1;xp = (x(3:end)-x(1:end-2)) / 2;
  xp = [x(1:3)*[-3;4;-1]/2, xp, x(end-2:end)*[1;-4;3]/2];
  rs = xp./x;
  r = mean(rs);tt = linspace(min(t), max(t), 101);
  plot(t, x, '-o', tt, x(1)*exp(r*tt),'-')
  title(sprintf('数值微分平均, r=%f', r))figure
  plot(t, rs, '-o')
  title('人口增长率(每10年)变化曲线')
  

• 改进的指数增长模型

  我们考虑修改原本 r r r 为常数的假设，由数据可以看出人口增长率随年份增长基本保持下降的趋势，我们作出新的假设为： r ( t ) = r 0 − r 1 t r(t)=r_0-r_1t r(t)=r0​−r1​t

  由上述过程，我们已经得到了一组 r ( t ) r(t) r(t) 的值：变量rs ，可以借此建立 t~rs 线性最小二乘求得 r 1 , r 0 r_1,r_0 r1​,r0​ 的值。

  此时的形式解为：
  x ( t ) = x 0 e ( r 0 − r 1 t 2 2 ) x(t)=x_0e^{(r_0-r_1\frac{t^2}{2})} x(t)=x0​e(r0​−r1​2t2​)

  %% 具体变量链接上一个代码块
  p = polyfit(t,rs,1);
  r0 = p(2); r1 = -p(1);
  f3 = @(t) x(1)*exp(r0*t-r1*t.^2/2);
  figure
  tt = linspace(min(t), max(t), 101);
  plot(t, x, 'o', tt, f3(tt) ,'-')
  legend('数据点', '模型计算结果', 'location', 'best')
  grid on
  title(['用改进的指数增长模型计算的美国人口, sse=' num2str(sum((x-f3(t)).^2))])
  

2 Logistic模型（阻滞增长模型）

在前面的基础上对问题进行深化。

人口增长到一定数量后增长率下降到主要原因是自然资源、环境条件等因素对人口增长起着阻滞作用，随着人口增加，阻滞作用越来越明显，logistic模型就说考虑到这些因素的基础上对指数增长模型的基本假设进行修改得到的。

首先，资源和环境对人口增长的阻滞作用主要体现在增长率 r r r 上， r r r 随人口数量 x x x 增长而下降，我们可以获得关系 r ( x ) r(x) r(x)，为简单，我们取线性函数 r ( x ) = a + b x r(x)=a+bx r(x)=a+bx，同时我们给予 a , b a,b a,b 实际含义：

• r r r 为 x = 0 x=0 x=0 时的理论增长率，即 r ( 0 ) = r r(0)=r r(0)=r，于是 a = r a=r a=r；

• x m x_m xm​ 为最大人口数量，当 x = x m x=x_m x=xm​ 后人口不再增长，也即有 r ( x m ) = r + b x m = 0 ⇒ b = − r x m r(x_m)=r+bx_m=0\Rightarrow b=-\frac{r}{x_m} r(xm​)=r+bxm​=0⇒b=−xm​r​ 于是有
  r ( x ) = r − r x m x = r ( 1 − x x m ) r(x)=r-\frac{r}{x_m}x=r(1-\frac{x}{x_m}) r(x)=r−xm​r​x=r(1−xm​x​)
  进一步有微分表达式：
  d x d t = r x ( 1 − x x m ) \frac{dx}{dt}=rx(1-\frac{x}{x_m}) dtdx​=rx(1−xm​x​)

我们可以得到最终的logistic模型：
x ( t ) = x m 1 + ( x m x 0 − 1 ) e − r t x(t)=\frac{x_m}{1+(\frac{x_m}{x_0}-1)e^{-rt}} x(t)=1+(x0​xm​​−1)e−rtxm​​

% 数据点
x = [3.9 5.3 7.2 9.6 12.9 17.1 23.2    31.4 38.6 50.2 62.9 76.0 92.0 105.7 122.8 131.7 150.7 179.3 203.2   226.5   248.7   281.4   308.7];
x(end) = [];
t = 0:length(x)-1;xp = (x(3:end)-x(1:end-2)) / 2;
xp = [x(1:3)*[-3;4;-1]/2, xp, x(end-2:end)*[1;-4;3]/2];
rs = xp./x;p = polyfit(x, rs, 1);
r0 = p(2); xm = -r0/p(1);
plot(x, rs, 'o', x, polyval([-r0/xm r0], x), '*')x0 = x(1);
funLogis1 = @(t) xm ./ (1+(xm/x0-1)*exp(-r0*t));
figure
tt = linspace(min(t), max(t), 101);
plot(t,x,'o',tt,funLogis1(tt),'*')
grid on
title(['用Logistic模型(线性拟合参数)计算的美国人口, sse=' num2str(sum((x-funLogis1(t)).^2))])

Chapter2 优化模型

2.1 线性规划

2.2 整数规划

2.3 0-1规划

Chapter3 评价模型

3.1 理想数法（TOPSIS）

基本步骤

Step1 【确定决策矩阵】 根据决策目标、备选方案、属性来确定决策矩阵

Step2 【决策矩阵标准化】 对决策矩阵进行标准化（a.区分效益型属性和费用型属性；b.标准化矩阵【归一化、模一化、标准化、最大化】）

Step3 【确定属性权重】 确定信息熵 E j E_{j} Ej​ 、区分度 F j F_{j} Fj​ 及属性权重 ω j \omega_{j} ωj​（各个属性对决策目标的重要程度）
信 息 熵 ： E j = − k ∑ i = 1 m r i j ∗ l n i j , k = 1 l n m , j = 1 , 2 , . . . , n 区 分 度 ： F j = 1 − E j , 0 ≤ F j ≤ 1 属 性 权 重 ： w j = F j ∑ j = 1 n F j , j = 1 , 2 , . . . , n 信息熵：E_{j}=-k\sum_{i=1}^{m}r_{ij}*ln_{ij},k=\frac{1}{ln_{m}},j=1,2,...,n\\ 区分度：F_{j}=1-E_{j},0\leq F_{j}\leq 1\\ 属性权重：w_{j}=\frac{F_{j}}{\sum_{j=1}^{n}F_{j}},j=1,2,...,n 信息熵：Ej​=−ki=1∑m​rij​∗lnij​,k=lnm​1​,j=1,2,...,n区分度：Fj​=1−Ej​,0≤Fj​≤1属性权重：wj​=∑j=1n​Fj​Fj​​,j=1,2,...,n

Step4 【综合评价】 求出加权决策矩阵的正负理想解，再求得各方案对正理想解的相对接近度 C j + ， 0 < F j < 1 C^{+}_{j}，0\lt F_{j}\lt 1 Cj+​，0<Fj​<1
C j + = S − S + + S − C^{+}_{j}=\frac{S^{-}}{S^{+}+S^{-}} Cj+​=S++S−S−​

案例分析

#理想数法Topsis# 输入权重矩阵，数值已经统一为效益型
m=matrix(c(1/25,1/18,1/12,9,7,5,7,7,5),ncol=3)
# 构造函数以便模一化
fun1<-function(x) {sqrt(sum(x^2))}
# 模一化矩阵
m.1<-sweep(m,2,apply(m,2,fun1),"/")fun2=function(x) (-1/log(3))*sum(x*log(x))
E<-apply(m.1,2,fun2) # 计算信息熵(对列定义)
w<-(1-E)/sum(1-E) # 计算权重
m.11<-sweep(m.1,2,w,"*") # 权重矩阵加权(对列定义)a<-NULL
b<-NULL
#构造最优方案与最劣方案
for (i in 1:3) a[i]=max(m.11[,i])
for (i in 1:3) b[i]=min(m.11[,i])
fun3=function(x,y) sqrt(sum((x-y)^2))
s2<-apply(m.11,1,fun3,y=b) # 负理想解
s1<-apply(m.11,1,fun3,y=a) # 正理想解
#计算方案见距离
c=s2/(s2+s1)
#各方案的评价得分
c

3.2 层次分析法（AHP）

层次分析法基本步骤

Step1 【建立层次结构模型】 将目标问题分解为目标层、准则层、方案层

Step2 【构建成对比较矩阵】

• 准则层可有子层，其中第一层的指标数： 2 ≤ n ≤ 9 2\leq n\leq 9 2≤n≤9

• 某一层占上一层的【权重个数】取决于该层的指标数（例如，准则层占目标层的权重数为4，取决于准则层的四个指标）

• 求某一层占上一层的权重时是在求下层占上一层某一个属性的比重，于是所需要的【矩阵个数】取决所需构建的成对比较矩阵的数量取决于上一层的指标数（例如，求方案层占准则层的权重时需要4个矩阵，是因为该层有四个属性指标，会求得4个权重向量）

  1. 准则层对目标层的成对比较矩阵
  2. 方案层对准则层的成对比较矩阵

Step3 【计算每层权重】

1. 代数理论

2. 简化手动算法

   列归一->行求和->行归一

Step4 【计算综合得分】

• 矩阵相乘法
• 连接线法，对于每一个方案 A i A_i Ai​，取其连线和各自权重相乘相加得到综合评分

案例分析

#层次分析AHP
m=matrix(c(1,1/2,1/3,2,1,1/2,3,2,1),ncol=3) #输入成对比较矩阵
eig<-eigen(m)  #计算特征值与特征向量
ma<-max(Re(eig$val)) #最大特征根
n=3 #矩阵维数
R=c(0,0,0.52,0.89,1.12,1.26,1.36,1.41,1.46) #临界表，固定的
ci=(ma-3)/(n-1)
cr=ci/R[n]  #一致性检验
if (cr<0.1) print("一致性检验通过") else print("错误，请重设比较矩阵")w<-eig$vec[,which.max(Re(eig$val))]
w1<-as.matrix(Re(w/sum(w)),ncol=3)    #计算权重并化为矩阵m1=matrix(c(1/25,1/18,1/12,9,7,5,7,7,5),ncol=3) #输入权重矩阵，数值已经统一为效益型
a=m1%*%w1 #计算方案得分

3.3 决策树

3.4 排名问题

原理：

设 A A A 为得分矩阵，假定矩阵 A A A 有特征值 λ \lambda λ 和特征向量 S S S ，行和即为总分，相当于 S 1 = A ∗ e ⃗ S_1=A*\vec{e} S1​=A∗e ，同理可以有 k k k 次加权，进而有：
{ S k = A k ∗ e ⃗ l i m k → ∞ A k ∗ e ⃗ λ k = S ⇒ l i m k → ∞ S k λ k = S \begin{cases} S_k=A^k*\vec{e}\\ lim_{k\rightarrow\infty}\frac{A^k*\vec{e}}{\lambda^k}=S \end{cases} \space\Rightarrow lim_{k\rightarrow\infty}\frac{S_k}{\lambda^k}=S {Sk​=Ak∗e limk→∞​λkAk∗e ​=S​ ⇒limk→∞​λkSk​​=S
容易看出当 k → ∞ k\rightarrow\infty k→∞ 时， S k S_k Sk​ 将趋向 λ k S \lambda^kS λkS，而 λ k \lambda^k λk 并不影响向量中各分量的相对大小，于是可以用 A A A 的最大特征根 λ \lambda λ 对应的特征向量 S S S 等价视为 S k S_k Sk​，也即用得分矩阵的特征向量 S S S 来作为排名的依据

基本步骤

Step1 【列出分数矩阵】

Step2【求最大特征值对应的特征向量】

• 直接算法

  1. 人工极限逼近法

     myfunc<-function(x,n){i<-1A<-xwhile(i<n){A<-A%*%Ai<-i+1}A
     }
     #取一个较大的n，人工调试到值接近不变，求出的A则是S的一个近似
     myfunc(A,100)
     

  2. 直接求对应特征向量，利用 R R R 或 M A T L A B MATLAB MATLAB 中的 eigen(A) 函数

     eigen(A)
     

• 简化手动算法（列归一；行求和；（取转置）；（行归一））

3.5 聚类分析

系统聚类

• 首先要去明确该算法的原理做法和目的：系统聚类法的做法是开始时把每个样品作为一类，然后把最靠近的样品（即距离最小的群品）首先聚为小类，再将已聚合的小类按其类间距离再合并，不断继续下去，最后把一切子类都聚合到一个大类

总体步骤

Step1 【最短距离计算】

• 我们使用 dist(wine) 函数进行两两样本间的距离求取

Step2 【系统聚类】

• 我们利用 hclust(d) 函数来进行系统聚类
• 此时可利用 plot(hclust(d)) 来给出一个聚类的可视化图像

Step3 【结果矩阵展现】

• 若有输出结果矩阵如下所示：

        [,1] [,2]      [,3]
  [1,]   -7   -8  105.1565
  [2,]   -4   -5  138.1701
  [3,]  -28  -30  189.6000
  [4,]   -3    1  194.3796
  [5,]  -15  -31  196.4032
  [6,]  -10  -23  234.9086
  [7,]   -6  -29  235.2408
  [8,]  -27    3  244.4392
  

  我们对该结果矩阵做以下解读：

  第一行：将第7个观测点和第8个观测点聚为一类
  第二行：将第4个观测点和第5个观测点聚为一类
  …
  第四行：将第3个观测点聚到第一行聚的那一类

  …

案例分析

数据源：葡萄酒分类问题.csv

解答：

wine<-read.csv(file.choose())
d<-dist(wine)
hclust(d)->HC
plot(HC)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VYkxyDTm-1624167571105)(F:\University\课程整合\大二下\数学模型\数学模型.assets\image-20210429201254808.png)]

cbind(wine.hc$merge,wine.hc$height)

      [,1] [,2]      [,3][1,]   -6  -10 0.4835287[2,]   -3   -4 0.5280152[3,]   -7   -8 0.5444263[4,]   -5   -9 0.5692100[5,]   -1    2 0.8718371[6,]   -2    4 1.1131936[7,]    1    3 2.3153617[8,]    6    7 4.5578504[9,]    5    8 9.3714140

我们对上述的结果给出解读：

将元素 6 和 10 归为Ⅰ类；将元素 3 和 4 归为Ⅱ类； 将元素 7 和 8 归为Ⅲ类； 将元素 5 和 9 归为Ⅳ类； 将元素 1 归入Ⅱ类；

将元素 2 归入Ⅳ类；将 Ⅰ 和 Ⅲ 归为Ⅰ类；将 Ⅳ 和 Ⅰ 归为Ⅰ类；将 Ⅱ 和 Ⅰ 归为Ⅰ类

于是我们最终可以得出以下分类，笔者按照 h e i g h t = 1.5 height=1.5 height=1.5 来分界：

类别	分类
Ⅰ类	6，10
Ⅱ类	1，3，4
Ⅲ类	7，8
Ⅳ类	2，5，9

---

KMS聚类

随机聚类法

总体步骤

Step1 【数据读入】

• 对首列是文字的csv数据文件读入时采用read.table(file.choose())

Step2 【KMS聚类】

• kmeans(eco,4)函数进行KMS聚类

  参数解析：eco 为待聚变量； 4 4 4 为分类的类别数

Step3 【排名分类可视化】

• sort(kmeans(eco,4)$clust) 函数进行排名

案例分析

数据源：中国各省经济发展水平.csv

解答：

fdc<-read.table("clipboard")

kmeans(fdc,5)->km

sort(km$cluster)

内蒙古   江西   湖南   贵州   云南   宁夏   新疆   河北   辽宁 黑龙江   江苏   湖北   陕西   甘肃   山西 1      1      1      1      1      1      1      2      2      2      2      2      2      2      3 吉林   安徽   山东   河南   广西   重庆   四川   青海   北京   上海   天津   浙江   福建   广东   海南 3      3      3      3      3      3      3      3      4      4      5      5      5      5      5

由结果分析：

KMS聚类分析结果，根据参数可将以上30个城市的房地产情况分为5类：

类别	城市
Ⅰ	内蒙古，江西，湖南，贵州，云南，宁夏，新疆
Ⅱ	河北，辽宁，黑龙江，江苏，湖北，陕西，甘肃
Ⅲ	山西，吉林，安徽，山东，河南，广西，重庆，四川，青海
Ⅳ	北京，上海
Ⅴ	天津，浙江，福建，广东，海南

---

主成分分析（PCA）

参考博客：https://blog.csdn.net/zhongkejingwang/article/details/42264479

来源：综合评价问题指标量多（存在相关性）

目标：利用相关性，去除重叠信息，减少指标量，保留主要信息，进行数据降维处理

本质：利用数据本身的变异系数去筛选变量

基本步骤

Step1 【数据读入】

Step2 【函数求主成分】

• princomp(data,cor=T) 主成分函数

Step3 【累加每一行的成分然后进行排名】

• rowSum(data$scores) 函数对每一行的成分进行排名，

Step4 【排名可视化（数据框形式展现）】

• 注意，要总数量+1-排名数才对，因为正负的原因比如原本的-4，-2，1，8，的第三个应该是排名第二的，此时需要5-3=2，5=4+1，3为原排名数

案例分析

数据源：工业部门各行业指标_主成分

解答：

ind<-read.csv(file.choose())
hy<-ind[,1]
ind<-ind[,-1]

target<-princomp(ind,cor=T)

结果概览：

summary(target,loadings=T)

Importance of components:Comp.1    Comp.2    Comp.3     Comp.4     Comp.5     Comp.6      Comp.7
Standard deviation     1.7620762 1.7021873 0.9644768 0.80132532 0.55143824 0.29427497 0.179400062
Proportion of Variance 0.3881141 0.3621802 0.1162769 0.08026528 0.03801052 0.01082472 0.004023048
Cumulative Proportion  0.3881141 0.7502943 0.8665712 0.94683649 0.98484701 0.99567173 0.999694778Comp.8
Standard deviation     0.0494143207
Proportion of Variance 0.0003052219
Cumulative Proportion  1.0000000000Loadings:Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
X1  0.477  0.296  0.104         0.184         0.758  0.245
X2  0.473  0.278  0.163 -0.174 -0.305        -0.518  0.527
X3  0.424  0.378  0.156                      -0.174 -0.781
X4 -0.213  0.451         0.516  0.539 -0.288 -0.249  0.220
X5 -0.388  0.331  0.321 -0.199 -0.450 -0.582  0.233
X6 -0.352  0.403  0.145  0.279 -0.317  0.714
X7  0.215 -0.377  0.140  0.758 -0.418 -0.194
X8         0.273 -0.891        -0.322 -0.122

screeplot(target,type="line")

Y<-rowSums(target$scores)
Y_rk<-14-rank(Y)
target_last<-data.frame(Y,Y_rk,hy)

            Y Y_rk   hy
1   4.9096093    2 冶金
2  -1.2532354    7 电力
3  -2.0708453   11 煤炭
4   2.3634281    3 化学
5   6.1667330    1 机器
6  -1.9458162    9 建材
7  -1.5359765    8 森工
8   1.0951279    4 食品
9   0.3171423    5 纺织
10 -2.0439633   10 缝纫
11 -2.9789630   12 皮革
12  0.1246927    6 造纸
13 -3.1479336   13 文教

最终该工业部门各行业的指标排名如上结果所示

附录

最小二乘法

参考回答：https://www.zhihu.com/question/37031188

---

---

1. 为行文通顺，文中没有具体案例提到的所有函数中的参数都只是代指 ↩︎

2. 没有明确指出的情况下，默认显著性水平为 0.05 ↩︎

【数学模型】模型整合相关推荐

1. 计划行为理论和技术接受模型整合模型图形_音乐与语言加工的二元模型（dual modal）...

   随着科技研究手段的进步,随着心理学主导的理论从强调外在的行为向强调内在的认知转变,音乐与语言之间关系的对比探讨,成为神经心理学研究的热点,而且业已成为经典的研究范式,吸引了语言学.音乐学.神经科学等跨 ...

2. 模型整合之模型堆叠——详细理解Stacking model

   详细理解Stacking model 如果你得到了10个不一样的model,并且每个model都各有千秋,这个时候你该怎么选?想必你一定是很为难吧,但通过集成方法,你可以轻松的将10个model合成为 ...

3. java 3d模型插件_3D模型整合插件 Kitbasher V1.2 支持3DS MAX 2012~2018

   Kitbasher V1.2 插件是可以在3DS MAX中将多个3D模型组合成一个完整的模型,多种参数可以设置. This plugin for 3ds Max gives users the pow ...

4. 无需多个模型也能实现知识整合？港中文MMLab提出「烘焙」算法，全面提升ImageNet性能...

   视学算法专栏 转载自:机器之心 作者:葛艺潇 来自港中文 MMLab 的研究者提出一种烘焙(BAKE)算法,为知识蒸馏中的知识整合提供了一个全新的思路,打破了固有的多模型整合的样式,创新地提出并尝试了 ...

5. TensorFlow与PyTorch模型部署性能比较

   TensorFlow与PyTorch模型部署性能比较 前言 2022了,选 PyTorch 还是 TensorFlow?之前有一种说法:TensorFlow 适合业界,PyTorch 适合学界.这种说 ...

6. AI 框架部署方案之模型部署概述

   0 概述 模型训练重点关注的是如何通过训练策略来得到一个性能更好的模型,其过程似乎包含着各种"玄学",被戏称为"炼丹".整个流程包含从训练样本的获取(包括数据采 ...

7. 一文读懂目标检测模型（附论文资源）

   来源: 大数据文摘 本文共1443字,建议阅读5分钟. 本文为你详细介绍目标检测,并分享资源大礼包,为你的目标检测入门打下基础. 这是一份详细介绍了目标检测的相关经典论文.学习笔记.和代码示例的清单, ...

8. “烘焙”ImageNet：自蒸馏下的知识整合

   ©作者|葛艺潇 学校|香港中文大学博士生 研究方向|图像检索.图像生成等 最先进的知识蒸馏算法发现整合多个模型可以生成更准确的训练监督,但需要以额外的模型参数及明显增加的计算成本为代价.为此,我们提出 ...

9. sam服务器是什么_使用SAM CLI将机器学习模型部署到无服务器后端

   sam服务器是什么 介绍 (Introduction) Over the last year at CMRA, we have been incorporating more machine lear ...

最新文章

1. 关于C语言中printf函数“输出歧视”的问题
2. mybatisplus 结果_springboot整合mybatisPlus 乐观锁的实现
3. python propresql mysql_Python中操作mysql的pymysql模块详解
4. ptmalloc内存分配释放
5. ABAP里的OAuth2.0 Standard Package
6. 尾调用优化 java_为什么JVM仍然不支持尾调用优化？
7. 114. 二叉树展开为链表 golang
8. matlab电压稳定极限,电力系统电压稳定性的Matlab建模分析
9. python ftp timeout_python - FTP文件传输期间Python数据通道超时 - 堆栈内存溢出
10. c语言中参数的传递方式是,C语言函数的参数及传递方式
11. 步步为营！高手教你如何有效使用深度学习解决实际问题
12. 如何在钉钉上开发自己的应用_对企业来说无代码开发平台是否安全
13. oracle 常用语句汇总
14. Atitit 延迟绑定架构法attilax总结
15. mysql表不存在但实际存在_历史上有哪些实际上并不存在的人物但很多人相信他存在的？...
16. python制作网页挂机_python使用ip代理抓取网页
17. 加速Web开发的9款知名HTML5框架
18. 微信小程序下拉刷新、下拉加载下一页操作逻辑
19. 需求分析之矩阵分析法
20. 1941. Scary Martian Word

热门文章

1. Billu_b0x靶机
2. FPGA：三大协议（IIC、UART、SPI）之IIC
3. 普安特：猫咪在什么状况下容易患湿疹？
4. 你的电脑该如何选择？-涵子的个人想法
5. 剪切板经常出问题，无法复制等现象。
6. 【大数据数仓项目集群配置 一】
7. 使用级联有序调度窗口生成无序调度
8. PageHelper版本差异造成的Interceptor和dialect问题。
9. 网安云新品速递 | 渗透测试服务，助力企业业务安全发展
10. 【Busybox】Busybox源码分析-01 | 源码目录结构和程序入口