夏至的个人博客

数学建模笔记(司守奎)

发表于 更新于
本文字数： 25k 阅读时长 ≈ 23 分钟

数学建模笔记（司守奎）

线性规划

  1. 定义： 在一组线性条件限制下，求一线性目标函数最大或最小的问题
  2. 标准的matlab形式：\(min f^T x,\) \(s.t. \begin{cases} A\cdot x \le b\\ Aeq \cdot x = beq \\ lb \le x \le ub \end{cases}\) 其中\(f,x,b,beq,lb,u\)为列向量，\(A，Aeq\)为矩阵
  3. matlab相关求解命令为 \([x,val]=linprog(f,A,b,Aeq,beq,lb,ub)\)

整数规划

  1. 定义： 数学规划中的变量（部分或全部）限制为整数时，称为整数规划
  2. 分类： 根据变量整数是部分还是全部可分为纯整数规划和混合整数规划
  3. 求解方法：
    1. 分支定界法（纯或混）
    2. 割平面法（纯或混）
    3. 隐枚举法（01规划）
      • 过滤隐枚举法
      • 分支隐枚举法
    4. 匈牙利法（指派问题）
    5. 蒙特卡洛法（各种类型）

0-1型整数规划

  1. 定义：整数规划的一种特殊情况，变量仅取值0或1 应用范围可分为以下几种情况

相互排斥的约束条件

也就是说题目中的某一种条件只有一个量，如果给了这个1，同类的都为0，较为典型的问题有运输问题，只用一种方式运输，用火车运了，其他的运输栏都为0，其约束条件可进一步简化为 \[y_i=\begin{cases} 1,第i个元素起约束作用\\ 0,第i个元素不起作用，i=1,2.....,m\\ \end{cases}\] \[a_{i1}x_1+...+a_{in}x_n\le b_i+(1-y_i)M,i=1,2...,m,\\ y_1+\dots+y_m=1 \] 由约束条件很容易看出，当\(y_i\)等于1,就只有这个约束起作用，其他的都是多余的

固定费用的问题

在讨论线性规划时，有些问题要求固定费用，这种问题可以通过改变为混合整数规划来解决，数学模型可表示为 \(y_i\epsilon\le x_i\le y_i M\) 其中\(\epsilon\)为充分小的正常数；M为充分大的正常数，表明\(x_i>0\)时，\(y_i\)必须为1，\(x_i=0\)时\(y_i\)必须为0,\([x_i]\)表示采用i方式生产时产量，\(y_i\)表示是否用第i种方式生产

指派问题

指派问题描述的是分配n个人去做n件事情，每个人做且仅做一件事情，且分配第i个人去做第j件事情，花费\(C_{ij}\)单位时间，求如何分配使总时间最小，这类问题的关键就是要求出分配矩阵，数学形式可表现为 \[x_{ij}=\begin{cases} 1,第i人做第j项工作\\ 0,第i人做第j项工作\\ \end{cases}\] 数学模型为：\[min \sum_{i=1}^N \sum_{j=1}^N c_{ij}x_{ij}\\ s.t. \begin{cases} \sum_{i=1}^N x_{ij}=1,i=1,2,...,n\\ \sum_{j=1}^N x_{ij}=1,j=1,2,...,n\\ x_{ij}=0 or 1,i,j=1,...,n \end{cases}\]

蒙特卡洛法（随机取样法）

  1. 又被称为计算机随机模拟法，它是基于对大量数据的统计结果来实现一些确定性问题的计算
  2. 使用该方法必须使用计算机生成相关分布的随机数

整数线性规划的计算机求解

  1. 整数规划的求解用Lingo等专用软件比较方便，对于整数线性规划也可以用matlab的intlinprog函数求解，但其的缺点是必须把所有的决策变量化为一维决策变量，变量替换后，约束条件很难写出，最好用lingo
  2. matlab求解混合整数线性规划的命令是
1
[x,fval]=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)

对应以下数学模型 \[min_x f^Tx,\\ s.t. \begin{cases} x(intcon)为整数\\ A \cdot x \le b,\\ Aeq \cdot x =beq,\\ lb \le x \le ub \\ \end{cases}\] 式中：\(f,x,intcon,b,beq,lb,ub为列向量；A，Aeq为矩阵\)

非线性规划

非线性规划模型

  1. 定义：如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题

  2. 通过投资决策问题归纳非线性规划数学模型的一般形式 总资金A元，投资第i个项目花\(a_i\)元，预计可收益\(b_i\)元 选择最佳投资方案 投资决策变量 \[x_i=\begin{cases} 1,决定投资第i个项目\\ 0,决定不投资第i个项目\\ \end{cases}\] 则该模型可用下列数学模型表示 \[max\, Q = \frac{\sum_{i=1}^n b_i x_i} {\sum_{i=1}^n a_i x_i}\\ s.t. \begin{cases} 0 < \sum_{i=1}^n a_i x_i \le A ,\\ x_i(1-x_i)=0,i=1,...,n\\ \end{cases}\]

  3. 根据2中例题，非线性规划问题可进一步概括为： \(min\,f(x)\\ s.t. \begin{cases} h_j(x)\le0,j=1,2,...,q\\ g_i(x)=0,i=1,2,...,p\\ \end{cases}\) 其中\(x=[x_1,...,x_n]^T\)为模型的决策变量，\(f\)为目标函数，\(g_i和h_j\)为约束函数，\(g_i(x)=0\)为等式约束，\(h_j(x)\le0\)为不等式约束

  4. 对一个实际问题，要将其规为非线性规划问题时，一般要注意以下几点

    1. 确定供选方案
    2. 提出追求目标
    3. 给出价值标准
    4. 寻求限制条件

  5. 线性规划与非线性规划的区别：线性规划最优解只能在可行域的边界上达到（特别是顶点），而非线性规划最优解可在可行域任一点达到

  6. 非线性规划的matlab表示 \[minf(x)\\ s.t. \begin{cases} A\cdot x \le b,\\ Aeq \cdot x=beq,\\ c(x)\le0\\ ceq(x)=0,\\ lb\le x \le ub \end{cases}\] 式中的\(f(x)\)为标量函数，\(A,b,Aeq,beq,lb,ub\)为相应维数的矩阵和向量，\(c(x),ceq(x)\)为非线性向量函数 matlab命令为

    1
    2
    3
    [x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
    # x返回决策变量x的取值，fval返回目标函数取值，fun是M文件自定义函数f(x),x0是x的初始值
    nonlcon是用M文件定义的c(x)ceq(x)，options定义优化参数

无约束问题的Matlab解法

  1. 在matlab工具箱中，用于求无约束极小值的函数有fminunc和fminsearch，用法分别为
1
2
[x,fval]=fminunc(fun,x0,options)
[x,fval]=fminsearch(fun,x0,options) #只能求初始值附近的一个极小值点

约束极值问题

  1. 定义：带有约束条件的极值问题，也叫规划问题

二次规划

  1. 定义：若某非线性规划的目标函数为自变量x的二次函数，约束条件又全为线性的，称这种规划为二次规划
  2. Matlab中二次规划的数学模型可表述为 \[min\, \frac{1}{2}x^THx+f^Tx,\\ s.t. \begin{cases} Ax\le b\\ Aeq \cdot x=beq,\\ lb\le x\le ub \end{cases} H为实对称矩阵，\\ \] matlab求解二次规划的命令为 
    1
    [x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)

罚函数法

  1. 利用罚函数法可以将非线性规划问题的求解转化为求一系列无约束极值的问题，也把这种方法叫做序列无约束最小化技术
  2. 罚函数求解非线性规划问题的思想是利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转换为无约束非线性规划问题，主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法
  3. 外罚函数法： 考虑问题：\[minf(x)\\ s.t.\begin{cases} g_i(x) \le 0,i=1,...,r\\ h_j(x) \ge 0,j=1,...,s\\ k_m(x) =0,m=1,...,t \end{cases}\]取一个充分大的数M>0,构造函数\(P(x,M)=f(x)+M\sum_{i=1}^nmax(g_i(x),0)-M\sum_{j=1}^nmin(h_j(x),0)+M\sum_{m=1}^n |k_m(x)|\),则以增广目标函数\(P(x,M)\)为目标函数的无约束极值问题\(minP(x,M)\)的最优解也是原问题的最优解
    1. 如果非线性规划问题要求实时算法，可以使用罚函数算法，但计算精度较低
    2. 如果不要求实时算法，要求高精度,可以使用lingo或matlab的fmincon命令求解

matlab求约束极值问题

  1. 在matlab工具箱中，用于求解约束最优化问题的函数有fminbnd,fmincon,quadprog,fseminf,fminimax函数

  2. fminbnd函数 

    1
    2
    3
    [x,fval]=fminbnd(fun,x1,x2,options) 
    # 用于求单变量非线性函数在[x1,x2]上极小值
    # 返回极小点x和函数的极小值

  3. fseminf函数 用于求下列模型 \[minf(x),\\ s.t. \begin{cases} A \cdot x \le b,\\ Aeq \cdot x =beq,\\ lb\le x \le ub\\ c(x)\le0\\ ceq(x)\le0\\ K_i(x,w_i)\le0,1\le i\le n \end{cases}\\其中c(x),ceq(x)为向量函数，K_i(x,w_i)为标量函数,w_1...为附加变量\] 

    1
    [x,fval]=fseminf(fun,x0,ntheta,seminfcon,A,b,Aeq,beq,lb,ub)

    fun定义目标函数f(x),x0为x初始值，ntheta是半无穷约束\(K_i(x,w_i)\)个数，函数seminfcon用于定义非线性不等式约束\(c(x)\),非线性等式约束\(ceq(x)\)和半无穷约束\(K_i(x,w_i)\)的函数，seminfcon有两个输入参量x,s,s是推荐的采样步长 可以不使用

  4. fminimax函数 用于求下列模型 \[min_xmax_iF_i(x),\\ s.t. \begin{cases} A \cdot x \le b,\\ Aeq \cdot x =beq,\\ lb\le x \le ub\\ c(x)\le0\\ ceq(x)=0\\ \end{cases}\] matlab命令为 

    1
    [x,fval]=fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)

