数学建模学习笔记_数学建模笔记 csdn

数学建模学习笔记

现学现卖，随缘更新QwQ
主要根据b站大师兄的视频整理而成，有不懂的可以去看原视频

---

一、 层次分析法

1.1 矩阵的一致性及其检验

对于矩阵$A$，$a_{ij}$表示第i个因素的重要程度是第j个因素的A[i][j]倍。
若矩阵A满足：$a_{ij}$ * $a_{jk}$ = $a_{ik}$，则称A为 一致矩阵，否则称为 不一致矩阵
对于一致矩阵A，有如下 三条性质：

1. $rank(A)$ = 1，且唯一非零特征根为n
2. A的任意列向量都是对于特征根n的特征向量

上述性质告诉我们，只要满足两两行/列 成倍数关系，就是一致性矩阵。
对于不一致矩阵，其最大特征根 ${\lambda }_{max}$>n，$\lambda$ 与n相差越大，其不一致程度越大。由此，给出一致性检验指标$CI$ = $\frac{{\lambda }_{max}-n}{n-1}$ 其中 ${\lambda }_{max}$表示最大特征值。CI越大，表示一致性越差，越接近0则一致性越好。
由于 $CI$受随机因素影响，阶数n越大，该因素越显著。为了衡量$CI$的大小，引入$RI$表示随机一致性指标，并以 一致性比例 $CR$=$\frac{CI}{RI}$ 作为判断依据。若$CR$<0.1，则认为判断矩阵的一致性可以接受，否则需要对判断矩阵进行修正，若不满足，则需构造出更多的倍数关系，直至$CR$ < 0.1

1.2 权重计算

层次分析法的最终目的就是获得不同要素的权重，以便于后续进行加权平均。具体地，先将每一列进行 归一化，也即将$a_{ij}$替换为$\frac{a_{ij}}{{\sum }_{i=1}^n{a}_{ij}}$；
然后，对于非一致矩阵，将$W_i=\frac{{\sum }_{j=1}^n{a}_{ij}}{n}$作为第i个要素的权重。这种求权重的方法称为算数平均求权重。由于一致矩阵各列成比例，所以只需取第一列即可，不需要进行算术平均。
另一种求权重的方式为：将最大特征值对应的特征向量进行归一化，第i个分量的值即为$W_i$。这种方法称为特征值法求权重。
这种方法较于前一种方法更加方便使用，因而在实际比赛和科研中建议使用特征值法。

1.3 具体流程

1. 根据专家意见、问卷调查等方式获得判断矩阵$A$，并对$A$进行一致性检验，若$CR$<0.1，认为一致性检验通过，否则构造新的判断矩阵
2. 将判断矩阵各个列进行归一化，随后对每一行求算数平均，所得即为权重

层次分析法是一种主观确定权重的方式，后续还会有客观确定权重的方式——熵权法

二、模糊综合评测

一般的对立集合具有排中律，即“非此即彼”，而对于不满足排中律的情况，则需要引入隶属函数衡量每个元素对于不同集合的隶属度。换句话说，就是对于一个论域 $U$ 构造一个函数$f:U\to [0,1]$，以评判$U$中每一个元素的隶属度。与层次分析法一样，这种方法也属于主观评价方法。

2.1 隶属函数

