文章目录
个人在数学建模过程中的学习梳理的笔记,编程软件和语言以RStudio/R为主。
供学习参考。部分内容待补充。
Chapter1 预测模型
1.1 线性回归模型
总体步骤
注1
step1【 建立模型】
- 使用内置函数
lm(Y~X,data=bp)
函数
step2【 F、t、拟合优度检验】
- 在建立模型后,我们需要对该模型进行检验,该检验包括:从整体效果上来看,该模型是否能够接受?从每个解释变量上来看,是否必要?同时,此模型的拟合效果又如何?
- 为讨论上述问题,需要使用
summary(lm.bp)
函数对模型进行以下检验:- 模型整体效果评价(回归方程显著性检验):F检验( H 0 : β 1 = β 2 = . . . = 0 H_0:\beta_1=\beta_2=...=0 H0:β1=β2=...=0): p < α p<\alpha p<α 时检验通过;
- 解释变量个体必要评价(回归系数显著性检验):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<α 时检验通过;
- 拟合优度检验: 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 【模型优化】
-
模型的优化一般有以下几种情况:
- 在进行summary检验时,三类检验发现效果不好,可以考虑剔除某个变量(谨慎);
- 模型未通过异方差和自相关检验(即残差不是同方差或残差自相关);
- 某些变量可能通过了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λ−1ln(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 K2CO3 溶液并受溴化物水解的有机磷浓度 x 2 x_2 x2 以及土壤内溶于 K 2 C O 3 K_2CO_3 K2CO3 溶液但不溶于溴化物的有机磷 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 ‘ ’ 1
Residual 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
若显著性水平2
α
=
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 ‘ ’ 1
Residual 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 rho
data: 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 rho
data: 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 test
data: 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)=x0ert
需要对已知指数增长模型
x
(
t
)
=
x
0
e
r
t
x(t)=x_0e^{rt}
x(t)=x0ert 中的参数
r
r
r 进行参数估计:
参数 r r r 的估计
-
线性最小二乘法
对指数增长模型 x ( t ) = x 0 e r t x(t)=x_0e^{rt} x(t)=x0ert 两边取对数,我们有:
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−r1t
由上述过程,我们已经得到了一组 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)=x0e(r0−r12t2)%% 具体变量链接上一个代码块 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=−xmr 于是有
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−xmrx=r(1−xmx)
进一步有微分表达式:
d x d t = r x ( 1 − x x m ) \frac{dx}{dt}=rx(1-\frac{x}{x_m}) dtdx=rx(1−xmx)
我们可以得到最终的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+(x0xm−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∑mrij∗lnij,k=lnm1,j=1,2,...,n区分度:Fj=1−Ej,0≤Fj≤1属性权重:wj=∑j=1nFjFj,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个权重向量)
- 准则层对目标层的成对比较矩阵
- 方案层对准则层的成对比较矩阵
Step3 【计算每层权重】
-
代数理论
-
简化手动算法
列归一->行求和->行归一
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∗elimk→∞λ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【求最大特征值对应的特征向量】
-
直接算法
-
人工极限逼近法
myfunc<-function(x,n){ i<-1 A<-x while(i<n){ A<-A%*%A i<-i+1 } A } #取一个较大的n,人工调试到值接近不变,求出的A则是S的一个近似 myfunc(A,100)
-
直接求对应特征向量,利用 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.999694778
Comp.8
Standard deviation 0.0494143207
Proportion of Variance 0.0003052219
Cumulative Proportion 1.0000000000
Loadings:
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