图与网络模型及方法

图的基本概念与数据结构

  • 基本概念 图可由顶点集V和边集E表示，G=(V,E)，若各边加上方向则为有向图，否则为无向图，另外赋权图是指每条边都有一个(或多个)实数对应的图，该实数称为这条边的权，赋权图在实际情况中非常有用
  • 数据结构 图的数据结构主要有两种方法
    1. 邻接矩阵表示法
    2. 稀疏矩阵表示法

最短路问题

两个指定顶点之间的最短路径

问题简述：给定一个连接若干个城镇的铁路网络，在这个网络的两个指定城镇间，找一条最短铁路线 模型抽象：该问题我们可以构建一个赋权图G=(V,E,W),其中顶点集V表示各个小城镇，E为边的集合，邻接矩阵W表示各顶点间的距离，若顶点间无边，则w=∞，则问题可转换为求赋权图G中指定两个顶点\(u_0,v_0\)间具有最小权的路，这条路称为\(u_0,v_0\)间的最短路，它的权称为二者间的距离，记为\(d(u_0,v_0)\) 常用方法：

  1. 狄克斯特拉(Dijkstra)算法:基本思想为按距\(u_0\)从近到远为顺序，依次求得\(u_0\)到G的各顶点的最短路和距离，直至\(v_0\)（或直至G所有顶点）

两个指定顶点间最短路的数学规划模型

假设有向图有n个顶点，现要求从顶点\(v_1\rightarrow v_n\)的最短路，仍然用E表示弧的集合，设\(W=(w_{ij})_{n×n}\)为邻接矩阵，其分量为：\[w_{ij}=\begin{cases} 弧v_iv_j的权值，v_iv_j\in E\\ ∞，其他， \end{cases}\] 决策变量为\(x_{ij}\),当\(x_{ij}=1\),说明弧\(v_iv_j\)位于顶点\(v_1\)到顶点\(v_n\)的最短路上，否则\(x_{ij}=0\),其数学规划表达式为\[min\sum_{v_iv_j\in E}w_{ij}x_{ij},\\ s.t. \begin{cases} \sum_{j=1,v_iv_j\in E}^n x_{ij}-\sum_{j=1,v_jv_i\in E}^n x_{ji}=\begin{cases} 1,i=1\\ -1,i=n\\ 0,i\ne1,n \end{cases}\\ x_{ij}=0或1 \end{cases} \]

每对顶点之间的最短路径

计算赋权图各队顶点之间的最短路径，调用方法有两种

  1. Dijkstra算法 每次以不同的顶点作为起点，用Dijkstra算法算出该起点到其余顶点的最短路径，反复执行n-1次操作，就可得到从每一个顶点到其他顶点的最短路径，时间复杂度为\(O(n^3)\)
  2. Floyd算法 递推产生一个矩阵序列\(A_1,...,A_k,...,A_n\),其中矩阵\(A_k\)的第i行第j列元素\(A_k(i,j)\)表示从顶点\(v_i\)到顶点\(v_j\)的路径经过顶点序号不大于k的最短路径长度，计算时用迭代公式\(A_k(i,j)=min(A_{k-1}(i,j),A_{k-1}(i,k)+A_{k-1}(k,j))\),k是迭代次数，i,j,k=1,2,...n,最后当k=n时，\(A_n\)即是各顶点间的最短通路值

最小生成树问题

基本概念

连通的无圈图叫做树，记为T，其度为1的顶点称为叶子节点，显然有边的树至少有两个叶子顶点 若图G=(V(G),E(G))和树T=(V(T),E(T))满足\(V(G)=V(T),E(T)\subset E(G)\)，则称T是G的生成树，图G连通的充要条件是G有生成树

最小生成树

欲修筑连接n个城市的铁路，已知i城与j城之间的铁路造价为\(c_{ij}\),设计一个线路图使总造价最低，这个问题的数学模型是在连通赋权图上求权最小的生成树，赋权图具有最小权的生成树叫做最小生成树 常用算法： 构造连通赋权图G=(V,E,W)的最小生成树，设置两个集合P和Q,其中P用于存放G的最小生成树的顶点，集合Q存放G的最小生成树的边。令集合P的初值为\(P={V_1}\)(假设构造最小生成树时，从顶点\(v_1\))出发，集合Q的初值为\(Q=\varnothing\)

  1. prim算法 从所有\(p\in P,v\in V-P\)的边中，选取具有最小权值的边pv，将顶点v加入集合P中,将边pv加入集合Q中，如此重复，直到P=V时，最小生成树构造完毕，这是集合Q包含了最小生成树的所有边。
  2. Kruskal算法 选\(e_1\in E(G)\)，使得\(e_1\)是权值最小的边, 若\(e_1,e_2,...,e_i,e_{i+1}\)已选好，则从E(G)-{\(e_1,e_2,...,e_i\)}中选取\(e_{i+1}\),使得\({e_1,e_2,...,e_i,e_{i+1}}\)中无圈，\(e_{i+1}\)是E(G)-\({e_1,e_2,...,e_i}\)中权值最小的边，直到选得\(e_{|V|-1}\)为止