对于模糊集合$A=“年轻”，U=(0,120)$，定义隶属函数${\mu }_A$：
μ A = { 1 ， 0 < x < 20 40 − x 20 ， 20 ≤ x ≤ 40 0 ， 40 < x < 120 \mu_A =

$$\begin{cases} 1，0<x<20\\ \displaystyle\frac{40-x}{20}，20≤x≤40\\ 0，40<x<120 \end{cases}$$ μA​=⎩ ⎨ ⎧​1，0<x<202040−x​，20≤x≤400，40<x<120​
可以看到，这个函数对U中每一个元素都给出了一个对于$A$的隶属度，越大则越符合。
模糊集合的表示方法包括扎德表示法，序偶表示法，向量表示法等。
序偶表示法：$\{(x_1,A(x_1)),(x_2,A(x_2))...(x_n,A(x_n)\}$
向量表示法： A = { A ( x 1 ) , A ( x 2 ) . . . A ( x n ) } A=\{A(x_1),A(x_2) ... A(x_n)\} A={A(x1​),A(x2​)...A(xn​)}
当U为无限集时，定义$A={\int }_{x\in U}\frac{{\mu }_A(x)}{x}dx$

上面例子中$A=“年轻”$属于极小型，因为值越小隶属度应当越高。对于这类集合，其隶属函数应形如：
μ A = { 1 ， x < a b − x b − a ， a ≤ x ≤ b 0 ， x > b \mu_A =

$$\begin{cases} 1，x<a\\ \displaystyle\frac{b-x}{b-a}，a≤x≤b\\ 0，x>b \end{cases}$$ μA​=⎩ ⎨ ⎧​1，x<ab−ab−x​，a≤x≤b0，x>b​
相应地，对于极大型 ，其形式为：
$${\mu }_A=\{\begin{array}{l}0，x<a\\ \frac{x-b}{b-a}，a\le x\le b\\ 1，x>b\end{array}$$
对于中间型，其左侧为极大型函数，右侧为极小型函数，不作赘述。
上面的构造方法称为梯形型，事实上还存在k次抛物型、柯西型、正态型等等，在需要时可以自行查阅。

2.2 隶属函数的确定方法

确定隶属函数有模糊统计、F分布、三分法等，至少应该满足如下条件：

1. 极小型集合的下界a，小于a的元素均不属于其他集合（极大型同理）
2. 增长趋势应符合主观经验

2.3 模糊综合评测

对于一级模糊综合评价，遵循如下步骤：

1. 确定因素集$U$，例如 {工作业绩、工作态度、沟通能力}，评语集$V$，例如 {好，较好，中，差，很差}
2. 确定各因素（$U$中各个元素）的权重，若无数据可采用层次分析法，若给出数据则使用熵权法的TOPSiS，也可以不确定权重
3. 对每个$u_i$，确定其对于每个$v_j$的隶属度，对指标$u_i$的评判记作： $R=[r_{i1},r_{i2}\dots r_{in}]$ ，其中$r_{ij}$表示$u_i$对于$v_j$的隶属度。
4. 对每一列进行加权平均（即左乘行向量）,取数值最大的评语作为最后综合评判结果

若将评语集更改为方案集，上述流程可以判断出哪个方案最优！ 此时也可称评语集不带有评价色彩。
类似地，对于因素集中、指标过多时，可根据相关性将因素归纳成一个个小集合，从而进行模糊综合评价的嵌套，也即多级模糊综合评测。

三、熵权法

信息熵是衡量混乱程度的量，其定义为:
$$H(X)=\sum _{i=1}^n[p(x_i)I(x_i)]=-\sum _{i=1}^n[p(x_i)ln(p(x_i))]$$
其中$I(x_i)=-ln(p(x_i))$表示$x_i$的信息量。根据公式，可以将信息熵理解为对信息量的期望值。信息熵越大，所当前状态的混乱程度最大，也即掌握的信息最少（信息的本质就是减少混乱程度，增大确定性），以该指标衡量对象的可靠性越差，信息有效性也就越小。

3.1 正向化与标准化

在应用熵权法之前，要先对数据进行正向化处理。指标大致可分为四类：极小型（越小越好）、极大型（越大越好）、中间型（越接近越好），区间型（落在区间内最好）。我们的目标是将其余三种类型转换为极大型。

• 将极小型转化为极大型，公式为$max-x$或$\frac{1}{x}(x>0)$，$max$表示该指标下的最大值；
• 中间型转化为极大型，公式为${\tilde{x}}_i=1-\frac{\left|x_i-x_{best}\right|}{M}$，其中 M = m a x { ∣ x i − x b e s t ∣ } M=max\{\left|x_i-x_{best}\right| \} M=max{∣xi​−xbest​∣}表示指标下的最大偏差；
• 区间型转化为极大型，方式类似于中间型隶属函数，${\tilde{x}}_i=\{\begin{array}{l}1-\frac{a-x_i}{M}，x_i<a_1\\ 1，a\le x_i\le b\\ 1-\frac{x_i-b}{M}，x_i\ge b\end{array}$，其中 M = m a x { a − m i n , m a x − b } M=max\{a-min,max-b\} M=max{a−min,max−b}为最值到边界的最大距离。
• 对于极大型，若存在负元素，则可将每一个元素替换为$x-min$

这样，三类数据都可以用极大型指标的方式呈现，且均非负。
随后，为了消除量纲的影响，对每个数据进行标准化。设标准化后的矩阵为$Z$，则$Z$满足：$z_{ij}=\frac{x_{ij}}{\sqrt{{\sum }_{i=1}^n{x}_{ij}^2}}$
对每一列再做归一化，得到的矩阵$P$称为比重矩阵

3.2 流程

1. 对数据进行正向话、标准化与归一化，最终得到矩阵$P$
2. 令信息熵$e_j=-\frac{1}{\ln n}{p}_{ij}\ln (p_{ij})(j=1,2,\dots ,m)$，注意这里新增了$-\frac{1}{\ln n}$，以将结果规范到[0,1]之间
3. 由此得到信息效用值$d_j=1-e_j$，再将信息效用值归一化，得到熵权 $W_j=d_j/{\sum }_{j=1}^m{d}_j(j=1,2,\dots ,m)$

四、TOPSIS优劣解距离法

该方法用于衡量某一个方案的优劣，具体思想是将每个衡量对象的指标标准化后分析其与最优解与最劣解的距离，进而得出最优解。

1. 将数据化为$n\times m$的矩阵$A$，n表示待评测对象个数，m表示指标个数。将$A$进行正向化、标准化；
2. 选出各个指标下最大值$M_i$与最小值$m_i$（也即每一列最大值），组成最优解向量$M$与最劣解向量$m$；
3. 对每个待衡量对象（也即$A$的每个行向量）计算其与$M$与$m$的欧氏距离，分别记作$D_i^{+}$与$D_i^{-}$，定义该对象的得分$S_i=\frac{D_i^{-}}{D_i^{+}+D_i^{-}}$。$S_i$越大，说明该对象越优，反之则越劣。

如果要考虑权重，则需将距离计算公式改为$\sqrt{{\sum }_{i=1}^m{w}_i(x_i-y_i{)}^2}$的形式即可，权重的获取方法如上述。

五、灰色关联分析

灰色关联分析用于处理已知信息较少的问题，通过分析大致趋势的贴合程度得出各个影响因素的重要程度。

5.1 前置定义

• 参考序列（母序列）：能反映系统行为特征的数据序列，记作$x_0$，类似于“因变量”。
• 比较序列（子序列）：影响系统行为的因素组成的数据序列，记作$x_i(i=1,2,\dots n)$，类似于“自变量”。
• 两级最小差$a=mi{n}_{s,t}\left|x_0(t)-x_s(t)\right|$
• 两级最大差$b=ma{x}_{s,t}\left|x_0(t)-x_s(t)\right|$

在这里，我们引入新的去除量纲的方法，即对每个指标中的元素除以该指标的均值。

5.2 流程

1. 对数据进行预处理（正向化，去除量纲等等），并计算出两级最大/小差；
2. 对每个元素定义$\gamma (x_0(k),x_i(k))=\frac{a+\rho b}{\left|x_0(k)-x_i(k)\right|+\rho b}$，其中$\rho$称为分辨系数，一般取值为0.5；
3. 将各指标曲平均值作为该序列的灰色关联度，即$\gamma (x_0,x_i)=\frac{1}{n}\sum _{k=1}^n\gamma (x_0(k),x_i(k))$，值越大，表明该因素与系统行为关联越大。

具体操作中有一些excel的小技巧，详见 大师兄的讲解视频

六、线性规划

数学规划是一类针对有限制的关系寻求最优解的方法。若关系及其限制均为线性，则称这类规划为线性规划。

6.1 基本定义

• 决策变量：即对最优解有影响的因素
• 目标函数：形如$y=a_1{x}_1+a_2{x}_2+\dots a_n{x}_n$，规划目标即为求得y的最大/小值
• 约束条件：形如$b_{p1}{x}_{p1}+b_{p2}{x}_{p2}+\dots b_{pn}{x}_{pn}\le (或\ge )c$，表示对决策变量的限制条件

针对这类问题，一般采取拉格朗日乘数法，由于实际应用中可以直接调用函数求解，故不作介绍。

6.2 标准型

调用MATLAB要求对输入数据进行标准化。引入线性规划的标准型：
min ⁡ f T x , s . t . { A ⋅ x ≤ b A e q ⋅ x = b e q l b ≤ x ≤ u b \min f^T x,s.t.

$$\begin{cases} A \cdot x\leq b\\ Aeq \cdot x = beq\\ lb \leq x \leq ub \end{cases}$$ minfTx,s.t.⎩ ⎨ ⎧​A⋅x≤bAeq⋅x=beqlb≤x≤ub​
其中$f^T$为系数向量$[c_1,c_2,\cdots c_n]$
在MATLAB中，求解函数格式为：

[x,val] = linprog(f, A, b. Aeq, beq, lb, ub)
% 若某项约束不存在则在对应位置填"[]"
1
2

其中 x 为目标函数取得最小值时各个变量的取值，val 表示目标函数的最小值
由于最开始给定的数据往往不满足上述形式，所以往往需要先进行预处理。

在实际建模中，我们也常常能见到多个目标函数的情况，例如，投资问题中，我们往往希望风险尽可能低，收益尽可能高。这时，我们可以定义它们的加权平均为新的目标函数。也可以通过将其余目标化为限制条件，进而保证目标函数只有一个，比如认为只要风险小于一定程度即认为可行，将收益最大化。

6.3 整数线性规划

整数线性规划，即部分或全部决策因素只能为整数的一类问题，在MATLAB中也有直接对应的函数，即：

[x,val] = intlinprog(f, intcon,  A, b, Aeq, beq, lb, ub)
1

注意到，intlinprog 相较于 linprog 新增了 intcon 这一参数，用来指定哪些变量为整数，例如决策变量$x_1,x_2,x_3$中第一个和第三个为整数，则 intcon=[1,3]

这里引入一个特殊的整数规划模型——0-1整数规划，即变量只能取0或1。落实到函数中，只需要令上述函数中lb=0,ub=1即可。
下面给出对应例题，依次为 0-1背包问题 与 指派问题：


对于整数规划，一个常见误解是只要对线性规划所得解进行简单取整就可以了。事实上，整数线性规划虽然与一般线性规划有许多联系，但绝不能粗暴地以取整相联系。
不过，根据这个思路，我们可以在一般线性规划的最优解附近寻找整数可行解，这种方法称为 完全枚举法；也可以通过将最优解中非整数元素$x_i$作为对整数域的分割点，分别对 { y i ∈ Z ∣ y i ≥ ⌊ x i ⌋ + 1 } \{ y_i \in \mathbb{Z}|y_i \geq \lfloor x_i\rfloor + 1 \} {yi​∈Z∣yi​≥⌊xi​⌋+1} 和 { y i ∈ Z ∣ y i ≤ ⌊ x i ⌋ } \{ y_i \in \mathbb{Z}|y_i \leq \lfloor x_i\rfloor \} {yi​∈Z∣yi​≤⌊xi​⌋}两个集合进行讨论，直到所有领域都讨论完 或 找到整数最优解为止。这种方法称为分枝定界法。

6.4 一点补充

对于 非线性规划，一些典型的处理方法包括 蒙特卡洛模拟，即用大量随机试验得出最终结果，我们将会在后续予以介绍。

这道题可以采用0-1整数规划的思维，对应的方程组·为：
{ max ⁡ y = ∑ i = 1 10 c i x i ∑ i = 1 10 w i x i ≤ 30 x i ∈ { 0 , 1 } ， i = 1 , 2 , ⋯ 10 \left \{

$$\begin{array}{c} \max y=\sum\limits_{i=1}^{10}c_i x_i \\ \sum\limits_{i=1}^{10} w_i x_i \leq 30 \\ x_i \in \{ 0,1 \}，i=1,2, \cdots 10 \end{array}$$ \right. ⎩ ⎨ ⎧​maxy=i=1∑10​ci​xi​i=1∑10​wi​xi​≤30xi​∈{0,1}，i=1,2,⋯10​
该题也是动态规划的经典例题，可以思考如何设计状态转移方程。

七、非线性规划

若存在某些限制 或 目标函数 不满足线性形式，即并非所有因素的幂次均为1，则称这类问题为非线性规划。在MATLAB中，同样有一个求解函数—— fmincon。

[x,val]=fmincon(@f, x0, A, b, Aeq, beq, lb, ub, @nonlcun, option)
1

对于新引入的变量：

• @f:目标函数，定义在外部
• x0:初始值，由于求出来的只能是局部最优解，所以必须设定
• option：求解方法，有interior-point（内点法）、sqp（序列二次规范法）、active-set（有效集法）、trust-region-reflective（信赖域反射算法）
• @nonlcun调用了一个定义在外部的非线性部分的约束（也可以定义在脚本内部）

虽然该算法只能求出局部最优解，但是如果设立多个初始值，则可能得到全局最优解（类似于蒙特卡洛的思想）。具体实现参看下图（对不起我太懒了）

八、图论与最短路算法

关于这部分内容，由于信息竞赛出身，就不赘述了。这里列出一条大纲，可以根据大纲配合网络资料学习 （以后有时间再回来填坑）

1. 图的定义与表示方法
2. 最短路模型与求解算法：
   （1） Dijkstra算法：非负边权单源最短路，贪心思想
   （2） Bellman_Ford算法：能解决存在负边权的模型并判断是否有负权环（此时无最短路），利用n-1次松弛求出最优解，复杂度较大
   （3）Floyd算法：三重循环给出无向图任意两点间的最短距离

在MATLAB中，已经有了对应的求解函数——shortestpath

[P, d] = shortestpath (G, start, end, ['method', algorithm])
1

• G：要求解的图
• start：起始节点
• end：目标节点
• [‘method’，algorithm]
• P：最短路径经过的节点序列
• d：最短距离

九、网络流模型

碎碎念：学信息竞赛的时候就被这种算法的魅力折服了，因为它真的好神奇，但也是真的难……
求解算法的实现就是简单的模板，难就难在如何判断某个问题能否划归为网络流问题，
关于这部分的算法原理以及C++实现，参看我之前写的一篇文章，包含最大流、最小割与费用流。
在MATLAB中，求解最大流的函数为 graphmax

[Maxflow, FlowMatrix, Cut] = graphmaxflow(matrix, vs, vt)
1

• matrix： 要求解的有向图的稀疏矩阵
• vs；源点
• vt：汇点
• Maxflow：最大流
• FlowMatrix：包含每条边所有流量的稀疏矩阵
• Cut：计算发点与终点的最小切合后链接起点的逻辑行向量，如果有多个解表示为一个矩阵

对于费用流，不断地在残量网络上进行最短路算法即可。

十、微分方程

微分方程与其说是一个模块，不如说是一种思想，一类问题。面对变化趋势与函数值关联的问题时，我们往往需要建立微分方程来进行刻画，例如：传染病模型，弹簧振动方程等等。
(后面有机会再填坑吧）

十一、插值与拟合

插值与拟合都是分析数据大致分布的方式，数据少而精确，往往用插值；数据多而低精度，往往用拟合。
插值与拟合在MATLAB中都有相应的算法，需要根据具体需要选择合适的拟合方式。
（后面有机会再回来填坑）

十二、时间序列

时间序列，简而言之，就是按时间顺序排列的数据。这里探讨的就是如何根据时间序列对未来进行预测

12.1 分类及其结合

根据时间对数据的影响，可分为以下几类：

• 长期趋势变动：长期上升或下降
• 季节变动：仅受季节因素影响
• 循环变动：存在周期性
• 随机变动：与时间无关

我们常见的数据往往是这四类数据以 加法 或 乘法 的方式结合起来，例如：若“羽绒服销量”呈现同比增长的趋势，则应将长期趋势变动程度 与 季节变动程度 相乘。

12.2 时间序列预测方法

某一时刻的数据受最近的数据影响程度较大，而受相对久远的数据影响较小。根据这一思想，要想构建函数进行分析预测，需要引入与时间有关的权值函数 $f(t)$。
长期趋势变动类型中，$f(t)$需要对 t 单调递增。最简单的思路是取 待预测时间点的前n项 取平均，其余项不做考虑，该方法称为 移动平均法（MA），此时
f n ( t ) = { 1 n , t 0 − t ≤ n 0 , t 0 − t > n f_n(t) =

$$\begin{cases} \frac{1}{n},t_0-t \leq n\\ 0,t_0-t>n \end{cases}$$ fn​(t)={n1​,t0​−t≤n0,t0​−t>n​
同样地，可以对移动平均法得到的数据进行二次或更多次移动平均。
这种方式显然是十分粗暴的。一般来说，$f(t)$是随t增大而平滑上升的。于是，我们构建这样一个递推式:
$${\hat{x}}_{t+1}=\alpha x_t+\alpha (1-\alpha ){\hat{x}}_t$$

• ${\hat{x}}_{t+1}$ ：对第t+1个时刻的预测值
• $x_t$表示t时刻的观测值。此时，越接近当下的数据项，其对应的权重越大

这种方法称为简单指数平滑法，$\alpha$被称为平滑系数。这种方法无法分析具有趋势的情况，只是单纯地调整了权重。

12.3 趋势模型

在此基础上，我们引入霍特线性趋势模型，其具体方程组如下：
{ l t = α x t + ( 1 − α ) ( l t − 1 + b t − 1 ) （水平平滑方程） b t = β ( l t − l t − 1 ) + ( 1 − β ) b t − 1 （趋势平滑方程） x ^ t + h = l t + h b t （预测方程） \left\{

$$\begin{array}{} l_t=\alpha x_t + (1 - \alpha )(l_{t-1} + b_{t-1})（水平平滑方程） \\ b_t=\beta(l_t - l_{t-1}) + (1 - \beta)b_{t-1}（趋势平滑方程）\\ \hat x_{t+h} = l_t + hb_t （预测方程） \end{array}$$ \right . ⎩ ⎨ ⎧​lt​=αxt​+(1−α)(lt−1​+bt−1​)（水平平滑方程）bt​=β(lt​−lt−1​)+(1−β)bt−1​（趋势平滑方程）x^t+h​=lt​+hbt​（预测方程）​

• $t$：当前期
• $h$：预测超前期数
• $x_t$：第t期的实际观测值
• $l_t$：时刻t的预估水平
• $b_t$：时刻t的预测趋势
• $\alpha$：水平的平滑参数
• $\beta$：趋势的平滑参数

我们已经能够看到，这里在继承简单平滑法（水平平滑）的基础上，又引入了对未来趋势的分析（趋势平滑），平滑系数则起到了平衡“惯性”（旧有预测）与“变化”（新增信息）的作用。
然而，这个模型依旧是不够理想的，且往往会出现预测过度的情况，所以后人又增加了“阻尼系数”的概念，用以抑制未来某一时刻的趋势，这一模型称为阻尼趋势模型。
{ l t = α x t + ( 1 − α ) ( l t − 1 + ϕ b t − 1 ) （水平平滑方程） b t = β ( l t − l t − 1 ) + ( 1 − β ) ϕ b t − 1 （趋势平滑方程） x ^ t + h = l t + ( ϕ + ϕ 2 + ⋯ + ϕ h ) b t （预测方程） \left\{

$$\begin{array}{} l_t=\alpha x_t + (1 - \alpha )(l_{t-1} + \phi b_{t-1})（水平平滑方程） \\ b_t=\beta(l_t - l_{t-1}) + (1 - \beta)\phi b_{t-1}（趋势平滑方程）\\ \hat x_{t+h} = l_t + (\phi + \phi ^2 + \cdots +\phi ^h)b_t （预测方程） \end{array}$$ \right . ⎩ ⎨ ⎧​lt​=αxt​+(1−α)(lt−1​+ϕbt−1​)（水平平滑方程）bt​=β(lt​−lt−1​)+(1−β)ϕbt−1​（趋势平滑方程）x^t+h​=lt​+(ϕ+ϕ2+⋯+ϕh)bt​（预测方程）​

• $t$：当前期
• $h$：预测超前期数
• $x_t$：第t期的实际观测值
• $l_t$：时刻t的预估水平
• $b_t$：时刻t的预测趋势
• $\alpha$：水平的平滑参数
• $\beta$：趋势的平滑参数
• $\varphi$：阻尼系数，$0<\varphi \le 1$

容易发现，当阻尼系数为1时，该模型退化为霍特线性模型。

12.4 季节性（周期性）模型

季节性的预测不同于上个模型的环比推导，更注重同期的比较。
{ l t = α ( x t − s t − m ) + ( 1 − α ) l t − 1 , （水平平滑方程） s t = γ ( x t − l t − 1 ) + ( 1 − γ ) s t − m , （季节平滑方程） x ^ t + h = l t + s t + h − m ( k + 1 ) , k = ⌊ h − 1 m ⌋ , （预测方程）

$$\begin{cases} l_t=\alpha (x_t-s_{t-m})+(1 - \alpha )l_{t-1},（水平平滑方程）\\ s_t = \gamma (x_t - l_{t-1}) + (1 - \gamma )s_{t-m},（季节平滑方程）\\ \hat x_{t+h} = l_t + s_{t+h-m(k+1)},k=\lfloor \frac{h-1}{m} \rfloor ,（预测方程）\\ \end{cases}$$ ⎩ ⎨ ⎧​lt​=α(xt​−st−m​)+(1−α)lt−1​,（水平平滑方程）st​=γ(xt​−lt−1​)+(1−γ)st−m​,（季节平滑方程）x^t+h​=lt​+st+h−m(k+1)​,k=⌊mh−1​⌋,（预测方程）​

• $m$：周期长度
• $\alpha$：水平平滑参数
• $\gamma$：季节平滑参数
• $h$：预测超前期数
• ${\hat{x}}_{t+h}$：第h期的预测值

与上一个模型复合起来，得到的模型可以更好地刻画增长趋势。
霍特-温特季节性加法模型，应用于较为稳定的、有季节性特征的情况
{ l t = α ( x t − s t − m ) + ( 1 − α ) ( l t − 1 + b t − 1 ) , （水平平滑方程） b t = β ( l t − l t − 1 ) + ( 1 − β ) b t − 1 , （趋势平滑方程） s t = γ ( x t − l t − 1 − b t − 1 ) + ( 1 − γ ) s t − m , （季节平滑方程） x ^ t + h = ( l t + h b t ) s t + h − m ( k + 1 ) , k = ⌊ h − 1 m ⌋ , （预测方程）

$$\begin{cases} l_t=\alpha (x_t-s_{t-m})+(1 - \alpha )(l_{t-1} + b_{t-1}),（水平平滑方程）\\ b_t=\beta(l_t-l_{t-1})+(1 - \beta)b_{t-1},（趋势平滑方程）\\ s_t = \gamma (x_t - l_{t-1} - b_{t-1}) + (1 - \gamma )s_{t-m},（季节平滑方程）\\ \hat x_{t+h} = (l_t + hb_t) s_{t+h-m(k+1)},k=\lfloor \frac{h-1}{m} \rfloor ,（预测方程）\\ \end{cases}$$ ⎩ ⎨ ⎧​lt​=α(xt​−st−m​)+(1−α)(lt−1​+bt−1​),（水平平滑方程）bt​=β(lt​−lt−1​)+(1−β)bt−1​,（趋势平滑方程）st​=γ(xt​−lt−1​−bt−1​)+(1−γ)st−m​,（季节平滑方程）x^t+h​=(lt​+hbt​)st+h−m(k+1)​,k=⌊mh−1​⌋,（预测方程）​
霍特-温特季节性阻尼趋势模型，应用于稳定性不太强的情况。
$$\{\begin{array}{l}{l}_t=\alpha \frac{x_t}{s_{t-m}}+(1-\alpha )(l_{t-1}+\varphi b_{t-1}),（水平平滑方程）\\ b_t=\beta (l_t-l_{t-1})+(1-\beta )\varphi b_{t-1},（趋势平滑方程）\\ s_t=\gamma \frac{x_t}{l_{t-1}+\varphi b_{t-1}}+(1-\gamma )s_{t-m},（季节平滑方程）\\ {\hat{x}}_{t+h}=(l_t+(\varphi +{\varphi }^2+\cdots +{\varphi }^h)b_t)s_{t+h-m(k+1)},k=\lfloor \frac{h-1}{m}\rfloor ,（预测方程）\end{array}$$

12.5 时间序列的平稳性

平稳性指性质与观测时间无关的序列，因此，具有趋势或季节性的时间序列不是平稳的。值得注意的是，仅具有周期性行为的时间序列是平稳的，因为周期长度不固定，我们预先不知道周期的波峰或波谷在哪。
若时间序列 { x t } \{x_t \} {xt​}满足以下三个条件：
(1)$E(x_t)=E(x_{t-s})=u$（均值为固定常数)
(2)$Var(x_t)=Var(x_{t-s}={\sigma }^2$（方差存在且为常数）
(3)$Cov(x_t,x_{t-s})={\gamma }_s$（协方差只和间隔s有关，与t无关）
则称$x_t$为 协方差稳定 。
要判断序列是否平稳，我们一般通过检查自相关系数ACF，定义:$${\rho }_s=\frac{cov(x_1,x_{t-s})}{\sqrt{Var(x_t)}\sqrt{Var(x_{t-s})}}=\frac{{\gamma }_s}{{\gamma }_0},其中{\gamma }_0=Cov(x_1,x_t)=Var(x_t)$$
这个量可以用来衡量间隔s的一对变量的相关系数。然而，由于某一个变量可能受其前面多个元素影响，对于$x_t={\varphi }_{s1}{x}_{t-1}+\cdots {\varphi }_{ss}{x}_{t-s}+e_t$,其中\phi _{ss}就是我们关注的偏自相关系数PACF。

十三、旅行商问题

有一个旅行商需要拜访每一个小镇后再回到家里，小镇之间两两可达且需要一定路费，该如何规划路线才能在总花费最少呢？
将这个问题转化为图论问题，即在一张带权无向图中，遍历所有点后返回的最小代价。

一个简单的贪心思路是：先随便找一个圈，再思考如何对这个圈进行改良。对于初始圈，只需要给出这个点的排列就可以了，用$w(i,j)$指代边权。
具体地，对于排列$P$，若存在 $i,j$ 满足: $$w(i,i+1)+w(j,j+1)<w(i,j)+w(i+1,j+1)$$ 则断开 $(i,i+1)$ 与 $(j,j+1)$，连接 $(i,j)$ 与 $(i+1,j)$，重新构筑起一个圈。体现在排列$P$上，即将$p_{i+1},p_{i+2},\cdots p_j$ 翻转（或称颠倒顺序），形成新的排列$P^{\prime }$。这种算法称为改良圈算法。在实际应用中，为了增大正确性，往往需要设立多个初始圈执行算法。

当然，例子中的旅行商若增大到多个，那么模型又会有新的变化，这个我们日后再讨论。

十四、聚类分析入门

聚类分析，即将数据进行分类，“离得近”的认为是一组，“离得远”的集合彼此分开。这里定义“远近”所需的“距离”往往需要根据我们自己的定义。

14.1 相关性的刻画——距离及其他

将数据视作一个个向量，再对向量定义各种各样的“距离”。下面列举一些常见距离：

• 绝对值距离（曼哈顿距离）：$\sum _{k=1}^n\left|x_k-y_k\right|$
• 欧氏距离：$\sqrt{\sum _{k=1}^n(x_k-y_k{)}^2}$
• Minkowski距离：$[\sum _{k=1}^n(x_k-y_k{)}^p{]}^{\frac{1}{p}}$
• 切比雪夫距离：$\underset{1\le k\le p}{\max }\left|x_k-y_k\right|$

这里再引入两个概念。第一个为 (线性)相关系数 $$r=\rho (X_i,X_j)=\frac{{\sum }_{k=1}^n(x_{ki}-{\bar{x}}_i)(x_{kj}-{\bar{x}}_j)}{\sqrt{{\sum }_{k=1}^n(x_{ki}-{\bar{x}}_i{)}^2}\sqrt{{\sum }_{k=1}^n(x_{kj}-{\bar{x}}_j{)}^2}}$$
可以看到，分子上为协方差 $cov(X_i,X_j)$，也即两个平移后向量的内积；分母上为两个向量各自标准差相乘 $\sigma (X_i)\cdot \sigma (X_j)$，也即两个平移后的向量对应模长相乘。由此也容易得知，$r\in [-1,1]$
r>0，对应正相关；r<0，对应负相关；r=0，认为不具备线性相关性。
注意：相关系数仅用于判断线性相关程度，不能说明两个向量不具备其他类型的相关性

第二个概念为夹角余弦，即未经平移后的向量直接求夹角余弦：
$$cos⟨X_i,X_j⟩=\frac{{\sum }_{k=1}^n(x_{ki}-{\bar{x}}_i)(x_{kj}-{\bar{x}}_j)}{\sqrt{{\sum }_{k=1}^n{x}_{ki}^2}\sqrt{{\sum }_{k=1}^n{x}_{kj}^2}}$$

在进行聚类后，我们还需要分析两个类之间的相关性，所以还需要对两个类（设为 $G_p$,$G_q$）之间的距离$D(G_p,G_q)$进行定义。

• 重心距离：定义类的重心$\bar{X}=\frac{1}{n}{\sum }_{k=1}^n{X}_k$，分别取两个类的重心${\bar{X}}_p$,${\bar{X}}_q$，将这两个向量的距离作为两个类之间的距离
  $$D(G_p,G_q)=d({\bar{X}}_p,{\bar{X}}_q)$$
• 平均距离：枚举两个类中的向量，将所有可能的距离算出后取平均
  $$D(G_p,G_q)=\frac{1}{||G_p||\cdot ||G_q||}\sum _{X\in G_p,Y\in G_q}d(X,Y)$$
• 最短距离法：取两个类间最短的距离作为类间距离
• D ( G p , G q ) = min ⁡ X ∈ G p , Y ∈ G q { d ( X , Y ) } D(G_p,G_q) = \min_{X \in G_p,Y \in G_q} \{ d(X,Y) \} D(Gp​,Gq​)=X∈Gp​,Y∈Gq​min​{d(X,Y)}

在实际算法中，我们可以认为孤点看作类。

14.2 聚类分析算法

最短距离系统聚类法

该算法遵循一个简单的思路：每次取一对距离最近的类，将这两个类合并，重复上述过程直至最后仅剩一类。
我们可以通过设置迭代次数改变最终类簇数目，设置迭代次数为$k$，则最终剩余 $n-k$ 个类簇。

K_means算法

该算法流程为：

1. 给定k个初始类簇中心点，将所有数据点与最近的类簇中心点归为一类；
2. 计算k个类簇的重心及其与所处类簇的中心点的距离，分两种情况：
   (1) 若这些距离足够小 或 迭代次数达到限制，则停止迭代；
   (2) 若这些距离较大，则将这k个重心作为新的类簇中心点，重复上述过程，直至满足第一种情况。

注意到该算法初值的设定直接影响最终结果，而这种预设又与数据集本身无关，于是对于初值设定做如下改进：

1. 随机选取某一样本作为第一个聚类中心
2. 计算每个样本与当前已有聚类中心的最短距离，该距离越大说明其作为新聚类中心的概率越大。计算结束后，用轮盘发（依据概率大小进行抽选）选出下一个聚类中心
3. 重复步骤二，直至选出k个聚类中心

DBSCAN算法

该算法根据密度进行聚类，定义核心点为半径E范围内数据点个数不小于规定值M的点。具体流程为：

1. 筛选出数据集中的核心点，并将核心点与半径E内的点合并为一类
2. 若核心点附近存在其他核心点，则这两个点合并成同一集合
3. 重复上述流程，直至所有核心点均被访问

在这个过程中，E与M的取值直接决定了聚类效果

十五、BP神经网络

神经网络，哇，一听就非常的高大上（bushi）
实际上，算法领域的神经网络算法就是模拟生物学上的神经网络构建出的模型，只要了解了其中本质就会发现，这个算法的难度比想象中要小很多。
下面各个量的表示方法可能不够标准，大家多包容QwQ。

15.1 组成及结构

生物学上，神经网络最基本的单位就是神经元，单个神经元需要具备以下功能：接受信息后产生电信号，在电信号强度超过阈值时对外传播信息。从宏观来看，一条反射弧包括感受器，效应器以及中间的信息处理部分。这种描述是不够严谨的，却反映出神经网络算法的大致结构。

整体上看，神经网络分为输入层（对应感受器）、隐藏层（对应信息处理部分）、输出层（对应效应器）。不难发现，这种构造与函数不谋而合，这三层依次对应 原像、对应关系 与 像。

对于隐藏层的每一个神经元，其所接受的“电信号”为$\sum x_i{w}_{ij}$，$x_i$表示某个神经元传递的信号，$w_{ij}$表示从信号来源 i 到本神经元 j 之间路径上的权值。
除去接收的信息，阈值（或称偏置，用 $\theta$ 表示）与处理传出信号的方式仅与每个神经元本身有关。这里引入激活函数的概念，$y_i=f(\sum x_i{w}_{ij}-\theta )$即为输出信号。常用的激活函数为两种，第一种为S型函数： $$f(x)=\frac{1}{1+e^{-\alpha x}}$$
取值范围为$(0,1)$；
第二种为双极S型函数：
$$f(x)=\frac{2}{1+e^{-\alpha x}}-1$$
取值范围为$(-1,1)$。
这里附图：

15.2 学习方式——反馈机制

（P.S. 不一定对，有待商榷，参考着看）

对原有网络进行优化的过程即为学习，相应地，将已经成型的网络用于处理数据即为工作。对于学习过程，按照是否提供训练集分为 监督学习 与 无监督学习，BP神经网络属于前者，通过预测值与实际值对照后的反馈来调节。
具体地，定义误差函数：$$e=\frac{1}{2}\sum _{i=1}^q(y_i-z_i{)}^2=\frac{1}{2}(Y-Z{)}^T(Y-Z),Y=[y_1,y_2,\cdots y_q{]}^T$$
这里用 $Z$ 表示神经网络求出的结果，$Y$表示实际结果，我们的目标就是让误差函数尽可能地小。由于我们希望通过学习可以调整 $w_i$ ，所以我们建立误差函数与权值的直接关系：
e = 1 2 ∑ i = 1 q ( y i − z i ) 2 = 1 2 ∑ i = 1 q ( f ( ∑ j = 1 p w j i x j − θ i ) − z i ) 2

$$\begin{align*} e &= \frac{1}{2}\sum_{i=1}^q (y_i - z_i)^2 \\ &=\frac{1}{2}\sum_{i=1}^q ( f(\sum_{j=1}^p w_{ji} x_j - \theta_i) - z_i)^2 \end{align*}$$ e​=21​i=1∑q​(yi​−zi​)2=21​i=1∑q​(f(j=1∑p​wji​xj​−θi​)−zi​)2​
对误差函数求偏导则有：
$$\begin{array}{rl}\frac{\partial e}{\partial w_{ji}}& =\frac{\partial e}{\partial y_i}\frac{\partial y_i}{\partial w_{ji}}\\ & =\frac{1}{2}\cdot 2(y_i-z_i)\cdot f^{{}^{\prime }}(\sum _{j=1}^p{w}_{ji}{x}_j-{\theta }_i)\cdot x_j\end{array}$$
若$f(x)=\frac{A}{1+e^{-\frac{x}{B}}}$，则$f^{{}^{\prime }}(x)=\frac{f(x)[A-f(x)]}{AB}$，代入公式，令$S_i={\sum }_{j=1}^p{w}_{ji}{x}_j-{\theta }_i$则：
$$\frac{\partial e}{\partial w_{ji}}=(y_i-z_i)\frac{f(S_i)[A-f(S_i)]}{AB}\cdot x_j={\delta }_i{x}_j$$

特殊地，对于输入层，由于没有输入信号，所以不存在激活函数及阈值；对于输出层，如果希望进行聚类分析，则采用（双极）S型函数，如果是求数值解，则令$f({\sum }_{j=1}^p{w}_{ji}{x}_j-{\theta }_i)={\sum }_{j=1}^p{w}_{ji}{x}_j$,也即 ${\theta }_i=0,f(x)=x$。

我们的目标是让误差函数$e$尽可能小，一个类似于牛顿迭代法的思路是，不断沿着当前点的切线方向移动，直至移动到某一极值点，这种方法也即梯度下降法。具体地，对于某一个权值$w_{ji}$，有：
$$w_{ji}=w_{ji}-\eta \frac{\partial e}{\partial w_{ji}}+\alpha \Delta w_{ji}$$

• $\eta \frac{\partial e}{\partial w_{ji}}$：调整量，其中$\eta$为学习率，实际计算时$\frac{\partial e}{\partial w_{ji}}={\delta }_i{x}_j=(y_i-z_i)\frac{f(S_i)[A-f(S_i)]}{AB}\cdot x_j$
• $\alpha \Delta w_{ji}$：平滑项，$\alpha$ 为平滑系数， $\Delta w_{ji}$ 表示上一次 $w_{ji}$的变化量。

这样每个权值就完成里其对应的修改。
讨论完权值后，对阈值$\theta$的调整相对简单
∂ e ∂ θ i = ∂ e ∂ y i ∂ y i ∂ θ i = ( y i − z i ) ⋅ ( − f ′ ( S i ) ) = − δ i

$$\begin{aligned} \frac{\partial e}{\partial \theta_i} &= \frac{\partial e}{\partial y_i} \frac{\partial y_i}{\partial \theta_i} \\ &=(y_i-z_i)\cdot (-f^{'}(S_i))\\ &= - \delta _i \end{aligned}$$ ∂θi​∂e​​=∂yi​∂e​∂θi​∂yi​​=(yi​−zi​)⋅(−f′(Si​))=−δi​​
相应地，有：
$${\theta }_i={\theta }_i-\eta \frac{\partial e}{\partial {\theta }_i}+\alpha \Delta {\theta }_i$$
各参数意义同上。

15.3 流程

学习流程中，对于给定的一组数据：

1. 初始化权值与阈值，阈值范围为$[-1,1]$内随机数，权重取0附近的随机数即可；

2. 从输入层开始，逐层计算各个节点的输入与输出；

3. 根据实际值与预测值计算误差函数值$e$，若$e$足够小则结束迭代；

4. 若误差较大，则从输出层开始反向更新权值与阈值。对于第k层与第k+1层之间的权值，$\frac{\partial e}{\partial w_{ij}}=\frac{\partial e}{\partial x_j(k+1)}\frac{\partial x_j(k+1)}{\partial w_{ij}}$，其中$\frac{\partial x_j(k+1)}{\partial w_{ij}}={\delta }_j(k+1)$， 而对于$\frac{\partial e}{\partial x_j(k+1)}=\frac{\partial e}{\partial }$，考虑矩阵偏导
   $$\frac{\partial e}{\partial X_{k+1}}=\frac{\partial e}{\partial X_{k+2}}\frac{\partial f(W\cdot X_{k+2}-{\Theta }_{k+1})}{\partial X_{k+1}}$$