网络最大流问题

基本概念

  1. 网络：给一个有向图D=(V,A),其中A为弧集，在V中指定一点，称为发点(记为\(v_s\)),另一点称为收点(记为\(v_t\))，其余点为中间点，对于每一条弧\((v_i,v_j)\in A\)，对应有一个\(c(v_i,v_j)\ge 0(或简写为c_{ij})\),称为弧的容量，通常把这样的有向图叫做一个网络，记为D=(V,A,C),其中C={\(c_{ij}\)}
  2. 流：定义在弧集合A的一个函数\(f={f_{ij}={f(v_i,v_j)}}\)，并称\(f_{ij}\)为弧\((v_i,v_j)\)上的流量
  3. 最大流问题的线性规划表示

\[max\quad v(f),\\ s.t.\begin{cases} \sum_{j:(v_i,v_j)\in A}f_{ij}-\sum_{j:(v_j,v_i\in A)}f_{ji}=\begin{cases} v(f),i=s\\ -v(f),i=t\\ 0,i\ne s,t, \end{cases}\\ 0\le f_{ij} \le c_{ij},\forall(v_i,v_j)\in A \end{cases}\]

最小费用最大流问题

最小费用最大流

给定网络D=(V,A,C),每一弧\((v_i,v_j)\in A\)上，除了已给容量\(c_{ij}\)外，还给了一个单位流量的费用\(b(v_i,v_j)\ge 0(b_{ij})\)，所谓最小费用最大流问题就是求一个发点\(v_s\)到收点\(v_t\)的最大流，使流的总输送费用\(\sum_{(v_i,v_j)\in A}b_{ij}f_{ij}\)取最小值，最小费用最大值问题可以归结为两个线性规划问题，首先用一个线性规划模型求出最大流量\(v(f_{max})\)，再用一个线性规划模型求出最大流对应最小费用 \[min \sum_{(v_i,v_j)\in A}b_{ij}f_{ij}\\ s.t. \begin{cases} 0\le f_{ij} \le c_{ij},\forall(v_i,v_j)\in A\\ \sum_{j:(v_i,v_j)\in A}f_{ij}-\sum_{j:(v_j,v_i\in A)}f_{ji}=d_i\\ d_i=\begin{cases} v(f_{max}),i=s\\ -v(f_{max}),i=t\\ 0,i\ne s,t,\\ \end{cases}\\ \end{cases} \] \(v_{f_{max}}\)上节求得最大流流量

Matlab的图论工具箱

旅行商(TSP)问题

  1. 问题描述：一推销员准备前往若干城市推销商品，然后回到驻地，如何为他设计一条最短的旅行路线(从驻地出发，经过每个城市恰好一次，最后返回驻地)，这个问题就被称为旅行商问题
  2. 数学抽象：用图论术语表示，就是在一个赋权完全图中，找出一个有最小权的Hamilton圈，这个圈就叫最优圈
  3. 数学模型：设城市的个数为n，\(d_{ij}\)是两个城市i与j之间的距离，\(x_{ij}=0 or 1\)(1表示走过城市i到城市j的路，0表示没有选择走这条路)，则有 \[min \sum_{i\ne j}d_{ij}x_{ij}\\ s.t. \begin{cases} \sum_{j=1}^n x_{ij}=1,i=1,2,\dots,n,(每个点只有一条边出去)\\ \sum_{i=1}^n x_{ij}=1,i=1,2,\dots,n,(每个点只有一条边进去)\\ \sum_{i,j\in s}x_{ij}\le|s|-1,2\le|s|\le n-1,s\subset{1,2,\dots,n},即s为{1,2,\dots,n}的真子集\\ x_{ij}\in {0,1},i,j=1,2,\dots,n,i\ne j \end{cases} \]

计划评审方法和关键路线法

这两个方法是网络分析的重要组成部分，已经合并为一种方法，国外称PERT/CPM，国内称为统筹方法

计划网络图

计划网络图中通常用圆圈表示事件，用箭线表示工作，用这种方法画出的网络图为计划网络图 虚作业用虚箭线表示，表示工时为0，不消耗任何资源的虚构作业，作用只是为了正确表述工作的前向后继关系 关键路线: 在计划网络图中，称从初始事件到最终事件的由各项工作连贯组成的一条路为路线，具有累计作业时间最长的路线称为关键路线

时间参数

  1. 事件时间参数
    1. 事件的最早时间
    2. 事件的最迟时间
  2. 工作时间参数
    1. 工作的最早可能开工时间与最早可能完工时间
    2. 最迟必须开工时间和最迟必须完工时间
  3. 时差
    1. 工作总时差
    2. 工作单时差

计划网络图的计算

  1. 建立计划网络图
  2. 写出相应的规划问题 \[min\sum_{i\in V}x_i\\ s.t\begin{cases} x_j\ge x_i+t_{ij},(i,j)\in A,i,j\in V\\ x_i\ge 0,i\in V\\ \end{cases}\] V为所有时间集合，A为所有作业集合，\(x_i\)为事件i的开始时间，\(t_{ij}\)为作业(i,j)的计划时间
  3. 问题求解
  4. 将关键路线看为最长路 如果将关键路线看成最长路,则可以按照求最短路的方法(将求极小改为求极大)求出关键路线 设\(x_{ij}\)为0-1变量，当作业(i,j)位于关键路线上取1，否则取0，数学规划问题写成 \[max\sum_{(i,j)\in A}t_{ij}x_{ij}\\ s.t.\begin{cases} \sum_{j:(i,j)\in A}x_{ij}-\sum_{j:(j,i)\in A}x_{ji}=\begin{cases} 1,i=1,\\ -1,i=n,\\ 0,i\ne1,n,\\ \end{cases}\\ x_{ij}=0 or 1,(i,j)\in A \end{cases}\]

插值与拟合

若仅已知函数\(f(x)\)在某区间\([a,b]\)上一系列点上的值\(y_i=f(x_i),i=0,1,\dots,n\)，当需要在这些节点间的x上的函数值时，有两种方法：插值法，拟合法，二者都是根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上完全不同，需要我们根据实际情况考虑

插值方法

在工程和数学应用中，若在平面上给定一组离散点列，要求一条曲线，把这些点按次序连接起来，称为插值，若已知n+1点\((x_i,y_i)(i=0,1,...,n)\),下面详细介绍各种插值函数

分段线性插值

简单来说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记住\(I_n(x_i)=y_i\),且\(I_n(x)\)在每个小区间\([x_i,x_{i+1}]\)是线性函数\((i=0,1,...,n-1)\) \(I_n(x)\)可以表示\(I_n(x)=\sum_{i=0}^ny_il_i(x)\),其中 \[l_i(x)=\begin{cases} \frac{x-x_{i-1}}{x_i-x_{i-1}},x\in[x_{i-1},x_i],i\ne 0,\\ \frac{x-x_{i+1}}{x_i-x_{i+1}},x\in[x_i,x_{i+1}],i\ne n,\\ 0,else\\ \end{cases}\] \(l_n(x)\)有良好的收敛性，即对于\(x\in[a,b]\),有\(\lim_{n \to \infty}I_n(x)=f(x)\),用\(I_n(x)\)计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n无关，但n越大，分段越多，插值误差越小

拉格朗日插值多项式

拉格朗日(Lagrange)插值的基函数为\[l_i(x)=\frac{(x-x_0)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)}{(x_i-x_0)...(x_i-x_{i-1})(x_i-x_{i+1})...(x_i-x_n)}\\ =\prod_{j=0,j\ne i}^n \frac{x-x_j}{x_i-x_j},i=0,1,...,n\\ l_i(x)是n次多项式，满足l_i(x_j)=\begin{cases} 0,j\ne i\\ 1,j=i \end{cases}\] 综上，拉格朗日插值函数\[L_n(x)=\sum_{i=0}^ny_il_i(x)=\sum_{i=0}^ny_i(\prod_{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j})\]

样条插值

对于一些要求插值曲线具有较高的光滑程度，有连续曲率的数学物理问题，这时需要用样条插值

  1. 样条函数的概念 数学上将具有一定光滑性的分段多项式称为样条函数，具体地说，给定区间[a,b]的一个分划：\(\delta:a=x_0<x_1<\dots<x_{n-1}<x_n=b\) 如果函数S(x)满足
    1. 在每个小区间\([x_i,x_{i+1}](i=0,1,...,n-1)\)上S(x)是m次多项式
    2. S(x)在[a,b]上具有m-1阶连续导数 则称S(x)为关于分划\(\delta\)的m次样条函数，其图形为m次样条曲线
  2. 三次样条插值 利用样条函数进行插值，称为样条插值，三次样条插值函数就是把上面S(x)中的m=3,但是仅有这些条件是不足与确定所有待定参数，还需要额外的边界条件 常用的三次样条函数的边界条件有三种类型：
    1. \(S'(a)=y_0',S'(b)=y_n'\),有这种边界条件建立的样条插值函数称为\(f(x)\)的完备三次样条插值函数，若\(f'(x)\)不知道，我们可以要求S'(x)与f'(x)在端点处近似相等，这时以\(x_0,x_1,x_2,x_3\)为节点做一个三次Newton插值多项式,有这种边界条件建立的三次样条称为f(x)的Lagrange三次样条插值函数
    2. \(S"(a)=y_0",S"(b)=y_n",when y_0"=y_n"=0\),称为自然边界条件
    3. \(S'(a+0)=S'(b-0),S"(a+0)=S"(b-0)\),称为周期条件

Matlab插值工具箱

  1. 一维插值函数
1
2
3
4
5
6
y=interp1(x0,y0,x,'method')
method 指定插值方法，默认为线性插值，其值可为
'nearest' 最近项插值
'linear' 线性插值
'spline' 立方线条插值
'cubic' 立方插值
  1. 三次样条插值
1
2
pp=csape(x0,y0,conds,valconds);y=fnval(pp,x);
csape的返回值是pp形式，要求插值点函数值必须调用fnval，conds指定插值边界条件，valconds指定边界导数值

  1. 二维插值 若节点是二维的，插值函数就是二元函数，即曲面

    1. 插值节点为网格节点
    1
    2
    z=interp2(x0,y0,z0,x,y,'method')
    x0,y0分别为m,n维向量,表示节点,z0为n×m矩阵，表示节点值,x,y为一维数组，表示插值点，z表示得到的插值
    1. 插值节点为散乱节点
    1
    2
    ZI=griddata(x,y,z,XI,YI)
    x,y,z指数据点横纵竖坐标,XI,YI给定的插值点横纵坐标,ZI网格(XI,YI)的函数值

曲线拟合的线性最小二乘法

线性最小二乘法

曲线拟合：已知一组(二维)数据，即平面上n个点\((x_i,y_i),i=1,2,...,n,x_i\)互不相同，寻求一个函数(曲线)\(y=f(x)\),使\(f(x)\)在某种准则下雨所有数据点最为接近，即曲线拟合的最好。

线性最小二乘法是解决曲线拟合最常用的方法，基本思路是，令\(f(x)=a_1r_1(x)+a_2r_2(x)+\dots+a_mr_m(x)\),式中:\(r_k(x)\)为事先选定好的一组线性无关的函数，\(a_k\)为待定系数\((k=1,2,...,m<n)\)

拟合准则是使\(y_i(i=1,2,...,n)\)与\(f(x_i)\)的距离\(\delta_i\)的平方和最小，称为最小二乘准则

最小二乘法的Matlab实现

  1. 解方程组方法
  2. 多项式拟合法

最小二乘优化

在无约束优化问题中，有些重要的特殊情形，比如目标函数有若干个函数的平方和构成，这类函数可以写成\(F(x)=\sum_{i=1}^mf_i^2(x),x\in R^n,x=[x_1,...,x_n]^T,m\ge n\) 把极小化这类函数的问题\(min F(x)=\sum_{i=1}^mf_i^2(x)\)称为最小二乘优化问题 matlab也提供一些函数优化此问题,有lsqlin,leqcurvefit,leqnonlin,lsqnonneg,也可以直接调用工具箱里命令cftool,该命令给出一维数据拟合的交互式环境

曲线拟合和函数逼近

曲线拟合是指已知一组离散数据\({(x_i,y_i),i=1,...,n}\)，选择一个较简单的函数\(f(x)\)(如多项式),在一定准则(如最小二乘准则)下，最接近这些数据 如果已知一个较为复杂的连续函数\(y(x),x\in[a,b]\),要求选择一个较简单的函数\(f(x)\),在一定准则下最接近\(y(x)\)，就是函数逼近 与曲线拟合的最小二乘准则对应，函数逼近常用的一种准则是最小平方逼近，即\(J=\int_{a}^{b} [f(x)-y(x)]^2\, {\rm d}x\)达到最小，与曲线拟合一样，选一组函数构造\(f(x)\),代入使J达到最小

微分方程建模(感觉不太重要)

对于生活中的实际问题，我们往往需要将其化为微分方程的定解问题来求解，大致可分为以下几步：

  1. 根据实际要求要求确定要研究的量(自变量，未知函数，必要的参数)并确定坐标系
  2. 找出这些量所满足的基本规律(物理的，几何的，化学的或生物的等)
  3. 运用这些规律列出方程和定解条件 列方程常见方法：
    • 按规律直接列方程
    • 微元分析法与任意区域取积分方法
    • 模拟近似法

Matlab求微分方程的符号解

用matlab求解常微分方程的符号解，首先定义符号向量，然后调用命令dsolve，dsolve调用格式为

1
2
[y1,...,yN]=dsolve(eqns,conds,Name,Value)
eqns符号微分方程(组)，conds，初值条件或边值条件，Name,Value可选的成对参数

数理统计(感觉不太重要)

数理统计研究对象主要是受随机因素影响的数据，面对一批数据进行分析建模，首先要掌握参数估计和假设检验这两个数理统计中最基本方法，给定数据符合一定分布要求后，才能建立回归分析和方差分析等数学模型

时间序列

简介

将预测对象按照时间顺序排列起来，构成一个所谓的时间序列，从所构成的这一组时间序列过去的变化规律，推断今后变化的可能性及变化趋势，变化规律，就是时间序列预测法

时间序列法其实也是一种回归模型，其基于原理是：一方面承认事物发展的延续性，运用过去时间序列的数据进行统计分析就能推测事物发展趋势，另一方面又充分考虑到偶然因素影响产生的随机性

其优点是简单易行，便于掌握，能充分运用原时间序列的各项数据，计算速度快，对模型参数有动态确定能力，精度较好，采用组合的时间序列或将时间序列和其他模型组合效果更好

其缺点是不能反映事物内在联系，不能分析两个因素间相关关系，只适用于短期预测

确定性时间序列分析方法

时间序列预测技术就是通过对预测目标自身时间序列的处理，来研究其变化趋势，一个时间序列往往是一下几类变化形式的叠加或耦合

  1. 长期趋势变动，他指时间序列朝着一定的方向持续上升或下降，或停留在某一水平上的倾向，反映了客观事物的主要变化趋势。
  2. 季节变动
  3. 循环变动，通常是周期为一年以上，由非季节因素引起的涨落起伏波形相似的波动
  4. 不规则变动，通常分为突然变动和随机变动

通常用\(T_t\)表示长期趋势项，\(S_t\)表示季节变动趋势项，\(C_t\)表示循环变动趋势项，\(R_t\)表示随机干扰项，常见确定性时间序列模型类型如下 1. 加法模型： \(y_t=T_t+S_t+C_t+R_t\) 2. 乘法模型： \(y_t=T_t\cdot S_t\cdot C_t\cdot R_t\) 3. 混合模型： \[y_t=T_t\cdot S_t+R_t\\ y_t=S_t+T_t\cdot R_t\] \(y_t\)为观测目标的观测记录，均值\(E(R_t)=0\),方差\(Var(R_t)=\sigma^2\) 如果在预测时间范围以内，无突然变动且随机变动的方差较小，且有理由认为过去和现在的演变趋势将继续发展到未来时，可用一些经验方法进行预测

移动平均法

设观测序列为\(y_1,...,y_T\),取移动平均的项数N<T,一次移动平均值计算公式为： \[M_t^{(1)}=\frac{1}{N}(y_t+y_{t-1}+\dots+y_{t-N+1})\\=\frac{1}{N}(y_{t-1}+\dots+y_{t_N})+\frac{1}{N}(y_t-y_{t-N})=M_{t-1}^{(1)}+\frac{1}{N}(y_t-y_{t-N})\]

二次移动平均值计算公式为 \[M_t^{(2)}=\frac{1}{N}(M_t^{(1)}+\dots+M_{t-N+1}^{(1)})=M_{t-1}^{(2)}+\frac{1}{N}(M_t^{(1)}-M_{t-N}^{(1)})\] 当预测目标的基本趋势是在某一水平上下波动时，可用一次移动平均方法建立预测模型,即 \[\hat{y_{t+1} }=M_t^{(1)}=\frac{1}{N}(y_t+\dots+y_{t-N+1}),t=N,N+1,\dots,T\] 其预测标准误差为\[S=\sqrt{\frac{\sum_{t=N+1}^T(\hat{y_t}-y_t)^2}{T-N}}\]

最近N期序列值的平均值作为未来各期的预测结果，一般N取值为:\(5\le N\le 200\),当历史序列的基本趋势变化不大且序列中随机变动成分较多时，N的取值应较大一些,否则N的取值应小一些，在有确定的季节变动周期资料中，移动平均的项数应取周期长度，选择最佳N值的一个有效方法是编辑若干模型的预测误差，预测误差最小值为好

当预测目标基本趋势和某一线性模型吻合时，常用二次移动平均法，但序列同时存在线性趋势和周期波动时，可用趋势移动平均法建立预测模型： \(\hat{y_{T+m}}=a_T+b_Tm,m=1,2,\dots\),式中 \(a_T=2M_T^{(1)}-M_T^{(2)};b_T=\frac{2}{N-1}(M_T^{(1)}-M_T^{(2)})\)

指数平滑法

一般来说历史数据对未来值的影响是随时间间隔增长而递减的，而上述移动平均法的权数都是两端项权数小，中间项权数大，不符合上诉规律，更好的方法是对各期观测值依时间顺序进行加权平均作为预测值，指数平滑法就满足这一要求，且具有简单递推形式

指数平滑法根据平滑次数的的不同，又分为一次指数平滑法，二次指数平滑法，和三次指数平滑法等

一次指数平滑法 1. 预测模型

设时间序列为\(y_1,y_2,...,y_t,...,\alpha\)为加权系数,\(0<\alpha<1\),一次指数平滑公式为: \(S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}=S_{t-1}^{(1)}+\alpha(y_t-S_{t-1}^{(1)})\) 以这种平滑值进行预测，就是一次指数平滑法，预测模型为\(\hat{y_{t+1}}=S_t^{(1)}=\alpha y_t+(1-\alpha)\hat{y_t}\)

  1. 加权系数的选择

    加权系数\(\alpha\)的大小决定了在新预测值中新数据和原预测值所占的比重，\(\alpha\)越大，新数据所占比重越大，原预测值所占比重越小，也可将上式改为\(\hat{t_{t+1}}=\hat{y_t}+\alpha(y_t-\hat{y_t})\)

    这可看出新预测值是根据预测误差对原预测值进行修正得到的，\(\alpha\)越大，修正幅度越大 其选取规则如下：

    • 如果时间序列波动不大，比较平稳，则α可取小一点以减少修正幅度，使预测模型能包含较长时间序列的信息
    • 如果时间序列具有迅速且明显的变动倾向，则α应取大一点，使预测模型灵敏度高一些，以便迅速跟上数据的变化

  2. 初始值的确定

    用一次平滑法预测，除了选择合适的\(\alpha\)，还要确定初始值\(s_0^{(1)}\),初始值是预测者估计指定的

二次指数平滑法 一次指数平滑法虽然克服了移动平均法的缺点，但当时间序列变动出现直线趋势时，用该方法预测由明显滞后偏差，应用二次指数平滑，用滞后偏差的规律建立直线趋势模型，计算公式如下 \(\begin{cases} S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}\\ S_t^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_{t-1}^{(2)} \end{cases}\) 式中：\(S_t^{(1)}\)为一次指数的平滑值，\(S_t^{(2)}\)为二次平滑函数的平滑值

当时间序列\({y_t}\)从某时期开始具有直线趋势，可用直线趋势模型 \(\hat{y_{t+m}}=a_t+b_tm,m=1,2,\dots,\) \(\begin{cases} a_t=2S_t^{(1)}-S_t^{(2)},\\ b_t=\frac{\alpha}{1-\alpha}(S_t^{(1)}-S_t^{(2)}) \end{cases}\)进行预测

三次指数平滑法 当时间序列的变动表现为二次曲线趋势，需要用三次指数平滑法，其在二次指数平滑基础上，再进行一次平滑，计算公式为 \(\begin{cases} S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}\\ S_t^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_{t-1}^{(2)}\\ S_t^{(3)}=\alpha S_t^{(2)}+(1-\alpha)S_{t-1}^{(3)} \end{cases}\) \(S_t^{(3)}\)为三次指数平滑值

其预测模型为: \(\hat{y_{t+m}}=a_t+b_tm+C_tm^2,m=1,2,\dots,\) \(\begin{cases} a_t=3S_t^{(1)}-3S_t^{(2)}+S_t^{(3)},\\ b_t=\frac{\alpha}{2(1-\alpha)^2}((6-5\alpha)S_t^{(1)}-2(5-4\alpha)S_t^{(2)}+(4-3\alpha)S_t^{(3)})\\ c_t=\frac{\alpha^2}{2(1-\alpha)^2}(S_t^{(1)}-2S_t^{(2)}+S_t^{(3)}) \end{cases}\)

差分指数平滑法

我们也可以从数据变换角度解决一次指数平滑法的滞后偏差问题，差分方法就是改变数据变动趋势的简易方法

下面以一阶差分指数平滑法为例，其公式如下 \(\nabla y_t=y_t-y_{t-1} \\ \nabla \hat{y_{t+1}}=\alpha \nabla y_t+(1-\alpha)\nabla\hat{y_t} \\\hat{y_{t+1}}=\nabla\hat{y_{t+1}}+y_t,\nabla\)为差分记号

具有季节性特点的时间序列的预测

对于季节性时间序列的预测，要在数学上完全拟合其变化曲线是非常困难的，但预测的目的是为了找到时间序列的变化趋势，尽可能的做到精确，下面介绍季节系数法，计算步骤如下 1. 收集m年的每年各季度或各月份的时间序列样本数据\(a_{ij}\),i表示年份序号，j表示季度(月份)序号 2. 计算每年所有季度或月份算术平均值\(\bar{a}=\frac{1}{k}\sum_{i=1}^m\sum_{j=1}^n a_{ij},k=mn\) 3. 计算同季度或同月份数据的算术平均值\(\bar{a_{.j}}=\frac{1}{m}\sum_{i=1}^ma_{ij}\) 4. 计算季度(月份)系数\(b_j=\frac{\bar{a_{.j}}}{\bar{a}}\) 5. 预测计算，当时间序列按季度列出，先求出预测年份(下一年)的年加权平均\(y_{m+1}=\frac{\sum_{i=1}^mw_iy_j}{\sum_{i=1}^mw_i} \\y_i=\sum_{j=1}^na_{ij}\)为第i年年合计数，\(w_i\)为第i年权数,\(w_i=i\),再计算预测年份的季度平均值\(\bar{y_{m+1}}=y_{m+1}/n\),最后预测年份第j季度的预测值为\(y_{m+1,j}=b_j\bar{y_{m+1}}\)

平稳时间序列模型

这里的平稳指的是宽平稳，其特性是序列的统计特性不随时间Pinyin而变化，即均值和协方差不随时间平移而变化

时间序列的基本概念

平稳序列:设随机序列\({X_t,t=0,\pm1,\pm2,...}\)满足 - \(E(X_t)=\mu\)=常数 - \(\gamma_{t+k,t}=\gamma_k(k=0,\pm1,\pm2,...)\)与t无关，则称\(X_t\)为平稳随机序列,坚持平稳序列

平稳白噪声序列:设平稳序列\({\epsilon_t,t=0,\pm1,\pm2,...}\)的自协方差函数为\(\gamma_k=\sigma^2\delta_{k,0}=\begin{cases} 0,k\ne0\\ \sigma^2,k=0 \end{cases}\) 则称该序列为平稳白噪声序列

偏相关函数:考虑由\({X_{t-1},X_{t-2},...,X_{t-k}}\)对\(X_t\)的线性最小均方估计，即选择系数\(\psi_{k,1},\psi_{k,2},...,\psi_{k,k}\),使得\(min\quad\delta=E[(X_t-\sum_{j=1}^k\psi_{k,j}X_{t-j})^2]\),\({\psi_{k,k},k=1,2,\dots}\)称为\(X_t\)的偏相关函数

ARMA时间序列:ARMA时间序列可以分为三种类型 1. AR序列，即自回归序列(Auto Regressive Model) 设\({X_t,t=0,\pm1,\pm2,...}\)为零均值平稳序列，满足下列模型 \(X_t=\psi_1 X_{t-1}+\psi_2 X_{t-2}+\dots+\psi_p X_{t-p}+\epsilon_t\) \(\epsilon_t\)为零均值，方差为\(\delta_\epsilon^2\)的平稳白噪声，\(X_t\)为阶数为p的自回归序列，简称AR(p)序列，\(\psi=[\psi_1,\psi_2,...,\psi_p]^T\)为自回归参数向量,其分量为自回归系数 引入后移算子B辅助定义，算子B定义\(BX_t\equiv X_{t-1},B^kX_t\equiv X_{t-k}\) 记算子多项式\(\psi(B)=1-\psi_1B-\psi_2B^2-...-\psi_pB^p\) 则AR序列可改写为\(\psi(B)X_t=\epsilon_t\)

  1. MA序列，即移动平均序列(Moving Average Model) 设\({X_t,t=0,\pm1,\pm2,...}\)为零均值平稳序列，满足下列模型 \(X_t=\epsilon_t-\theta_1 \epsilon_{t-1}-\dots-\theta_q \epsilon_{t-q}\) \(\epsilon_t\)为零均值，方差为\(\delta_\epsilon^2\)的平稳白噪声，\(X_t\)为阶数为q的移动平均序列，简称MA(q)序列，\(\theta=[\theta_1,\theta_2,...,\theta_p]^T\)为移动平均参数向量,其分量为移动平均系数，引入线性后移算子B，再定义一个和上面类似多项式,则可MA序列改写为\(X_t=\theta(B)\epsilon_t\)

  2. ARMA序列，即自回归移动平均序列(Auto Regressive Moving Model) 设\({X_t,t=0,\pm1,\pm2,...}\)为零均值平稳序列，满足下列模型 \(X_t-\psi_1 X_{t-1}-\psi_2 X_{t-2}-\dots-\psi_p X_{t-p}=\epsilon_t-\theta_1 \epsilon_{t-1}-\dots-\theta_q \epsilon_{t-q}\) \(\epsilon_t\)为零均值，方差为\(\delta_\epsilon^2\)的平稳白噪声，\(X_t\)为阶数为p,q的移动平均序列，简称ARMA(p,q)序列,利用后移算子可表为\(\psi(B)(X_t-\mu)-\theta(B)\epsilon_t\)

    对于算子多项式\(\psi(B),\theta(B)\)通常还要做下列假定

    1. 二者无公共因子，且\(\psi_p\ne0,\theta_q\ne0\)
    2. \(\psi(B)=0\)的根全在单位圆外，称为模型平稳性条件
    3. \(\theta(B)=0\)的根全在单位圆外，称为模型可逆性条件

ARMA模型的构建和预报

在实际问题建模中，首先要进行模型的识别和定阶，即要判断AR(p),MA(q),ARMA(p,q)模型的类别，并估计阶数p,q，其实都归结到模型的定阶问题，当模型定阶后，就要对模型参数\(\psi,\theta\)进行估计，完成后还需要对模型进行检验，检验\(\epsilon_t\)是否为平稳白噪声

  1. ARMA模型的构建
    • ARMA模型定阶的AIC准则：选p,q使\(min \quad AIC-nln\hat{\delta_\epsilon^2}+2(p+q+1)\),若序列含未知均值参数，修正为2(p+q+2)
    • ARMA模型的参数估计：直接用Matalb工具箱即可
    • ARMA模型检验的\(\chi^2\)检验：给定显著性水平\(\alpha\)，查表得上\(\alpha\)分位数\(\chi_{\alpha}^2(L-r)\),当\(\chi^2>chi_{\alpha}^2(L)\)认为模型检验未通过
  2. ARMA序列的预报 时间序列的m步预报，是根据\({X_k,X_{k-1},...}\)的取值对未来k+m时刻的随机变量\(X_{k+m}(m>0)\)做出估计，估计量记作\(\hat{X_k}(m)\)，是\({X_k,X_{k-1},...}\)的线性组合
    • AR(p)序列的预报:\(\begin{cases} \hat{X_k(1)}=\psi_1X_k+\psi_2X_{k-1}+\dots+\psi_pX_{k-p+1}\\ \hat{X_k(2)}=\psi_1X_k(1)+\psi_2X_{k}+\dots+\psi_pX_{k-p+2}\\ ...\\ \hat{X_k(p)}=\psi_1X_k(p-1)+\psi_2X_{k}(p-2)+\dots+\psi_pX_{k}\\ \hat{X_k(m)}=\psi_1X_k(m-1)+\psi_2X_{k}(m-2)+\dots+\psi_p\hat{X_k(m-p)},m>p\\ \end{cases}\) 由此可见\(\hat{X_K}(m)(m\ge1)\)仅依赖于\(X_t\)的k时刻以前的p个时刻的值\(X_K,X_{k-1},\dots,X_{k-p+1}\),这是AR(p)序列预报的特点
    • MA(q)和ARMA(p,q)序列的预报 MA(q)序列 \(\hat{X_{k+1}^{(q)}}=\begin{bmatrix} \theta_1 &1&0&\cdots&0\\ \theta_2 &0&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \theta_{q-1} &0&0&\cdots&1\\ \theta_q &0&0&\cdots&0\\ \end{bmatrix}\hat{X_{k}^{(q)}}-\begin{bmatrix} \theta_1\\ \theta_2\\ \vdots\\ \theta_1 \end{bmatrix}X_{k+1}\) 递推初值可取\(\hat{X_{k_0}^{(q)}}=0(k_0)\)较小，因为模型可逆性保证递推式渐进稳定，即当n充分大后，初始误差影响可以逐渐消失 ARMA(p,q)序列 \[\hat{X_{k+1}^{(q)}}=\begin{bmatrix} -G_1 &1&0&\cdots&0\\ -G_2 &0&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ -G_{q-1} &0&0&\cdots&1\\ -G_q+\psi_q^* &\psi_{q-1}^* &\psi_{q-2}^*&\cdots&\psi_1^*\\ \end{bmatrix}\hat{X_{k}^{(q)}}\\+\begin{bmatrix} G_1\\ G_2\\ \vdots\\ G_q\\ \end{bmatrix}X_{k+1}+\begin{bmatrix} 0\\ 0\\ \vdots\\ 0\\ \sum_{j=q+1}^p\psi_j^*X_{k+q+1-j}\\ \end{bmatrix},G_j:X_t=\sum_{j=0}^∞G_j\epsilon_{t-j}\] 实际中，模型参数是未知的，若已建立了时间序列的模型，则理论模型的未知参数可用其估计替代，再用上面介绍的方法预报

时间序列的Matlab相关工具箱及命令

Matlab时间序列的相关指令在系统辨识，计量经济学，金融工具箱内

ARIMA序列与季节性序列

在实践中遇到的时间序列往往有三个特性，趋势性，季节性，非平稳性，一般采用差分方法或其他时间序列变化方法消除趋势性季节性，使变换后序列为平稳序列，并假设为ARMA序列，再用上面介绍方法研究。

ARIMA序列及其预报

对于非平稳序列，可借助差分运算使其平稳化，若\({X_t,t=0,\pm1,\pm2,...}\)为非平稳序列，若存在正整数d，使\(\nabla^dX_t=W_t\)，且\({W_t,t=0,\pm1,\pm2,...}\)是ARMA(p,q)序列，称\(X_t\)为ARIMA(p,d,q)序列

ARIMA序列的预报(d=1,d=2) 设\({X_t,t=0,\pm1,\pm2,...}\)是ARIMA(p,d,q)序列 1. \(d=1,\nabla X_t=W_t \to\hat{X_k(m)}=\hat{X_k}(m-1)+\hat{W_k}(m)=X_k+\sum_{j=1}^m\hat{W_k}(j)\) 2. \(d=2,\nabla^2 X_t=\hat{W_k}(m) \to\hat{X_k(m)}=X_k+\sum_{j=1}^m\hat{W_k}(j)(m+1-j)+m(X_k-X_{k-1})\)

季节性序列及其预报

由季节性因素或其他周期因素引起的周期性变化的时间序列，称为季节性时间序列，相应的模型为季节性模型，一般的，对周期为s的序列，先进行差分运算\(\nabla_sX_t=(1-B^s)X_t,\\\nabla_s^d=(1-B^s)^dX_t\),然后再进行ARIMA建模

支持向量机

支持向量机是数据挖掘中一项新技术，借助最优化方法解决机器学习问题的新工具，在模式识别等领域获得了广泛应用，其主要思想是找到一个超平面，使得它能够尽可能多地将两类数据点正确分开，同时使分开的两类数据点距离分类面最远

支持向量分类机的基本原理

根据给定的训练集，\(T={[a_1,y_1],[a_2,y_2],\dots,[a_t,y_t]}\in(\Omega\times Y)^l\),\(a_i\in\Omega=R^n,\Omega\)为输入空间，输入空间的每一个点\(a_i\)由n个属性特征组成，\(y_i\in Y={-1,1},i=1,\dots,l\) 寻找\(R^n\)上的一个实值函数\(g(x)\),以便用分类函数\(f(x)=sgn(g(x))\)推断任意一个模式x所对应的y值的问题为分类问题

线性可分支持向量分类机

考虑训练集T,若\(\exists \omega\in R^n,b\in R,\epsilon(>0)\)，使得对所有使\(y_i=1\)的\(a_i\)有\((\omega\cdots a_i)+b\ge\epsilon\),这里\((\omega\cdots a_i)\)表示向量\(\omega,a_i\)的内积，而对所有使\(y_i=-1\)的\(a_i\)有\((\omega\cdots a_i)+b\le-\epsilon\),则称训练集T线性可分，称相应的分类问题是线性可分的

规范超平面：空间\(R^n\)中超平面都可以写为\((\omega\cdot x)+b=0\)的形式，参数\((w,b)\)乘以任何一个非零参数后得到的是同一个超平面，定义满足条件\(\begin{cases}y_i[(\omega\cdot a_i+b)]\ge 0,\\min_{i=1,...,l}|(\omega\cdot a_i)+b|=1,i=1,...,l\\\end{cases}\)的超平面为训练集T的规范超平面，

当训练集T为线性可分时，存在唯一的规范超平面\((\omega\cdot x)+b=0\),使得\(\begin{cases}(\omega\cdot a_i)+b\ge 1,y_i=1\\(\omega\cdot a_i)+b\le -1,y_i=-1 \end{cases}\)，其中=正负1的被称为普通支持向量

支持向量机一般又可分为以下几类 - 线性支持向量机：除了普通支持向量分布在两个边界，其余的所有样本点分布在分类边界以外，这是构造的超平面是硬间隔超平面，若存在样本点在边界之间，需要进行软化处理，这是获得的就是软间隔超平面，软化方法就是加入一个调节因子

  • 非线性可分支持向量机：当两个凸包重合太大了，也就是软化的方法不适用的时候，采用映射的方法，换到另外一个高维空间域里进行分类。通过引进从输入空间 X 到另一个高维的 Hilbert 空间 H 的变换，空间H叫做特征空间。

  • C-支持向量机（非线性不可分）：映射到高维H空间之后还是不能直接可分，在H空间进行软化再分类。

该模型主要是利用图像学原理，通过对凸包的研究把求超平面的问题转化成求约束下的优化问题，从而利用拉格朗日算子和KTT条件来求解。

支持向量机的Matlab命令

Matlab支持向量机的命令有，训练支持向量机分类器的函数svmtrain,使用支持向量机分类的函数svmclassify,指定支持向量机函数使用的序列最小化参数函数svmsmoset

多元分析

多元分析是多变量的统计分析方法，是数理统计中应用广泛的一个重要分支

聚类分析

将认识对象进行分类是人类认识世界的一种重要方法，聚类分析作为一种定量方法，可以从数据分析的角度，给出一个更精确，细致的分类工具，聚类分析又被称为群分析，是对多个样本(或指标)进行定量分类的一种多元统计分析方法，对样本进行分类称为Q型聚类分析，对指标进行分类称为R型聚类分析

Q型聚类分析

  1. 样本的相似性度量 要用数量化方法对事物进行分类，就必须用数量化方法描述事物之间的相似程度，对每个样本点的相似程度，一般利用距离来衡量，在这里距离的选择有很多种，最常用的是闵氏距离，绝对值距离，欧几里得距离，切比雪夫距离都是其的一些特例，其中最常用的还是欧几里得距离，其主要优点是当坐标轴进行正交旋转时，欧氏距离是保持不变的，因此，若对原坐标系进行平移和旋转变换，则变换后样本点的距离和变换前完全相同，闵氏距离也有一些缺点，需要量纲一致，需要避免变量的多重相关性 针对其缺点，马氏距离做出了相应的改进

  2. 类与类间的相似性度量 度量方法如下：

    • 最短距离法
    • 最长距离法
    • 重心法
    • 类平均法
    • 离差平方和法

  3. 聚类图 Q型聚类结果可由一个聚类图表示出来，生成步骤如下

    1. 计算n个样本点两两之间的距离\({d_{ij}}\),记为矩阵\(D=(d_{ij})_{n\times n}\)
    2. 首先构造n个类，每一个类只包含一个样本点，每一类的平台高度均为0
    3. 合并距离最近的两类为新类，并且以这两类的距离值作为聚类图的平台高度
    4. 计算新类与当前各类的距离，若类的个数已经等于1，转入步骤5，否则回到步骤3
    5. 画聚类图
    6. 决定类的个数和类

  4. Matlab聚类分析的相关指令 pdist,linkage,cluster,zsore,dendrogram,clusterdata

R型聚类法

在实际工作中，变量聚类法的应用非常重要，人们希望能研究变量间的相似关系，按照变量的相似关系把它们聚合成若干个类，进而找出影响系统的主要因素 1. 变量相似性度量 在对变量进行聚类分析时，首先要确定变量的相似性度量，常用的变量相似性度量有相关系数和夹角余弦 各种定义的相似度量均需具有以下两个性质 1. \(|r_{jk}\le 1|,r_{jk}=r_{kj}\)，越接近1，二者越相关，越接近0，相似性越弱 2. 变量聚类法 类似系统聚类法，在变量聚类中常用的是最长距离法和最短距离法

主成分分析

主成分分析的主要目的是希望用较少的变量去解释原来资料中大部分变异，通常是选出比元素变量个数少，能解释大部分资料中变异的几个新变量，即所谓的主成分，并用以解释资料的综合性指标，本质上是一种降维方法

分析步骤

  1. 构建初始数据矩阵：矩阵的每一行表示一个样本，每一列表示一个原始指标，矩阵中每一个元素表示某一个样本在某一个指标下的得分
  2. 计算相关系数矩阵：对初始数据矩阵中的各个指标，求出相关系数矩阵
  3. 计算特征值和特征向量：计算样本相关系数矩阵的特征值和特征向量，计算完成和对特征值进行从大到小的排序
  4. 计算主成分贡献率和累积贡献率：每一个特征值对应一个主成分和贡献率，某特征值的贡献率=该特征值/所有特征值的和，该特征值对应的主成分即为对应的特征向量
  5. 选择保留的主成分：一般取累积贡献率超过80%的那些主成分作为最终结果，特征值对应的特征向量即为主成分的系数
  6. 利用系数分析主成分代表的含义
  7. 利用主成分进行后续的分析过程

不足与注意事项

不足：主成分的含义解释一般都带有一些模糊性，不像原始变量那样清楚确切，这是变量降维过程中必须付出的代价，只有当主成分个数远小于原始变量的个数才使用主成分分析

注意事项： 在评价类模型中，不能先用主成分分析后再对主成分进行评价。

主成分分析可以用于聚类。因为聚类结果的维度往往很高，这时可以通过主成分分析法对聚类结果降维，从而在二维或三维空间中作出聚类结果图

因子分析

因子分析可以看为主成分分析的推广，也是多元统计分析中常用的一种降维方式，他们的差别体现在： 1. 主成分分析将方差划分为不同的正交成分，而因子分析把方差划分为不同的起因因子 2. 主成分分析只是变量变换，因子分析需要构造因子模型 3. 主成分分析原始变量的线性组合表示新的综合变量，即主成分。因子分析潜在的假想变量和随机影响变量的线性组合表示原始变量

因子分析模型

数学模型：\(\mathbf{X-\mu=\Lambda F+\epsilon}\\\mathbf{X,\mu,\Lambda,F,\epsilon}=\begin{bmatrix}X_1\\X_2\\\vdots\\X_p\end{bmatrix},\begin{bmatrix}\mu_1\\\mu_2\\\vdots\\\mu_p\end{bmatrix},\begin{bmatrix}\alpha_{11}\quad\alpha_{12}\quad\cdots\quad\alpha_{1m}\\\alpha_{21}\quad\alpha_{22}\quad\cdots\quad\alpha_{2m}\\\vdots\qquad\vdots\quad\ddots\quad\vdots\\\alpha_{p1}\quad\alpha_{p2}\quad\cdots\quad\alpha_{pm}\end{bmatrix},\begin{bmatrix}F_1\\F_2\\\vdots\\F_p\end{bmatrix},\begin{bmatrix}\epsilon_1\\\epsilon_2\\\vdots\\\epsilon_p\end{bmatrix}\),在上式中 称\(F_1,F_2,...,F_p\)为公共因子，是不可观测的变量，其系数称为载荷因子，\(\epsilon_i\)是特殊因子，不能被前m个公共因子包含的部分

因子载荷矩阵的估计方法

用于估计\(\mathbf{\Lambda}\) 1. 主成分分析法 2. 主因子法 3. 最大似然估计法

因子旋转(正交变换)

对因子载荷矩阵进行旋转的目的是使因子载荷矩阵结构简化，使载荷矩阵每列或行的元素平方值向0和1两极分化，主要正交旋转法有：方差最大法，四次方最大法，等量最大法

因子得分

上述小节解决了用公共因子的线性组合表示一组观测变量的问题，如果利用因子做其他研究，比如作为自变量来做回归分析，对样本分类或平均，需要对公共因子测度，给出公共因子的值，其数学模型就是将上述数学模型的\(\mathbf{\mu}\)移到等式右边 常用方法有: 巴斯莱特因子得分(加权最小二乘法),回归方法

因子分析的步骤： 1. 选择分析的变量 2. 计算所选的原始变量的相关系数矩阵 3. 提出公共因子 4. 因子旋转 5. 计算因子得分

判别分析

判别分析是根据所研究的个体的观测指标来推断该个体所属类型的一种统计方法，用统计方法表述，就是已有q个总体\(X_1,X_2,...,X_q\),他们的分布函数为\(F_1(x),F_2(x),...,F_q(x)\)，对于给定样本X,要判断它来自哪个总体