现代优化算法 之 禁忌搜索算法
最后更新于：2022-04-01 16:07:13
这次是首次接触这个算法，看了一些资料，总结一下。
### 禁忌搜索算法简介
禁忌搜索算法是组合优化算法的一种，是局部搜索算法的扩展。禁忌搜索算法是人工智能在组合优化算法中的一个成功应用。禁忌搜索算法的特点是采用了禁忌技术。所谓禁忌就是禁止重复前面的工作。禁忌搜索算法用一个禁忌表记录下已经到达过的局部最优点，在下一次搜索中，利用禁忌表中的信息不再或有选择地搜索这些点。
禁忌搜索算法实现的技术问题是算法的关键。禁忌搜索算法涉及侯选集合、禁忌对象、评价函数、特赦规则、记忆频率信息等概念。
### 邻域
在组合优化中，距离的概念通常不再适用，但是在一点附近搜索另一个下降的点仍然是组合优化数值求解的基本思想。因此，需要重新定义邻域的概念。
一个组合优化问题可用三参数(D,F,f)表示，其中D表示决策变量的定义域，F表示可行解区域， F中的任何一个元素称为该问题的可行解， f表示目标函数。满足
f(x∗)=min{f(x)|x∈F}的可行解x∗称为该问题的最优解。
定义 1 对于组合优化问题(D,F,f)，D上的一个映射：
N:S∈D→N(S)∈2D
称为一个邻域映射，其中
2D
表示
D
的所有子集组成的集合，
N(S)
称为
S
的邻域，
S′∈N(S)
称为
S
的一个邻居。
### 侯选集合
侯选集合由邻域中的邻居组成。常规的方法是从邻域中选择若干个目标值或评价值最佳的邻居入选。
### 禁忌对象和禁忌长度
禁忌表中的两个主要指标是禁忌对象和禁忌长度。禁忌算法中，由于我们要避免一些操作的重复进行，就要将一些元素放到禁忌表中以禁止对这些元素进行操作，这些元素就是我们指的禁忌对象。禁忌长度是被禁对象不允许选取的迭代次数。一般是给被禁对象x一个数（禁忌长度） t ，要求对象x 在t 步迭代内被禁，在禁忌表中采用tabu(x)=t记忆，每迭代一步，该项指标做运算tabu(x)=t−1，直到tabu(x)=0时解禁。于是，我们可将所有元素分成两类，被禁元素和自由元素。禁忌长度t 的选取可以有多种方法，例如t=常数，或t=[n√]，其中n为邻域中邻居的个数；这种规则容易在算法中实现。
### 评价函数
评价函数是侯选集合元素选取的一个评价公式，侯选集合的元素通过评价函数值来选取。以目标函数作为评价函数是比较容易理解的。目标值是一个非常直观的指标，但有时为了方便或易于计算，会采用其他函数来取代目标函数。
### 特赦规则
在禁忌搜索算法的迭代过程中，会出现侯选集中的全部对象都被禁忌，或有一对象被禁，但若解禁则其目标值将有非常大的下降情况。在这样的情况下，为了达到全局最优，我们会让一些禁忌对象重新可选。这种方法称为特赦，相应的规则称为特赦规则。
### 记忆频率信息
在计算的过程中，记忆一些信息对解决问题是有利的。如一个最好的目标值出现的频率很高，这使我们有理由推测：现有参数的算法可能无法再得到更好的解。根据解决问题的需要，我们可以记忆解集合、被禁对象组、目标值集合等的出现频率。
频率信息有助于进一步加强禁忌搜索的效率。我们可以根据频率信息动态控制禁忌的长度。一个最佳的目标值出现的频率很高，有理由终止计算而将此值认为是最优值。
### 模型及求解
我们用禁忌搜索算法研究如下的两个问题：
（1）研究前几篇现代优化算法解决的问题。
（2）我方有三个基地，经度、纬度分别为（70,40）,（72,45），（68,48）。假设我方所有无人侦察机的速度都为1000 公里/小时。三个基地各派出一架飞机侦察敌方目标，怎样划分任务，才能使时间最短，且任务比较均衡。
### 问题（1）的求解
求解的禁忌搜索算法描述如下：
（1）解空间
解空间S 可表为{1,2,…,101,102}的所有固定起点和终点的循环排列集合，即
S={(π1,…,π102)|π1=1,(π2,…,π101)为2,3,…,101}的循环排列,，π102=102。其中每一个循环排列表示侦察 100个目标的一个回路， πij = 表示第i次侦察 j点。
（2）目标函数
目标函数为侦察所有目标的路径长度。我们要求
min f(π1,π2,…,π102)=∑i=1101wπiπi+1
（3）候选集合
定义2 任选序号u,v（u<v）交换u与v之间的顺序，此时的新路径为：
π1…πu−1πvπv+1…πu+1πuπv+1…π102
称为原路径二邻域的一个邻居。
定义 3 任选序号u,v 和 w，将u 和v之间的路径插到 w 之后，（设u<v<w）对应的新路径为：
π1…πu−1πv+1…πwπu…πvπw+1…π102
称为原路径三邻域的一个邻居。
如果要考虑当前解的全部二邻域（或三邻域）的邻居，将面临着太大的工作量。因此我们用随机选取的方法每次选取50 个邻居组成的集合作为侯选集合。而将省下的时间作更多次搜索，这样做同样可以保证较高的精确度，同时可以大大提高算法的效率。
（4）禁忌长度及禁忌对象
对 于 上述定义的二邻域中的邻居总数为C2100， 我们的禁忌长度取为
t=70≈C2100−−−−√ 。
我们把禁忌表设计成一个循环队列，初始化禁忌表H=Φ 。从候选集合C中选出一个向量x⃗ ，如果x⃗ ∉H ，并且H 不满，则把向量x添加到禁忌表中；如果H已满，则最早进入禁忌表的向量出列，向量x⃗ 进入到出列的位置。
（5）评价函数(类似启发式搜索的那种，大家有兴趣可以百度一下算法A∗，对评估函数的设计会更加理解)
可以用目标函数作为评价函数，但是这样每选取一个新的路径都得去计算总时间，计算量比较大。对于上述二邻域中的邻居作为侯选集合，每一个新路径中只有两条边发生了变化，因此将目标函数的差值作为评价函数可以极大地提高算法的效率。评价函数取为
Δf=(wπu−1πv+wπuπv+1)−(wπu−1πu+wπvπv+1)
**禁忌搜索算法的流程如下：**
STEP1:
选定一个初始解xnow及给以禁忌表H=Φ；
STEP2:
若满足停止条件，停止计算；否则，在xnow的邻域N(H,xnow)中选出满足禁忌要求的侯选集Can_N(xnow )，在Can _ N(xnow ) 中选一个评价值最佳的解xnext，xnow:=xnext，更新禁忌表H ，重复 STEP2。
利用 Matlab 程序求得，我们的巡航时间大约在41 小时左右，其中的一个巡航路径如下图所示
### 问题（2）的求解
对于这个问题，我们的基本想法是，先根据敌方基地的分布特点将敌方的基地大体划分在三个区域之内，并使三架侦察机分别对这三个区域的敌军基地进行侦察，求取各自的最短时间。然后对任务不均衡区域之中的点做适当调整。
我们解决问题的步骤如下：
（1）划分子图。要达到比较均衡，应使每架飞机的巡航时间基本相同，由于敌方目标的分布较均匀，可以将敌方目标的地理位置图分成面积基本相同的三部分。如过点（70，40）以斜率k1=45作一条斜线，过点（68，48）以斜率k2=4368作一个斜率，把基地所在的地区划分成三部分G1,G2,G3。
（2）再对每个子图Gi(i=1,2,3) 分别运用禁忌搜索算法求得其最短侦察路径ψi，和最短侦察时间 t(ψi)(i=1,2,3)。
（3）均衡任务。定义均衡率
η=min{t(ψ1),t(ψ2),t(ψ3)}max{t(ψ1),t(ψ2),t(ψ3)}×100%
若η接近于1，则上面划分的任务就可以接受。否则的话，根据t(ψi)(i=1,2,3)的大小用局部搜索算法，调整k1,k2的斜率，从而调整各分区内点的个数，直至任务达到均衡。
最后计算结果如下
两直线的斜率分别为k1=0.68,k2=0.6484。
均衡率η=97.57%
三架飞机侦察完所有目标所用的时间分别为
t(ψ1)=20.401小时，t(ψ1)=20.910小时，t(ψ1)=20.729小时
现代优化算法 之 遗传算法
最后更新于：2022-04-01 16:07:10
之前两篇转载的文章：
[遗传算法入门到掌握（一）](http://blog.csdn.net/u013007900/article/details/43523975)、[遗传算法入门到掌握（二）](http://blog.csdn.net/u013007900/article/details/43524005)
对遗传算法的数学推导讲解得非常详细，同时我也附带了一份[遗传算法的C语言实现](http://blog.csdn.net/u013007900/article/details/43524261)，这篇文章将要运用遗传算法对一个多项式求最小值，要求在(-8,8)间寻找使表达式达到最小的x，误差为0.001。
但是那篇文章仅仅讲解了关于本例的遗传算法的规则，并没有详细的算法过程。
这篇文章简介一下遗传算法的具体算法过程，并且用MATLAB实现遗传算法的代码，该算法将解决[模拟退火](http://blog.csdn.net/u013007900/article/details/50351523)一文中的例题。
### 遗传算法简介
遗传算法（Genetic Algorithms，简称 GA）是一种基于自然选择原理和自然遗传机制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目标的优化。遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终得到最优解或准最优解。它必须做以下操作：初始群体的产生、求每一个体的适应度、根据适者生存的原则选择优良个体、被选出的优良个体两两配对，通过随机交叉其染色 体的基因并随机变异某些染色体的基因后生成下一代群体，按此方法使群体逐代进化，直到满足进化终止条件。其实现方法如下：
- （1） 根据具体问题确定可行解域，确定一种编码方法，能用数值串或字符串表示 可行解域的每一解。
- （2） 对每一解应有一个度量好坏的依据，它用一函数表示，叫做适应度函数，适 应度函数应为非负函数。
- （3） 确定进化参数群体规模M 、交叉概率pc 、变异概率pc 、进化终止条件。
为便于计算，一般来说，每一代群体的个体数目都取相等。群体规模越大、越容易找到最优解，但由于受到计算机的运算能力的限制，群体规模越大，计算所需要的时间也相应的增加。进化终止条件指的是当进化到什么时候结束，它可以设定到某一代进化结束，也可能根据找出近似最优是否满足精度要求来确定。表1 列出了生物遗传概念在遗传算法中的对应关系。
表1 生物遗传概念在遗传算法中的对应关
| 生物遗传概念 | 遗传算法中的作用 |
|-----|-----|
| 适者生存 | 算法停止时，最优目标值的解有最大的可能被留住 |
| 个体 | 解 |
| 染色体 | 解的编码 |
| 基因 | 解中每一分量的特征 |
| 适应性 | 适应度函数值 |
| 种群 | 根据适应度函数值选取的一组解 |
| 交配 | 通过交配原则产生一组新解的过程 |
| 变异 | 编码的某一分量发生变化的过程 |
### 模型及算法
我们用遗传算法研究上次提到的问题。
- 种群大小：M=50
- 最大代数：G=1000
- 交叉率： pc=1 ，交叉概率为1 能保证种群的充分进化。
- 变异率： pm=0.1 ， 一般而言，变异发生的可能性较小。
（1） 编码策略
采用十进制编码，用随机数列ω1,ω2,…,ω102 作为染色体，其中 0<ωi<1（i=2,3,…,101）， ω1=0，ω102=1 ；每一个随机序列都和种群中的一个个体相对应，例如一个9 城市问题的一个染色体为
[0.23，0.82，0.45，0.74，0.87，0.11，0.56，0.69，0.78]
其中编码位置i代表城市i，位置i的随机数表示城市i在巡回中的顺序，我们将这些随 机数按升序排列得到如下巡回：
6→1→3→7→8→4→9→2→5
（2） 初始种群 本文中我们先利用经典的近似算法—改良圈算法求得一个较好的初始种群。即对于初始圈C=π1…πu−1πuπu+1…πv−1πvπv+1…π102，2≤u<v≤101，2≤πu<πv≤101，交换u与v之间的顺序，此时的新路径为：
π1…πu−1πvπv+1…πu+1πuπv+1…π102
记Δf=(dπu−1πv+dπuπv+1)−(dπu−1πu+dπvπv+1)，若Δf<0，则以新的路经修改旧的路经，直到不能修改为止。
（3） 目标函数
目标函数为侦察所有目标的路径长度，适应度函数就取为目标函数。我们要求
minf(π1,π2,…,π102)=∑i=1101dπiπi+1
（4） 交叉操作
我 们 的 交 叉 操 作 采 用 单 点 交 叉 。 设 计 如 下 ， 对 于 选 定 的 两 个 父 代 个 体f1=ω1ω2…ω102，f2=ω1ω2…ω102，我们随机地选取第t个基因处为交叉点，则经过交叉运算后得到的子代编码为s1 和s1， s1 的基因由f1的前t个基因和f2 的后102−t个基因构成， s2的基因由f2 的前t个基因和f1的后102−t个基因构成，例如：
f1=[0,0.14,0.25,0.27, ∥ 0.29,0.54,…,0.19,1]f2=[0,0.23,0.44,0.56, ∥ 0.74,0.21,…,0.24,1]
设交叉点为第四个基因处，则
s1=[0,0.14,0.25,0.27, ∥ 0.74,0.21,…,0.24,1]s2=[0,0.23,0.44,0.56, ∥ 0.29,0.54,…,0.19,1]
交叉操作的方式有很多种选择，我们应该尽可能选取好的交叉方式，保证子代能继承父代的优良特性。同时这里的交叉操作也蕴含了变异操作。
（5） 变异操作
变异也是实现群体多样性的一种手段，同时也是全局寻优的保证。具体设计如下，按照给定的变异率，对选定变异的个体，随机地取三个整数，满足1<u<v<w<102，
把u,v之间（包括u和v）的基因段插到w后面。
（6） 选择
采用确定性的选择策略，也就是说选择目标函数值最小的M个个体进化到下一代，这样可以保证父代的优良特性被保存下来。
### MATLAB 程序
~~~
clc,clear
load sj.txt %加载敌方100 个目标的数据
x=sj(:,1:2:8);x=x(:);
y=sj(:,2:2:8);y=y(:);
sj=[x y];
d1=[70,40];
sj0=[d1;sj;d1];
%距离矩阵d
sj=sj0*pi/180;
d=zeros(102);
for i=1:101
for j=i+1:102
temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));
d(i,j)=6370*acos(temp);
end
end
d=d+d';L=102;w=50;dai=100;
%通过改良圈算法选取优良父代A
for k=1:w
c=randperm(100);
c1=[1,c+1,102];
flag=1;
while flag>0
flag=0;
for m=1:L-3
for n=m+2:L-1
if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<d(c1(m),c1(m+1))+d(c1(n),c1(n+1))
flag=1;
c1(m+1:n)=c1(n:-1:m+1);
end
end
end
end
J(k,c1)=1:102;
end
J=J/102;
J(:,1)=0;J(:,102)=1;
rand('state',sum(clock));
%遗传算法实现过程
A=J;
for k=1:dai %产生0～1 间随机数列进行编码
B=A;
c=randperm(w);
%交配产生子代B
for i=1:2:w
F=2+floor(100*rand(1));
temp=B(c(i),F:102);
B(c(i),F:102)=B(c(i+1),F:102);
B(c(i+1),F:102)=temp;
end
-280-
%变异产生子代C
by=find(rand(1,w)<0.1);
if length(by)==0
by=floor(w*rand(1))+1;
end
C=A(by,:);
L3=length(by);
for j=1:L3
bw=2+floor(100*rand(1,3));
bw=sort(bw);
C(j,:)=C(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]);
end
G=[A;B;C];
TL=size(G,1);
%在父代和子代中选择优良品种作为新的父代
[dd,IX]=sort(G,2);temp(1:TL)=0;
for j=1:TL
for i=1:101
temp(j)=temp(j)+d(IX(j,i),IX(j,i+1));
end
end
[DZ,IZ]=sort(temp);
A=G(IZ(1:w),:);
end
path=IX(IZ(1),:)
long=DZ(1)
toc
xx=sj0(path,1);yy=sj0(path,2);
plot(xx,yy,'-o')
~~~
计算结果为 40 小时左右。其中的一个巡航路径如图2 所示。
现代优化算法 之 模拟退火
最后更新于：2022-04-01 16:07:08
### 现代优化算法
现代优化算法是 80 年代初兴起的启发式算法。这些算法包括禁忌搜索（tabu search），模拟退火（simulated annealing），遗传算法（genetic algorithms），人工神经网络（neural networks）。它们主要用于解决大量的实际应用问题。目前，这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标-求 NP-hard 组合优化问题的全局最优解。虽然有这些目标，但 NP-hard 理论限制它们只能以启发式的算法去求解实际问题。
启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（Ant Colony Algorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。
现代优化算法解决组合优化问题，如 TSP（Traveling Salesman Problem）问题，QAP（Quadratic Assignment Problem）问题，JSP（Job-shop Scheduling Problem）问题等效果很好。
这一章讲解模拟退火的算法过程，之前也介绍过一些简单的[模拟退火的思想](http://blog.csdn.net/u013007900/article/details/43525733)，上次是基于ACM-ICPC的思想进行介绍的，这次是详细的计算推导过程。
### 模拟退火
模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。
如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。假设材料在状态i之下的能量为E(i)，那么材料在温度T时从状态i进入状态j就遵循如下规律：
- （1）如果E(j)≤E(i)，接受该状态被转换。
- （2）如果E(j)>E(i)，则状态转换以如下概率被接受：
eE(i)−E(j)KT
其中K是物理学中的波尔兹曼常数，T是材料温度。
在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态i的概率满足波尔兹曼分布：
PT(x=i)=e−E(i)KT∑j∈Se−E(j)KT
其中 x 表示材料当前状态的随机变量, S 表示状态空间集合。
显然
limT→∞e−E(i)KT∑j∈Se−E(j)KT=1|S|
其中|S|表示集合S中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，
limT→0e−E(i)−EminKT∑j∈Se−E(j)−EminKT=limT→0e−E(i)−EminKT∑j∈Smine−E(j)−EminKT+∑j∉Smine−E(j)−EminKT=limT→0e−E(i)−EminKT∑j∈Smine−E(j)−EminKT=⎧⎩⎨⎪⎪1|Smin|0 i∈Smin other
其中Emin=minj∈SE(j)且Smin=i | E(i)=Emin。
上式表明当温度降至很低时,材料会以很大概率进入最小能量状态。
假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。
考虑这样一个组合优化问题:优化函数为 f:x→R+ ,其中 x∈S ,它表示优化问题的一个可行解, R+=y|y∈R,y>0，S表示函数的定义域。N(x)⊆S表示x 的一个邻域集合。
首先给定一个初始温度T0和该优化问题的一个初始解x(0)，并由x(0)生成下一个解x′∈N(x(0))，是否接受x′作为一个新解x(1)依赖于下面概率:
P(x(0)→x′)=⎧⎩⎨1e−f(x′)−f(x(0))T0 若f(x′)<f(x(0)) other
换句话说，如果生成的解 x’ 的函数值比前一个解的函数值更小,则接受 x(1)=x’ 作为一个新解。
否则以概率 e−f(x′)−f(x(0))T0 T0 接受 x’ 作为一个新解。
泛泛地说,对于某一个温度Ti和该优化问题的一个解x(k), 可以生成x’。接受x’ 作为下一个新解 x(k+1) 的概率为:
P(x(k)→x′)=⎧⎩⎨1e−f(x′)−f(x(k))T0 若f(x′)<f(x(k)) other(1)
在温度Ti下,经过很多次的转移之后,降低温度Ti ,得到Ti+1<Ti 。在Ti+1 下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问 题寻优的结果。
我们注意到,在每个Ti 下,所得到的一个新状态x(k+1)完全依赖于前一个状态 x(k) , 可以和前面的状态 x(0),…,x(k−1) 无关,因此这是一个马尔可夫过程。使用马 尔可夫过程对上述模拟退火的步骤进行分析,结果表明:从任何一个状态 x(k ) 生成 x’ 的 概率,在N(x(k))中是均匀分布的,且新状态x’被接受的概率满足式(1),那么经过有限次的转换,在温度Ti下的平衡态 xi 的分布由下式给出:
Pi(xi)=e−f(x−i)T∑j∈Se−f(xi)Ti
当温度T降为0时, xi 的分布为:
P∗i=⎧⎩⎨⎪⎪1|Smin|0 i∈Smin other
并且
∑xi∈SminPi=1
这说明如果温度下降十分缓慢,而在每个温度都有足够多次的状态转移,使之在每一个 温度下达到热平衡,则全局最优解将以概率 1 被找到。因此可以说模拟退火算法可以找 到全局最优解。
在模拟退火算法中应注意以下问题:
(1)理论上,降温过程要足够缓慢,要使得在每一温度下达到热平衡。但在计算 机实现中,如果降温速度过缓,所得到的解的性能会较为令人满意,但是算法会太慢, 相对于简单的搜索算法不具有明显优势。如果降温速度过快,很可能最终得不到全局最 优解。因此使用时要综合考虑解的性能和算法速度,在两者之间采取一种折衷。
(2)要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m 次的 转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定 为一个较小的值Te ,或连续几个温度下转换过程没有使状态发生变化算法就结束。
(3)选择初始温度和确定某个可行解的邻域的方法也要恰当。
### 应用举例
### 题目
已知敌方100 个目标的经度、纬度如表1 所示。
表1 经度和纬度数据表
| 经度 | 纬度 | 经度 | 纬度 | 经度 | 纬度 | 经度 | 纬度 |
|-----|-----|-----|-----|-----|-----|-----|-----|
| 53.7121 | 15.3046 | 51.1758 | 0.0322 | 46.3253 | 28.2753 | 30.3313 | 6.9348 |
| 56.5432 | 21.4188 | 10.8198 | 16.2529 | 22.7891 | 23.1045 | 10.1584 | 12.4819 |
| 20.1050 | 15.4562 | 1.9451 | 0.2057 | 26.4951 | 22.1221 | 31.4847 | 8.9640 |
| 26.2418 | 18.1760 | 44.0356 | 13.5401 | 28.9836 | 25.9879 | 38.4722 | 20.1731 |
| 28.2694 | 29.0011 | 32.1910 | 5.8699 | 36.4863 | 29.7284 | 0.9718 | 28.1477 |
| 8.9586 | 24.6635 | 16.5618 | 23.6143 | 10.5597 | 15.1178 | 50.2111 | 10.2944 |
| 8.1519 | 9.5325 | 22.1075 | 18.5569 | 0.1215 | 18.8726 | 48.2077 | 16.8889 |
| 31.9499 | 17.6309 | 0.7732 | 0.4656 | 47.4134 | 23.7783 | 41.8671 | 3.5667 |
| 43.5474 | 3.9061 | 53.3524 | 26.7256 | 30.8165 | 13.4595 | 27.7133 | 5.0706 |
| 23.9222 | 7.6306 | 51.9612 | 22.8511 | 12.7938 | 15.7307 | 4.9568 | 8.3669 |
| 21.5051 | 24.0909 | 15.2548 | 27.2111 | 6.2070 | 5.1442 | 49.2430 | 16.7044 |
| 17.1168 | 20.0354 | 34.1688 | 22.7571 | 9.4402 | 3.9200 | 11.5812 | 14.5677 |
| 52.1181 | 0.4088 | 9.5559 | 11.4219 | 24.4509 | 6.5634 | 26.7213 | 28.5667 |
| 37.5848 | 16.8474 | 35.6619 | 9.9333 | 24.4654 | 3.1644 | 0.7775 | 6.9576 |
| 14.4703 | 13.6368 | 19.8660 | 15.1224 | 3.1616 | 4.2428 | 18.5245 | 14.3598 |
| 58.6849 | 27.1485 | 39.5168 | 16.9371 | 56.5089 | 13.7090 | 52.5211 | 15.7957 |
| 38.4300 | 8.4648 | 51.8181 | 23.0159 | 8.9983 | 23.6440 | 50.1156 | 23.7816 |
| 13.7909 | 1.9510 | 34.0574 | 23.3960 | 23.0624 | 8.4319 | 19.9857 | 5.7902 |
| 40.8801 | 14.2978 | 58.8289 | 14.5229 | 18.6635 | 6.7436 | 52.8423 | 27.2880 |
| 39.9494 | 29.5114 | 47.5099 | 24.0664 | 10.1121 | 27.2662 | 28.7812 | 27.6659 |
| 8.0831 | 27.6705 | 9.1556 | 14.1304 | 53.7989 | 0.2199 | 33.6490 | 0.3980 |
| 1.3496 | 16.8359 | 49.9816 | 6.0828 | 19.3635 | 17.6622 | 36.9545 | 23.0265 |
| 15.7320 | 19.5697 | 11.5118 | 17.3884 | 44.0398 | 16.2635 | 39.7139 | 28.4203 |
| 6.9909 | 23.1804 | 38.3392 | 19.9950 | 24.6543 | 19.6057 | 36.9980 | 24.3992 |
| 4.1591 | 3.1853 | 40.1400 | 20.3030 | 23.9876 | 9.4030 | 41.1084 | 27.7149 |
我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000 公里/小时。
我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。
这是一个[旅行商问题](http://baike.baidu.com/link?url=uqnx_j59Qy_JKzQzfdQwnOiN0I3XzFGcxa_GfHci4WG53qVm02KwgZJddmEulnpl7a2TkxSlMzbhoukR_0UdTa)，旅行社问题又是NP完全问题，目前没有已知的算法可以解决。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，最后我方基地再重复编号为 102（这样便于程序中计算）。
距离矩阵D=(dij)102×102，其中dij 表示表示i,j两点的距离，i,j=1,2,…,102，这里D为实对称矩阵。则问题抽象成：
求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。
上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设A, B两点的地理坐标分别为(x1,y1),(x2,y2)，过 A, B两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O，以赤道平面为XOY平面，以0度经线圈所在的平面为XOZ平面建立三维直角坐标系。则 A, B两点的直角坐标分别为：
A(Rcosx1cosy1,Rsinx1cosy1,Rsiny1)B(Rcosx2cosy2,Rsinx2cosy2,Rsiny2)
其中R=6370为地球半径。
A, B两点的实际距离：
d=R arccos(OA→⋅OB→|OA→|⋅|OB→|)
化简得
d=Rarccos[cos(x1−x2)cos y1cos y2+sin y1sin y2]
### 算法描述
求解的模拟退火算法描述如下：
（1）解空间
解空间S可表为1,2,…,101,102的所有固定起点和终点的循环排列集合，即
S={(π1,…,π102)|π1=1,(π2,…,π101为为2,3,…,101的循环排列),π102=102}
其中每一个循环排列表示侦察100个目标的一个回路，πi=j表示在第i次侦察 j点，初始解可选为(1,2,…,102)，本文中我们使用[ Monte Carlo 方法](http://baike.baidu.com/link?url=cW0fyu21-nMWHWfwe0oKhpk25p4ep2Nnk0CwlSz70ihI-OnnKxuHwlN4ej_MvO1saIOcqlKn4-j0L4bGCk0AZq)求得一个较好的初始解。
（2）目标函数
此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求
minf(π1,π2,…,π102)=∑i=1101dπiπi+1
而一次迭代由下列三步构成：
（3）新解的产生
① 2 变换法
任选序号u,v（u<v）交换u与v之间的顺序，此时的新路径为：
π1…πu−1πvπv+1…πu+1πuπv+1…π102
② 3 变换法
任选序号u,v 和 w，将u 和v之间的路径插到 w 之后，（设u<v<w）对应的新路径为：
π1…πu−1πv+1…πwπu…πvπw+1…π102
（4）代价函数差
对于2 变换法，路径差可表示为
Δf=(dπu−1πv+dπuπv+1)−(dπu−1πu+dπvπv+1)
（5）接受准则
P={1exp(−∇f/T) Δf<0 Δf≥0
如果Δf<0，则接受新的路径。否则，以概率exp(−Δf/T)接受新的路径，即若exp(−Δf/T)大于 0 到1之间的随机数则接受。
（6）降温
利用选定的降温系数α 进行降温即：T←αT，得到新的温度，这里我们取
α=0.999。
（7）结束条件
用选定的终止温度e=10−30，判断退火过程是否结束。若T < e$，算法结束，输出当前状态。
MATLAB程序如下:
~~~
clc,clear
load sj.txt %加载敌方100 个目标的数据，数据按照表格中的位置保存在纯文本
文件sj.txt 中
x=sj(:,1:2:8);x=x(:);
y=sj(:,2:2:8);y=y(:);
sj=[x y];
d1=[70,40];
sj=[d1;sj;d1];
sj=sj*pi/180;
%距离矩阵d
d=zeros(102);
for i=1:101
for j=i+1:102
temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));
d(i,j)=6370*acos(temp);
end
end
d=d+d';
S0=[];Sum=inf;
rand('state',sum(clock));
for j=1:1000
S=[1 1+randperm(100),102];
temp=0;
for i=1:101
temp=temp+d(S(i),S(i+1));
end
if temp<Sum
S0=S;Sum=temp;
end
end
e=0.1^30;L=20000;at=0.999;T=1;
%退火过程
for k=1:L
%产生新解
c=2+floor(100*rand(1,2));
c=sort(c);
c1=c(1);c2=c(2);
%计算代价函数值
df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1));
%接受准则
if df<0
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;
elseif exp(-df/T)>rand(1)
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;
end
T=T*at;
if T<e
break;
end
end
% 输出巡航路径及路径长度
S0,Sum
~~~
计算结果为 44 小时左右。其中的一个巡航路径如图所示。
决策树学习 之 概念与构造算法
最后更新于：2022-04-01 16:07:06
### 分类（Classification）
分类预测的是分类标号，什么是分类标号呢？
分类标号是离散值，比如把一个群体分成屌丝、土豪、高富帅三类，分类标号则分别是屌丝、土豪、高富帅。再比如把土豪分为小土豪、大土豪、高级土豪，可以根据月收入来分，5000以下的为小土豪，5000到0000的为大土豪，10000以上的为高级土豪。所以，分类标号需要离散化。
分类是预测的一种，预测的是分类标号，即把没有分类标号的群体分到有分类标号的群体里，也就完成了分类，具体如何分类则需要分类算法。
不管使用什么算法，分类过程主要有两个步骤，如下图表示。
在分类算法中，比较基础的有朴素贝叶斯分类与贝叶斯网络两种分类算法。这两种算法都以贝叶斯定理为基础，可以对分类及决策问题进行概率推断。这个在之后再讨论。
在这一篇文章中，将讨论一种被广泛使用的归纳推理算法——决策树（decision tree）。它是一种逼近离散值函数的方法，对噪声数据有很良好的健壮性且能够学习析取表达式。相比贝叶斯算法，决策树的优势在于构造过程不需要任何领域知识或参数设置，因此在实际应用中，对于探测式的知识发现，决策树更加适用。
### 决策树（Decision Tree）
通俗来说，决策树分类的思想类似于找对象。现想象一个女孩的母亲要给这个女孩介绍男朋友，于是有了下面的对话：
> 女儿：多大年纪了？
母亲：26。
女儿：长的帅不帅？
母亲：挺帅的。
女儿：收入高不？
母亲：不算很高，中等情况。
女儿：是公务员不？
母亲：是，在税务局上班呢。
女儿：那好，我去见见。
这个女孩的决策过程就是典型的分类树决策。相当于通过年龄、长相、收入和是否公务员对将男人分为两个类别：见和不见。假设这个女孩对男人的要求是：30岁以下、长相中等以上并且是高收入者或中等以上收入的公务员，那么这个可以用下图表示女孩的决策逻辑。
上图完整表达了这个女孩决定是否见一个约会对象的策略，其中绿色节点表示判断条件，橙色节点表示决策结果，箭头表示在一个判断条件在不同情况下的决策路径，图中红色箭头表示了上面例子中女孩的决策过程。
这幅图基本可以算是一颗决策树，说它“基本可以算”是因为图中的判定条件没有量化，如收入高中低等等，还不能算是严格意义上的决策树，如果将所有条件量化，则就变成真正的决策树了。
有了上面直观的认识，我们可以正式定义决策树了：
决策树（decision tree）是一个树结构（可以是二叉树或非二叉树）。其每个非叶节点表示一个特征属性上的测试，每个分支代表这个特征属性在某个值域上的输出，而每个叶节点存放一个类别。使用决策树进行决策的过程就是从根节点开始，测试待分类项中相应的特征属性，并按照其值选择输出分支，直到到达叶子节点，将叶子节点存放的类别作为决策结果。
可以看到，决策树的决策过程非常直观，容易被人理解。目前决策树已经成功运用于医学、制造产业、天文学、分支生物学以及商业等诸多领域。知道了决策树的定义以及其应用方法，下面介绍决策树的构造算法。
### 决策树的构造
不同于贝叶斯算法，决策树的构造过程不依赖领域知识，它使用属性选择度量来选择将元组最好地划分成不同的类的属性。所谓决策树的构造就是进行属性选择度量确定各个特征属性之间的拓扑结构。
构造决策树的关键步骤是分裂属性。所谓分裂属性就是在某个节点处按照某一特征属性的不同划分构造不同的分支，其目标是让各个分裂子集尽可能地“纯”。尽可能“纯”就是尽量让一个分裂子集中待分类项属于同一类别。分裂属性分为三种不同的情况：
1、属性是离散值且不要求生成二叉决策树。此时用属性的每一个划分作为一个分支。
2、属性是离散值且要求生成二叉决策树。此时使用属性划分的一个子集进行测试，按照“属于此子集”和“不属于此子集”分成两个分支。
3、属性是连续值。此时确定一个值作为分裂点split_point，按照>split_point和<=split_point生成两个分支。
构造决策树的关键性内容是进行属性选择度量，属性选择度量是一种选择分裂准则，是将给定的类标记的训练集合的数据划分D“最好”地分成个体类的启发式方法，它决定了拓扑结构及分裂点split_point的选择。
属性选择度量算法有很多，一般使用自顶向下递归分治法，并采用不回溯的贪心策略。这里介绍[ID3](http://en.wikipedia.org/wiki/ID3_algorithm)和[C4.5](http://en.wikipedia.org/wiki/C4.5_algorithm)两种常用算法。
### ID3算法
从[信息论](http://en.wikipedia.org/wiki/Information_theory)知识中我们直到，期望信息越小，[信息增益](http://en.wikipedia.org/wiki/Information_gain)越大，从而纯度越高。所以ID3算法的核心思想就是以信息增益度量属性选择，选择分裂后信息增益最大的属性进行分裂。下面先定义几个要用到的概念。
设D为用类别对训练元组进行的划分，则D的[熵（entropy）](http://en.wikipedia.org/wiki/Entropy)表示为：
info(D)=−∑i=1mpilog2(pi)
其中pi表示第i个类别在整个训练元组中出现的概率，可以用属于此类别元素的数量除以训练元组元素总数量作为估计。熵的实际意义表示是D中元组的类标号所需要的平均信息量。
现在我们假设将训练元组D按属性A进行划分，则A对D划分的期望信息为：
infoA(D)=∑j=1v|Dj||D|info(Dj)
而信息增益即为两者的差值：
gain(A)=info(D)−info(D)
ID3算法就是在每次需要分裂时，计算每个属性的增益率，然后选择增益率最大的属性进行分裂。下面我们继续用SNS社区中不真实账号检测的例子说明如何使用ID3算法构造决策树。为了简单起见，我们假设训练集合包含10个元素：
我们也介绍一下这里面的专有名词：
- 日至密度、好友密度、是否使用真实的头像——条件属性
-
账号是否为真实——决策属性
其中s、m和l分别表示小、中和大。
设L、F、H和R表示日志密度、好友密度、是否使用真实头像和账号是否真实，下面计算各属性的信息增益。
info(D)=−0.7log20.7−0.3log20.3=0.879
infoL(D)=0.3∗(−03log203−33log233)+0.4∗(−14log214−34log234)+0.3∗(−13log213−23log223)=0.603
因此日志密度的信息增益是0.276。
用同样方法得到H和F的信息增益分别为0.033和0.553。
因为F具有最大的信息增益，所以第一次分裂选择F为分裂属性，分裂后的结果如下图表示：
在上图的基础上，再递归使用这个方法计算子节点的分裂属性，最终就可以得到整个决策树。
上面为了简便，将特征属性离散化了，其实日志密度和好友密度都是连续的属性。对于特征属性为连续值，可以如此使用ID3算法：
先将D中元素按照特征属性排序，则每两个相邻元素的中间点可以看做潜在分裂点，从第一个潜在分裂点开始，分裂D并计算两个集合的期望信息，具有最小期望信息的点称为这个属性的最佳分裂点，其信息期望作为此属性的信息期望。
ID3算法存在一些问题：
- (1)信息增益的计算依赖于特征数目较多的特征，而属性取值最多的属性并不一定最优。
- (2)ID3是非递增算法。
- (3)ID3是单变量决策树(在分枝节点上只考虑单个属性)，许多复杂概念的表达困难，属性相互关系强调不够，容易导致决策树中子树的重复或有些属性在决策树的某一路径上被检验多次。
- (4)抗噪性差，训练例子中正例和反例的比例较难控制。
由于ID3算法在实际应用中存在一些问题，于是Quilan提出了C4.5算法，严格上说C4.5只能是ID3的一个改进算法。
### C4.5算法
ID3算法存在一个问题，就是偏向于多值属性，例如，如果存在唯一标识属性ID，则ID3会选择它作为分裂属性，这样虽然使得划分充分纯净，但这种划分对分类几乎毫无用处。ID3的后继算法C4.5使用[增益率（gain ratio）](http://en.wikipedia.org/wiki/Information_gain_ratio)的信息增益扩充，试图克服这个偏倚。
C4.5算法首先定义了“分裂信息”，其定义可以表示成：
split_infoA(D)=−∑j=1v|Dj||D|log2(|Dj||D|)
其中各符号意义与ID3算法相同，然后，增益率被定义为：
gainratio(A)=gain(A)split_info(A)
C4.5选择具有最大增益率的属性作为分裂属性，其具体应用与ID3类似，不再赘述。
### 举例
### 代码
神经网络学习 之 BP神经网络
最后更新于：2022-04-01 16:07:03
上一次我们讲了[M-P模型](http://blog.csdn.net/u013007900/article/details/50066315)，它实际上就是对单个神经元的一种建模，还不足以模拟人脑神经系统的功能。由这些人工神经元构建出来的网络，才能够具有学习、联想、记忆和模式识别的能力。BP网络就是一种简单的人工神经网络。
本文具体来介绍一下一种非常常见的神经网络模型——反向传播(Back Propagation)神经网络。
### 概述
BP（Back Propagation）神经网络是1986年由Rumelhart和McCelland为首的科研小组提出，参见他们发表在Nature上的论文 [Learning representations by back-propagating errors](http://www.cs.toronto.edu/~hinton/absps/naturebp.pdf) 。
BP神经网络是一种按误差逆传播算法训练的多层前馈网络，是目前应用最广泛的神经网络模型之一。BP网络能学习和存贮大量的 输入-输出模式映射关系，而无需事前揭示描述这种映射关系的数学方程。它的学习规则是使用最速下降法，通过反向传播来不断 调整网络的权值和阈值，使网络的误差平方和最小。
### BP算法的基本思想
上一次我们说到，多层感知器在如何获取隐层的权值的问题上遇到了瓶颈。既然我们无法直接得到隐层的权值，能否先通过输出层得到输出结果和期望输出的误差来间接调整隐层的权值呢？BP算法就是采用这样的思想设计出来的算法，它的基本思想是,学习过程由信号的正向传播与误差的反向传播两个过程组成。
- 正向传播时，输入样本从输入层传入,经各隐层逐层处理后,传向输出层。若输出层的实际输出与期望的输出(教师信号)不符,则转入误差的反向传播阶段。
- 反向传播时，将输出以某种形式通过隐层向输入层逐层反传,并将误差分摊给各层的所有单元,从而获得各层单元的误差信号,此误差信号即作为修正各单元权值的依据。
这两个过程的具体流程会在后文介绍。
BP算法的信号流向图如下图所示
### BP网络特性分析——BP三要素
我们分析一个ANN时，通常都是从它的三要素入手，即
1)网络拓扑结构；
2)传递函数；
3)学习算法。
每一个要素的特性加起来就决定了这个ANN的功能特性。所以，我们也从这三要素入手对BP网络的研究。
### 3.1 BP网络的拓扑结构
上一次已经说了，BP网络实际上就是多层感知器，因此它的拓扑结构和多层感知器的拓扑结构相同。由于单隐层（三层）感知器已经能够解决简单的非线性问题，因此应用最为普遍。三层感知器的拓扑结构如下图所示。
一个最简单的三层BP：
### 3.2 BP网络的传递函数
BP网络采用的传递函数是非线性变换函数——[Sigmoid函数](http://blog.csdn.net/u013007900/article/details/50066315#t6)（又称S函数）。其特点是函数本身及其导数都是连续的，因而在处理上十分方便。为什么要选择这个函数，等下在介绍BP网络的学习算法的时候会进行进一步的介绍。
单极性S型函数曲线如下图所示。
f(x)=11+e−x
双极性S型函数曲线如下图所示。
f(x)=1−e−x1+e−x
### 3.3 BP网络的学习算法
BP网络的学习算法就是BP算法，又叫 δ 算法（在ANN的学习过程中我们会发现不少具有多个名称的术语）， 以三层感知器为例，当网络输出与期望输出不等时，存在输出误差 E ，定义如下
E=12(d−O)2=12∑ℓκ=1(dk−ok)2
将以上误差定义式展开至隐层，有
E=12∑ℓκ=1[dκ−f(netκ)]2=12∑ℓκ=1[dκ−f(∑mj=0ωjκyj)]2
进一步展开至输入层，有
E=12∑ℓκ=1dκ−f[∑mj=0ωjκf(netj)]2=12∑ℓκ=1dκ−f[∑mj=0ωjκf(∑nj=0υijχi)]2
由上式可以看出，网络输入误差是各层权值ωjκ、υij的函数，因此调整权值可改变误差 E。 显然，调整权值的原则是使误差不断减小，因此应使权值与误差的梯度下降成正比，即
Δωjκ=−η∂E∂ωjκj=0,1,2,…,m;κ=1,2,…,ℓ
Δυij=−η∂E∂υiji=0,1,2,…,n;j=1,2,…,m
对于一般多层感知器，设共有 h 个隐层，按前向顺序各隐层节点数分别记为 m1,m2,…,mh，各隐层输出分别记为 y1,y2,…,yh，各层权值矩阵分别记为 W1,W2,…,Wh,Wh+1，则各层权值调整公式为
输出层
Δωh+1jκ=ηδκh+1yhj=η(dκ−oκ)oκ(1−oκ)yκj(j=0,1,2,…,mh;κ=1,2,…,ℓ)
第 h 隐层
Δωhij=ηδhjyhi−1=η(∑lκ=1δoκωh+1jκyκj(1−ykjappa)yhi−1(i=0,1,2,…,m(h−1);j=1,2,…,mh)
按以上规律逐层类推，则第一隐层权值调整公式
Δω1pq=ηδ1qχp=η(∑m2r=1δ……2——rω2qr)y1q(1−y1q)χp(p=0,1,2,…,n;j=1,2,…,m1)
容易看出，BP学习算法中，各层权值调整公式形式上都是一样的，均由3个因素决定，即：
1. 学习率 η
1. 本层输出的误差信号δ
1. 本层输入信号 Y（或X）
其中输入层误差信号与网络的期望输出与实际输出之差有关，直接反应了输出误差，而各隐层的误差信号与前面各层的误差信号有关，是从输出层开始逐层反传过来的。
可以看出BP算法属于δ学习规则类，这类算法常被称为误差的梯度下降算法。δ学习规则可以看成是Widrow-Hoff(LMS)学习规则的一般化(generalize)情况。LMS学习规则与神经元采用的变换函数无关，因而不需要对变换函数求导，δ学习规则则没有这个性质，要求变换函数可导。这就是为什么我们前面采用Sigmoid函数的原因。
综上所述，BP三要素如下图所示。
下面我们会介绍BP网络的学习训练的具体过程。
### **BP网络的训练分解**
训练一个BP神经网络，实际上就是调整网络的权重和偏置这两个参数，BP神经网络的训练过程分两部分：
- 前向传输，逐层波浪式的传递输出值；
- 逆向反馈，反向逐层调整权重和偏置；
我们先来看前向传输。
### 前向传输（Feed-Forward前向反馈）
在训练网络之前，我们需要随机初始化权重和偏置，对每一个权重取[−1,1]的一个随机实数，每一个偏置取[0,1]的一个随机实数，之后就开始进行前向传输。
神经网络的训练是由多趟迭代完成的，每一趟迭代都使用训练集的所有记录，而每一次训练网络只使用一条记录，抽象的描述如下：
~~~
while 终止条件未满足：
for record:dataset:
trainModel(record)
~~~
首先设置输入层的输出值，假设属性的个数为100，那我们就设置输入层的神经单元个数为100，输入层的结点Ni为记录第i维上的属性值xi。对输入层的操作就这么简单，之后的每层就要复杂一些了，除输入层外，其他各层的输入值是上一层输入值按权重累加的结果值加上偏置，每个结点的输出值等该结点的输入值作变换
前向传输的输出层的计算过程公式如下：
Ij=∑i=1ωijoi+θj
oj=f(Ij)=11+eIj
对隐藏层和输出层的每一个结点都按照如上图的方式计算输出值，就完成前向传播的过程，紧接着是进行逆向反馈。
### 逆向反馈（Backpropagation）
逆向反馈从最后一层即输出层开始，我们训练神经网络作分类的目的往往是希望最后一层的输出能够描述数据记录的类别，比如对于一个二分类的问题，我们常常用两个神经单元作为输出层，如果输出层的第一个神经单元的输出值比第二个神经单元大，我们认为这个数据记录属于第一类，否则属于第二类。
还记得我们第一次前向反馈时，整个网络的权重和偏置都是我们随机取，因此网络的输出肯定还不能描述记录的类别，因此需要调整网络的参数，即权重值和偏置值，而调整的依据就是网络的输出层的输出值与类别之间的差异，通过调整参数来缩小这个差异，这就是神经网络的优化目标。对于输出层：
Ej=Oj(1−Oj)(Tj−Oj)
其中Ej表示第j个结点的误差值，Oj表示第j个结点的输出值，Tj记录输出值，比如对于2分类问题，我们用01表示类标1,10表示类别2，如果一个记录属于类别1，那么其T1=0，T2=1。
中间的隐藏层并不直接与数据记录的类别打交道，而是通过下一层的所有结点误差按权重累加，计算公式如下：
Ej=Oj(1−Oj)∑kEkWjk
其中Wjk表示当前层的结点j到下一层的结点k的权重值，Ek下一层的结点k的误差率。
计算完误差率后，就可以利用误差率对权重和偏置进行更新，首先看权重的更新：
其中λ表示表示学习速率，取值为0到1，学习速率设置得大，训练收敛更快，但容易陷入局部最优解，学习速率设置得比较小的话，收敛速度较慢，但能一步步逼近全局最优解。
更新完权重后，还有最后一项参数需要更新，即偏置：
至此，我们完成了一次神经网络的训练过程，通过不断的使用所有数据记录进行训练，从而得到一个分类模型。不断地迭代，不可能无休止的下去，总归有个终止条件
### 训练终止条件
每一轮训练都使用数据集的所有记录，但什么时候停止，停止条件有下面两种：
1. 设置最大迭代次数，比如使用数据集迭代100次后停止训练
1. 计算训练集在网络上的预测准确率，达到一定门限值后停止训练
### BP网络运行的具体流程
### 网络结构
输入层有n个神经元，隐含层有p个神经元，输出层有q个神经元。
### 变量定义
输入变量：
x=(x1,x2,…,xn)
隐含层输入变量：
hi=(hi1,hi2,…,hip)
隐含层输出变量：
ho=(ho1,ho2,…,hop)
输出层输入变量：
yi=(yi1,yi2,…,yiq)
输出层输出变量：
yo=(yo1,yo2,…,yoq)
期望输出向量：
do=(d1,d2,…,dq)
输入层与中间层的连接权值：
wih
隐含层与输出层的连接权值：
who
隐含层各神经元的阈值：
bh
输出层各神经元的阈值：
bo
样本数据个数：
k=1,2,…,m
激活函数：
f(⋅)
误差函数：
e=12∑o=1q(do(k)−yoo(k))2
### 第一步：网络初始化
给各连接权值分别赋一个区间(−1,1)内的随机数，设定误差函数e，给定计算精度值ε和最大学习次数M。
### 第二步：随机选取
随机选取第k个输入样本以及对应的期望输出
x(k)=(x1(k),x2(k),…,xn(k))do(k)=(d1(k),d2(k),…,dq(k))
### 第三部：隐含层计算
计算隐含层各神经元的输入和输出
hih(k)=∑i=1nwihxi(k)−bh(h=1,2,…,p)
hih(k)=f(hih(k))(h=1,2,…,p)
yio(k)=∑h=1pwhohoh(k)−bo(o=1,2,…,q)
yoo(k)=f(yio(k))(o=1,2,…,q)
### 第四步：求偏导数
利用网络期望输出和实际输出，计算误差函数对输出层的各神经元的偏导数δo(k)
### 第六步：修正权值
利用输出层各神经元的δo(k)和隐含层各神经元的输出来修正连接权值who(k)。
### 第七部：修正权值
利用隐含层各神经元的δh(k)和输入层各神经元的输入修正连接权值。
### 第八步：计算全局误差
E=12m∑k=1m∑o=1q(do(k)−yo(k))2
### 第九步：判断模型合理性
判断网络误差是否满足要求。
当误差达到预设精度或者学习次数大于设计的最大次数，则结束算法。
否则，选取下一个学习样本以及对应的输出期望，返回第三部，进入下一轮学习。
### BP网络的设计
在进行BP网络的设计是，一般应从网络的层数、每层中的神经元个数和激活函数、初始值以及学习速率等几个方面来进行考虑，下面是一些选取的原则。
### 1.网络的层数
理论已经证明，具有偏差和至少一个S型隐层加上一个线性输出层的网络，能够逼近任何有理函数，增加层数可以进一步降低误差，提高精度，但同时也是网络 复杂化。另外不能用仅具有非线性激活函数的单层网络来解决问题，因为能用单层网络解决的问题，用自适应线性网络也一定能解决，而且自适应线性网络的 运算速度更快，而对于只能用非线性函数解决的问题，单层精度又不够高，也只有增加层数才能达到期望的结果。
### 2.隐层神经元的个数
网络训练精度的提高，可以通过采用一个隐含层，而增加其神经元个数的方法来获得，这在结构实现上要比增加网络层数简单得多。一般而言，我们用精度和 训练网络的时间来恒量一个神经网络设计的好坏：
（1）神经元数太少时，网络不能很好的学习，训练迭代的次数也比较多，训练精度也不高。
（2）神经元数太多时，网络的功能越强大，精确度也更高，训练迭代的次数也大，可能会出现过拟合(over fitting)现象。
由此，我们得到神经网络隐层神经元个数的选取原则是：在能够解决问题的前提下，再加上一两个神经元，以加快误差下降速度即可。
### 3.初始权值的选取
一般初始权值是取值在(−1,1)之间的随机数。另外威得罗等人在分析了两层网络是如何对一个函数进行训练后，提出选择初始权值量级为s√r的策略， 其中r为输入个数，s为第一层神经元个数。
### 4.学习速率
学习速率一般选取为0.01−0.8，大的学习速率可能导致系统的不稳定，但小的学习速率导致收敛太慢，需要较长的训练时间。对于较复杂的网络， 在误差曲面的不同位置可能需要不同的学习速率，为了减少寻找学习速率的训练次数及时间，比较合适的方法是采用变化的自适应学习速率，使网络在 不同的阶段设置不同大小的学习速率。
### 5.期望误差的选取
在设计网络的过程中，期望误差值也应当通过对比训练后确定一个合适的值，这个合适的值是相对于所需要的隐层节点数来确定的。一般情况下，可以同时对两个不同 的期望误差值的网络进行训练，最后通过综合因素来确定其中一个网络。
### BP网络的局限性
BP网络具有以下的几个问题：
- (1)需要较长的训练时间：这主要是由于学习速率太小所造成的，可采用变化的或自适应的学习速率来加以改进。
- (2)完全不能训练：这主要表现在网络的麻痹上，通常为了避免这种情况的产生，一是选取较小的初始权值，而是采用较小的学习速率。
- (3)局部最小值：这里采用的梯度下降法可能收敛到局部最小值，采用多层网络或较多的神经元，有可能得到更好的结果。
### BP网络的改进
P算法改进的主要目标是加快训练速度，避免陷入局部极小值等，常见的改进方法有带动量因子算法、自适应学习速率、变化的学习速率以及作用函数后缩法等。 动量因子法的基本思想是在反向传播的基础上，在每一个权值的变化上加上一项正比于前次权值变化的值，并根据反向传播法来产生新的权值变化。而自适应学习 速率的方法则是针对一些特定的问题的。改变学习速率的方法的原则是，若连续几次迭代中，若目标函数对某个权倒数的符号相同，则这个权的学习速率增加， 反之若符号相反则减小它的学习速率。而作用函数后缩法则是将作用函数进行平移，即加上一个常数。
### BP网络实现
由于BP网络具有出色的非线性映射能力、泛化能力和容错能力，因此BP网络成了至今为止应用最广泛的人工神经网络。下图是Matlab下用BP网络做线性拟合的结果，效果很好。
~~~
% BP网络函数逼近实例
% 1.首先定义正弦函数，采样率为20Hz，频率为1Hz
k = 1; % 设定正弦信号频率
p = [0:0.05:4];
t = cos(k*pi*p) + 3*sin(pi*p);
plot(p, t, '-'), xlabel('时间'); ylabel('输入信号');
% 2.生成BP网络。用newff函数生成前向型BP网络，设定隐层中神经元数目为10
% 分别选择隐层的传递函数为 tansig，输出层的传递函数为 purelin，
% 学习算法为trainlm。
net =
newff(minmax(p),[10,10,1],{'tansig','tansig','purelin'},'trainlm');
% 3.对生成的网络进行仿真并做图显示。
y1 = sim(net,p); plot(p, t, '-', p, y1, '--')
% 4.训练。对网络进行训练，设定训练误差目标为 1e-5，最大迭代次数为300，
% 学习速率为0.05。
net.trainParam.lr=0.05;
net.trainParam.epochs=1000;
net.trainParam.goal=1e-5;
[net,tr]=train(net,p,t);
%5.再次对生成的网络进行仿真并做图显示。
y2 = sim(net,p);
plot(p, t, '-', p, y2, '--')
~~~
这是用C语言写的：用BP神经网络拟合函数：Y=sin(X)
~~~
#include "math.h"
#include "time.h"
#include "stdio.h"
#include "stdlib.h"
#include "ctype.h"
#define Ni 1
#define Nm 4
#define No 1
#define L 100
#define Enom 0.02
#define loopmax 100000
#define e 2.71828
double E;
double a,u,n;
double W1[Ni][Nm],D1[Ni][Nm],W2[Nm][No],D2[Nm][No];
double D22[Nm][No],D11[Ni][No];
double a1[Ni][Nm],a2[Nm][No];
double Pi[L][Ni],Pm[L][Nm],Po[L][No],T[L][No];
double Xm[L][Nm],Xo[L][No];
double Qm[L][Nm],Qo[L][No];
void proceed();
void proceedR();
void forQ();
void amend();
void initiate();
double newa(double a,double D);
double cal(double d);
double vcal(double d);
main()
{
long int i;
int flag;
char choice;
for(;;)
{
flag=0;
initiate();
for(i=0;;i++)
{
proceed();
if( E < Enom )
{
flag=1;
break;
}
if( i >= loopmax)
{
flag = -1;
break;
}
if(i%2500==0)
printf("第%10d轮误差：%20f,学习速率:%10f\n",i,E,a1[0][0]);
forQ();
amend();
}
if(flag>0)proceedR();
else printf("训练失败！\n");
for(;;)
{
choice=getchar();
printf("是否继续？(Y/N)\n");
choice=getchar();
choice=toupper(choice);
if(choice=='Y')break;
if(choice=='N')exit(0);
}
}
}
void initiate()
{
int i,j;
int random;
double x;
double step;
int stime;
long ltime;
ltime=time(NULL);
stime=(unsigned)ltime/2;
srand(stime);
a=0.02;
u=1;
n=1;
printf("本程序将用BP神经网络拟合函数：Y=sin(X)\n\n");
for( i=0; i<Nm; i++)
{
for( j=0; j<Ni; j++)
{
random=rand()%100-50;
x=random;
x=x/100;
W1[j][i]=x;
D11[j][i]=0;
D1[j][i]=0;
a1[j][i]=0.01;
}
for( j=0; j<No; j++)
{
random=rand()%100-50;
x=random;
x=x/100;
W2[i][j]=x;
D22[i][j]=0;
D2[i][j]=0;
a2[i][j]=0.01;
}
}
step=1.0/L;
for(i=0;i<L;i++)
{
x=i;
Pi[i][0]=x*step;
T[i][0]=sin(Pi[i][0]);
}
printf("初始化成功!\n\n下面将对神经网络进行训练请稍候。\n");
}
void proceed()
{
int i, j, k;
E=0 ;
for( i=0; i<L; i++ )
{
for( j=0; j<Nm; j++ )
{
Pm[i][j] = 0;
for( k=0; k<Ni; k++ )
{
Pm[i][j] = Pi[i][k] * W1[k][j] + Pm[i][j];
}
Xm[i][j] = cal( Pm[i][j] );
}
for( j=0; j<No; j++)
{
Po[i][j] = 0;
for( k=0; k<Nm; k++)
{
Po[i][j] = Xm[i][k] * W2[k][j] + Po[i][j];
}
Xo[i][j] = cal( Po[i][j] );
E = E + ( Xo[i][j] - T[i][j] ) * ( Xo[i][j] - T[i][j] ) / 2;
}
}
}
void forQ()
{
int i,j,k;
for( i=0; i<L; i++ )
{
for( j=0; j<No; j++)
{
Qo[i][j] = ( T[i][j] - Xo[i][j] )* vcal( Xo[i][j] );
}
for(j=0; j<Nm; j++)
{
Qm[i][j]=0;
for( k=0; k<No; k++)
{
Qm[i][j] = Qo[i][k] * W2[j][k] + Qm[i][j];
}
Qm[i][j] = Qm[i][j] * vcal( Xm[i][j] );
}
}
}
void amend()
{
int i,j,k;
double D;
for( i=0; i<Nm; i++)
{
for( j=0; j<Ni; j++)
{
D1[j][i]=0;
}
for( j=0; j<No; j++)
{
D2[i][j]=0;
}
}
for( i=0; i<Ni; i++)
{
for( j=0; j<Nm; j++)
{
for( k=0; k<L; k++)
{
D1[i][j] = Qm[k][j] * Pi[k][i] + D1[i][j];
}
D = D1[i][j] * D11[i][j] ;//为D11付初值
a1[i][j] = newa( a1[i][j] , D ); // a 付初值
W1[i][j] = W1[i][j] + a1[i][j] * ( n * D1[i][j] + ( 1 - n ) * D11[i][j] );
D11[i][j] = D1[i][j];
}
}
for( i=0; i<Nm; i++)
{
for( j=0; j<No; j++)
{
for( k=0; k<L; k++)
{
D2[i][j] = Qo[k][j] * Xm[k][i] + D2[i][j];
}
D = D2[i][j] * D22[i][j] ;//为D11付初值
a2[i][j] = newa( a2[i][j] , D );
W2[i][j] = W2[i][j] + a2[i][j] * ( n * D2[i][j] + ( 1 - n ) * D22[i][j] );
D22[i][j] = D2[i][j];
}
}
}
void proceedR()
{
int i, j;
float x;
double input,output;
char choice;
for(;;)
{
for(;;)
{
printf("在此输入需要计算的值(0,1):\n");
scanf("%f",&x);
input=(double)x;
if((input>=0)&(input<=1))break;
printf("注意输入值应介于0、1之间！\n");
for(;;)
{
choice=getchar();
printf("是否继续？(Y/N)\n");
choice=getchar();
choice=toupper(choice);
if(choice=='Y')break;
if(choice=='N')exit(0);
}
}
for(i=0;i<Nm;i++)
{
Pm[0][i]=0;
for( j=0; j<Ni; j++ )
{
Pm[0][i] = input* W1[j][i]+Pm[0][i] ;
}
Xm[0][i] = cal( Pm[0][i] );
}
for( i=0; i<No; i++)
{
Po[0][i] = 0;
for( j=0; j<Nm; j++)
{
Po[0][i] = Xm[0][j] * W2[j][i]+Po[0][i];
}
}
output=cal( Po[0][0] );
printf("输入值为%20f对应的结果为%f\n",input,output);
printf("输入值为%20f对应的正常结果为%f\n",input,sin(input));
for(;;)
{
choice=getchar();
printf("是否继续？(Y/N)\n");
choice=getchar();
choice=toupper(choice);
if(choice=='Y')break;
if(choice=='N')exit(0);
}
}
}
double newa(double a, double D)
{
if( D > 0 )
{
{
if(a<=0.04)
a = a * 2;
else a=0.08;
}
}
else
if ( D < 0)
{
if(a>=0.02)
{
a = a / 2;
}
else a=0.01;
}
return a;
}
double cal(double d)
{
d = - (d * u); // chushihua
d = exp( d );
d = 1 / ( 1 + d );
return d;
}
double vcal(double d)
{
return u * d * ( 1 - d );
}
~~~
神经网络学习 之 M-P模型
最后更新于：2022-04-01 16:07:01
### M-P模型的来源
所谓M-P模型，其实是按照生物神经元的结构和工作原理构造出来的一个抽象和简化了的模型。
下图是生物神经元结构。
大家可以查一查一些生物方面的书籍，了解一下这个神经元是如何工作的。我们可以概括出生物神经网络的假定特点：
1. 每个神经元都是一个**多输入单输出**的信息处理单元；
2. 神经元输入**分兴奋性输入**和**抑制性输入**两种类型；
3. 神经元具有**空间整合特性**和**阈值特性**；
4. 神经元输入与输出间有固定的**时滞**，主要取决于突触延搁
### M-P模型
按照生物神经元，我们建立M-P模型。为了使得建模更加简单，以便于进行形式化表达，我们忽略时间整合作用、不应期等复杂因素，并把神经元的突触时延和强度当成常数。下图就是一个M-P模型的示意图。
那么接下来就好类比理解了。我们将这个模型和生物神经元的特性列表来比较：
结合M-P模型示意图来看，对于某一个神经元j （注意别混淆成变量了，在这里j 只是起到标识某个神经元的作用），它可能接受同时接受了许多个输入信号，用 χi 表示。
由于生物神经元具有不同的突触性质和突触强度，所以对神经元的影响不同，我们用权值 ωij 来表示，其正负模拟了生物神经元中突出的兴奋和抑制，其大小则代表了突出的不同连接强度。
θj表示为一个阈值（threshold），或称为偏置（bias）。
由于累加性，我们对全部输入信号进行累加整合，相当于生物神经元中的膜电位（水的变化总量），其值就为：
net'j(t)=∑i=1nωijχi(t)−θj
神经元激活与否（外接专用水管流出与否）取决于某一阈值电平（水位高度），即只有当其输入总和超过阈值θj 时，神经元才被激活而发放脉冲，否则神经元不会发生输出信号。整个过程可以用下面这个函数来表示：
yj=f(netj)
yj表示神经元j的输出，函数f称为激活函数 ( Activation Function )或转移函数 ( Transfer Function ) ，net'j(t)称为净激活(net activation)。
若将阈值看成是神经元j的一个输入x0的权重w0j，则上面的式子可以简化为：
net'j(t)=∑i=0nωijχi(t)
yj=f(netj)
若用X表示输入向量，用W表示权重向量，即：
X=[χ0,χ1,…,χn]
M=⎡⎣⎢⎢⎢⎢⎢ω0jω1j⋮ωnj⎤⎦⎥⎥⎥⎥⎥
则神经元的输出可以表示为向量相乘的形式：
netj=XW
yj=f(netj)=f(XW)
若神经元的净激活net为正，称该神经元处于激活状态或兴奋状态(fire)，若净激活net为负，则称神经元处于抑制状态。
由此我们可以得到总结出M-P模型的6个特点：
1. 每个神经元都是一个多输入单输出的信息处理单元；
2. 神经元输入分兴奋性输入和抑制性输入两种类型；
3. 神经元具有空间整合特性和阈值特性；
4. 神经元输入与输出间有固定的时滞，主要取决于突触延搁；
5. 忽略时间整合作用和不应期；
6. 神经元本身是非时变的，即其突触时延和突触强度均为常数。
前面4点和生物神经元保持一致。
结合公式来看，输入χij的下标 i=1,2,...,n，输出 yj 的下标j 体现了第1个特点“多输入单输出”；
权重值 ωij 的正负体现了第2个特点中“突触的兴奋与抑制”；
θj 代表第3个特点 中的阈值，当 net'j(t)−Tj>0 时，神经元才能被激活；
为了简单起见，对膜电位的计算net'j(t) 并没有考虑时间整合，只考虑了空间整合，即只对每条神经末梢传来的信号根据权重进行累加整合，而没有考虑输入输出间的突触时延,体现了第5个特点。
这种“阈值加权和”的神经元模型称为**M-P模型 ( McCulloch-Pitts Model )**，也称为神经网络的一个**处理单元( PE, Processing Element )**。
### 常用激活函数
激活函数的选择是构建神经网络过程中的重要环节，下面简要介绍常用的激活函数。
### (1) 线性函数 ( Liner Function )
### (2) 斜面函数 ( Ramp Function )
### (3) 阈值函数 ( Threshold Function )
以上3个激活函数都属于线性函数，下面介绍两个常用的非线性激活函数。
### (4) S形函数 ( Sigmoid Function )
该函数的导函数：
### (5) 双极S形函数
该函数的导函数：
S形函数与双极S形函数的图像如下：
双极S形函数与S形函数主要区别在于函数的值域，双极S形函数值域是(−1,1)，而S形函数值域是(0,1)。
由于S形函数与双极S形函数都是可导的(导函数是连续函数)，因此适合用在BP神经网络中。（BP算法要求激活函数可导）
### 最简单的神经网络结构——感知器
在1958年,美国心理学家Frank Rosenblatt提出一种具有单层计算单元的神经网络,称为感知器(Perceptron)。它其实就是基于M-P模型的结构。我们可以看看它的拓扑结构图。
这个结构非常简单，如果你还记得前面所讲的M-P神经元的结构的话，这个图其实就是输入输出两层神经元之间的简单连接（如果忘了可以看看第一话的模型示意图）。
由第一话的(2)中我们知道输入层各节点的输入加权和
net'j=∑i=1nωijxi
我们一般采用符号函数来当作单层感知器的传递函数，即输出
oj=sgn(net′j−θj)=sgn(∑i=0nωijxi)=sgn(WTjX)
公式(2)可以进一步表达为：
oj={ 1,−1,WTjX>0WTjX<0
### 单层感知器的局限性
虽然单层感知器简单而优雅，但它显然不够聪明——它仅对线性问题具有分类能力。什么是线性问题呢？简单来讲，就是用一条直线可分的图形。比如，逻辑“与”和逻辑“或”就是线性问题，我们可以用一条直线来分隔0和1。
1）逻辑“与”的真值表和二维样本图如图2：
2）逻辑“或”的真值表如图3：
为什么感知器就可以解决线性问题呢？这是由它的传递函数决定的。这里以两个输入分量 x1 和 x2 组成的二维空间为例，此时节点 j 的输出为
oj={ 1,−1,ω1jx1+ω2jx2−θj>0ω1jx1+ω2jx2−θj<0
所以，方程
ω1jx1+ω2jx2−θj=0
确定的直线就是二维输入样本空间上的一条分界线。对于三维及更高维数的推导过程可以参考其他的Tutorials。
如果要让它来处理非线性的问题，单层感知器网就无能为力了。例如下面的“异或”，就无法用一条直线来分割开来，因此单层感知器网就没办法实现“异或”的功能。
### 多层感知器的瓶颈
既然一条直线无法解决分类问题，当然就会有人想到用弯曲的折线来进行样本分类。我们常常听到一句批评人笨的话“你这人脑袋就是不会转弯！”大意就是如此，脑袋会转弯的人才善于解决问题。所以，人们请来了单层感知器他哥——多层感知器来帮忙。所谓多层感知器，就是在输入层和输出层之间加入隐层，，以形成能够将样本正确分类的凸域。多层感知器的拓扑结构如下图所示。
我们可以比较一下单层感知器和多层感知器的分类能力：
由上图可以看出，随着隐层层数的增多，凸域将可以形成任意的形状，因此可以解决任何复杂的分类问题。实际上，Kolmogorov理论指出：**双隐层感知器就足以解决任何复杂的分类问题**。
多层感知器确实是非常理想的分类器，但问题也随之而来：隐层的权值怎么训练？对于各隐层的节点来说，它们并不存在期望输出，所以也无法通过感知器的学习规则来训练多层感知器。因此，多层感知器心有余而力不足，虽然武功高强，但却无力可施。
### ANN的低潮期
1966年，Minisky和Papert在他们的《感知器》一书中提出了上述的感知器的研究瓶颈，指出理论上还不能证明将感知器模型扩展到多层网络是有意义的。这在人工神经网络的历史上书写了极其灰暗的一章。对ANN的研究，始于1890年开始于美国著名心理学家W.James对于人脑结构与功能的研究，半个世纪后W.S.McCulloch和W.A.Pitts提出了M-P模型，之后的1958年Frank Rosenblatt在这个基础上又提出了感知器，此时对ANN的研究正处在升温阶段，《感知器》这本书的出现就刚好为这刚刚燃起的人工神经网络之火泼了一大盆冷水。一时间人们仿佛感觉对以感知器为基础的ANN的研究突然间走到尽头，看不到出路了。于是，几乎所有为ANN提供的研究基金都枯竭了，很多领域的专家纷纷放弃了这方面课题的研究。
### ANN研究的复苏和BP神经网络的诞生
所以说真理的果实总是垂青于能够忍受寂寞的科学家。尽管ANN的研究陷入了前所未有的低谷， 但仍有为数不多的学者忍受住寂寞，坚持致力于ANN的研究。在长达10年的低潮时期之间，相 继有一些开创性的研究成果被提出来，但还不足以激起人们对于ANN研究的热情。一直到上世 纪80年代，两个璀璨的成果诞生了：1982年美国加州理工学院的物理学家John J.Hopfield博 士的Hopfield网络和David E.Rumelhart以及James L.McCelland研究小组发表的《并行分布 式处理》。这两个成果重新激起了人们对ANN的研究兴趣，使人们对模仿脑信息处理的智能计 算机的研究重新充满了希望。
前者暂不讨论，后者对具有非线性连续变换函数的多层感知器的误差反向传播(Error Back Propagation)算法进行了详尽的分析，实现了 Minsky 关于多层网络的设想。Error Back Propagation算法的简称就是BP算法，以BP算法实现的多层感知器网络就是[BP网络](http://blog.csdn.net/u013007900/article/details/50118945)。
所以，BP网络本质上并不是一个新的网络，而是使用BP学习算法的多层感知器网络。
数值积分方法
最后更新于：2022-04-01 16:06:59
数值积分是工程师和科学家经常使用的基本工具，用来计算无法解析求解的定积分的近似解。
如：Φ(x)=∫x0t3et−1dt不存在Φ(x)的解析解，要求Φ(5)。
那么我们就要通过数值积分的方法来计算，数值积分的目的是，通过在有限个采样点上计算f(x)的值来逼近 f(x)在区间[a,b]上的定积分.
设a=x0<x1<…<xM=b. 称形如
且具有性质∫baf(x)dx=Q[f]+E[f]的公式为数值积分或面积公式。项E[f]称为积分的截断误差，值{xk}Mk=0称为面积节点， {wk}Mk=0称为权。
下面介绍几种常用的数值积分方法。
### 基于多项式差值的面积公式
通过M+1个等距点{(xk,f(xk))}Mk=0存在唯一的次数小于等于M的多项式PM(x)。当用该多项式来近似[a,b]上的f(x)时，PM(x)的积分就近似等于f(x)的积分，这类公式称为牛顿－科特斯公式。当使用采样点x0=a和xM=b时，称为闭型牛顿－科特斯公式。
### 左/中/右矩形公式、梯形公式
左/中/右矩形公式
梯形公式
图形如下
### 辛普森公式
#### 推导过程
若f(x)在[a,b]上有定义，将区间[a,b]分割成n等分（取n为偶数），则有a=x0<x1<…<xM=b，其中
xi=a+iΔx（∀i=0,1,…,n,Δx=b−an为步长）
这里我们想通过(x0,f(x0)),(x1,f(x1)),(x2,f(x2))三点抛物线g(x)=αx2+βx+γ来取代f(x)在[x0,x2]的定义，今儿求出它的近似积分值(如图)，最后用累加的方式求得f(x)在[a,b]上的近似积分。
由假设我们有
所以可得
∫baf(x)dx=∫x2x0+∫x4x2+⋯+∫xnxn−2f(x)dx
#### 误差估计
若令Sn=Δx3[f(x0)+4f(x1)+2f(x2)+4f(x3)⋯+2f(xn−2)+4f(xn−1)+f(xn)]且f(4)(x)在[a,b]上连续，则我们可以估计出辛普森公式的误差值为
#### 例题1
试用辛普森公式估计∫10e−x2dx，取n=6
解：
令f(x)=e−x2,Δx=16则
## 拉格朗日插值
### 分段线性插值
例如：函数f(x)=11+x2如果在区间[−5，5]上取11个等距节点:xk=−5+k(k=0,1,2,...,10),由Lagrange插值公式可得到f(x)的10次L10(x)。如图所示：
L10(x)仅在区间的中部能较好的逼近函数f(x),在其它部位差异较大，而且越接近端点，逼近效果越差。
可以证明，当节点无限加密时，Ln(x)也只能在很小的范围内收敛，这一现象称为Runge现象。它表明通过增加节点来提高逼近程度是不适宜的，因而不采用高次多项式插值。
#### 推导过程
已知函数f(x)在区间[xk,xk+1]的端点上的函数值yk=f(xk),yk+1=f(xk+1),求一个一次函数y=P1(x)使得yk=f(xk),yk+1=f(xk+1), 其几何意义是已知平面上两点(xk,yk),(xk+1,yk+1),求一条直线过该已知两点。
由直线的点斜式公式可知：
P1(x)=yk+yk+1−ykxk+1−xk(x−xk)
把此式按照yk和yk+1写成两项:
P1(x)=x−xk+1xk−xk+1yk+x−xkxk+1−xkyk+1
记
lk=x−xk+1xk−xk+1 lk+1=x−xkxk+1−xk
并称它们为一次插值基函数。该基函数的特点如下：
li(xk)={10k=ik!=i
对于i=0，
l0(x)=x−x1x0−x1,x∈[x0,x1]
其它点上，l0(x)=0;
对于i=1,2,…,n−1，
li(x)={x−xi−1xi−xi−1x−xi+1xi−xi+1x∈[xi−1,xi]x∈[xi,xi+1]
其它点上，li(x)=0; 对于i=n，
ln(x)=x−xn−1xn−xn−1,x∈[xn−1,xn]
其它点上，ln(x)=0. 于是，
P1(x)=∑k=0nyklk(x)
此表达式与前面的表达式是相同的,这是因为在区间[xk,xk+1]上, 只有lk(x),lk+1(x)是非零的，其它基函数均为零。
从而
P1(x)=yklk(x)+yk+1lk+1(x)
此形式称之为拉格朗日型插值多项式。其中, 插值基函数与yk、yk+1 无关，而由插值结点xk、xk+1 所决定。一次插值多项式是插值基函数的线性组合, 相应的组合系数是该点的函数值yk、yk+1 .
#### 误差估计
根据拉格朗日一次插值函数的余项,可以得到分段线性插值函数的插值误差估计: 对x∈[a,b],当x∈[xk,xk+1]时,
R(x)=12f′′(ξ)(x−xk)(x−xk+1)
则R(x)≤h28M，其中
h=max0≤k≤n−1|xk+1−xk|,M=maxx∈(a,b)f′′(x)
于是有下面的定理：
如果f(x)在[a,b]上二阶连续可微，则分段连续函数φ(x)的余项有以下误差估计：
|R(x)|=|f(x)−ϕ(x)|≤h28M
其中
h=max0≤k≤n−1|xk+1−xk|,M=maxx∈(a,b)f′′(x)
于是可以加密插值结点, 缩小插值区间, 使h减小, 从而减小插值误差。
#### 例题1
已知函数y=f(x)=11+x2在区间[0,5]上取等距插值节点(如下表)，
求区间上分段线性插值函数,并利用它求出f(4.5)近似值。
解:
在每个分段区间[k,k+1]上，
P1(x)=x−(k+1)k−(k+1)yk+x−(k)k+1−(k)yk+1=−yk(x−k−1)+yk+1(x−k)
于是
P1(4.5)=−0.05882×(4.5−5)+0.03846×(4.5−4)=0.04864
### 拉格朗日型二次插值多项式
#### 推导过程
已知函数y=f(x)在点xk−1,xk,xk+1上的函数值yk−1=f(xk−1),yk=f(xk),yk+1=f(xk+1), 求一个次数不超过二次的多项式P2(x), 使其满足,
P2(xk−1)=yk−1,P2(xk)=yk,P2(xk+1)=yk+1
其几何意义为:
已知平面上的三个点
(xk−1,yk−1),(xk,yk),(xk+1,yk+1),
求一个二次抛物线, 使得该抛物线经过这三点。
#### 插值基本多项式
有三个插值结点xk−1,xk,xk+1构造三个插值基本多项式，要求满足：
(1) 基本多项式为二次多项式；
(2) 它们的函数值满足下表：
因为lk−1(xk)=0,lk−1(xk+1)=0, 故有因子(x−xk)(x−xk+1), 而其已经是一个二次多项式, 仅相差一个常数倍, 可设
lk−1(x)=a(x−xk)(x−xk+1)
又因为
lk−1(xk−1)=1⇒a(xk−1−xk)(xk−1−xk+1)=1
得
a=1(xk−1−xk)(xk−1−xk+1)
从而
lk−1(x)=(x−xk)(x−xk+1)(xk−1−xk)(xk−1−xk+1)
同理得
lk(x)=(x−xk−1)(x−xk+1)(xk−xk−1)(xk−xk+1)
lk+1(x)=(x−xk−1)(x−xk)(xk+1−xk−1)(xk+1−xk)
#### 拉格朗日型二次插值多项式
由前述, 拉格朗日型二次插值多项式:
P2(x)=yk−1lk−1(x)+yklk(x)+yk+1lk+1(x)
P2(x)是三个二次插值多项式的线性组合,因而其是次数不超过二次的多项式,且满足:
P2(xi)=yi,(i=k−1,k,k+1)
#### 例题2
已知：
| xi | 10 | 15 | 20 |
|-----|-----|-----|-----|
| yi=lgxi | 1 | 1.1761 | 1.3010 |
利用此三值的二次插值多项式求lg(12)的近似值。
解：设x0=10,x1=15,x2=20,则:
故:
所以
#### 误差估计
我们在[a,b]上用多项式Pn(x) 来近似代替函数f(x), 其截断误差记作Rn(x)=f(x)−Pn(x)
当x在插值结点xi上时Rn(xi)=f(xi)−Pn(xi)=0,下面来估计截断误差：
定理1：
设函数y=f(x)的n阶导数y(n)=f(n)(x)在[a,b]上连续，
y(n+1)=f(n+1)(x)在(a,b)上存在；
插值结点为:a≤x0<x1<…<xn≤b,
Pn(x)是n次拉格朗日插值多项式；
则对任意x∈[a,b]有：
Rn(x)=1(n+1)!f(n+1)(ξ)ωn+1(x)
其中ξ∈(a,b), ξ依赖于x：ωn+1(x)=(x−x0)(x−x1)…(x−xn)
设maxa≤x≤b|fn+1(x)|≤Mn+1，则
|Rn(x)|≤1(n+1)!Mn+1ωn+1(x)
易知，线性插值的截断误差为:
R1(x)=12f′′(ξ)(x−x0)(x−x1)
二次插值的截断误差为:
R2(x)=16f(3)(ξ)(x−x0)(x−x1)(x−x2)
#### 例题3
在例2中，用lg10,lg15和lg20计算lg12.
P2(12)=1.0766,
e=|1.0792−1.0766|=0.0026
估计误差：
## 三次样条插值法
### 三次样条曲线原理
x:a=x0<x1<…<xM=by:y0 y1 … yn
#### 定义
样条曲线S(x)是一个分段定义的公式。给定n+1个数据点，共有n个区间，三次样条方程满足以下条件：
a. 在每个分段区间[xi,xi+1]（i=0,1,…,n−1，x递增），S(x)=Si(x) 都是一个三次多项式。
b. 满足S(xi)=yi（i=0,1,…,n）
c.S(x)，导数S′(x)，二阶导数S′′(x)在[a,b]区间都是连续的，即S(x)曲线是光滑的
所以n个三次多项式分段可以写作：
Si(x)=ai+bi(x−xi)+ci(x−xi)2+di(x−xi)3,i=0,1,…,n−1
其中ai,bi,ci,di代表4n个未知系数。
#### 求解
已知:
a. n+1个数据点[xi,yi],i=0,1,…,n
b. 每一分段都是三次多项式函数曲线进行逼近
c. 节点达到二阶连续
d. 左右两端点处特性（自然边界，固定边界，非节点边界）
根据定点，求出每段样条曲线方程中的系数，即可得到每段曲线的具体表达式。
插值和连续性：
, 其中 i=0,1,…,n−1
微分连续性：
, 其中i=0,1,…,n−2
样条曲线的微分式：
将步长hi=xi+1−xi带入样条曲线的条件：
a. 由Si(xi)=yi(i=0,1,…,n−1)推出
ai=yi
b. 由Si(xi+1)=yi+1(i=0,1,…,n−1)推出
ai+hibi+h2ici+h3idi=yi+1
c. 由S′i(xi+1)=S′i+1(xi+1)(i=0,1,…,n−2)推出
S′i(xi+1)=bi+2ci(xi+1−xi)+3di(xi+1−xi)2=bi+2cih+3dih2S′i+1(xi+1)=bi+1+2ci(xi+1−xi+1)+3di(xi+1−xi+1)2=bi+1
由此可得：
bi+2hici+3h2idi−bi+1=0
d. 由S′′i(xi+1)=S′′i+1(xi+1)(i=0,1,…,n−2)推出
2ci+6hidi−2ci+1=0
设mi=S′′i(xi)=2ci，则
a. 2ci+6hidi−ci+1=0可写为：
mi+6hidi−mi+1=0，推出
di=mi+1−mi6hi
b. 将ci,di带入 yi+hibi+h2ici+h3idi=yi+1可得：
bi=yi+1−yihi−hi2mi−hi6(mi+1−mi)
c. 将bi,ci,di带入bi+2hici+3hidi=bi+1(i=0,1,…,n−2)可得：
himi+2(hi+hi+1)mi+1+hi+1mi+2=6[yi+2−yi+1hi+1−yi+1−yihi]
#### 端点条件
由i的取值范围可知，共有n-1个公式， 但却有n+1个未知量m 。要想求解该方程组，还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。 选择不是唯一的，3种比较常用的限制如下。
##### a. 自由边界(Natural)
首尾两端没有受到任何让它们弯曲的力，即S′′=0。具体表示为m0=0和 mn=0
则要求解的方程组可写为：
##### b. 固定边界(Clamped)
首尾两端点的微分值是被指定的，这里分别定为A和B。则可以推出
将上述两个公式带入方程组，新的方程组左侧为
##### c. 非节点边界(Not-A-Knot)
指定样条曲线的三次微分匹配，即
S′′′o(x1)=S′′′1(x1)S′′′n−2(xn−1)=S′′′n−1(xn−1)
根据S′′′i(x)=6di和di=mi+1−mi6hi，则上述条件变为
h1(m1−m0)=h0(m2−m1)hn−1(mn−1−mn−2)=hn−2(mn−mn−1)
新的方程组系数矩阵可写为：
右下图可以看出不同的端点边界对样条曲线的影响：
#### 算法总结
假定有n+1个数据节点
(x0,y0),(x1,y1),(x2,y2),…,(xn,yn),
a. 计算步长hi=xi+1−xi(i=0,1,…,n−1)
b. 将数据节点和指定的首位端点条件带入矩阵方程
c. 解矩阵方程，求得二次微分值mi。
d. 计算样条曲线的系数：
其中i=0,1,…,n−1
e. 在每个子区间xi≤x≤xi+1中，创建方程
gi(x)=ai+bi(x−xi)+ci(x−xi)2+di(x−xi)3
单纯形法 — 求解线性规划
最后更新于：2022-04-01 16:06:56
目前，运用最广的线性规划方法就是著名的单纯形方法。这种方法是G.B.Dantzig在1947年提出的。几十年的实践证明，单纯形方法的确是一种使用方便、行之有效的重要算法。如今，它已经成为线性规划的中心内容。
单纯形法的基本思路是有选择地取（而不是枚举所有的）基本可行解，即是从可行域的一个顶点出发，沿着可行域的边界移到另一个相邻的顶点，要求新顶点的目标函数值不比原目标函数值差，如此迭代，直至找到最优解，或判定问题无界。
### 线性规划问题及其表示
线性规划问题可表示为如下形式：
变量满足约束条件(8.2)-(8.5)式的一组值称为线性规划问题的一个可行解。
所有可行解构成的集合称为线性规划问题的**可行区域**。
使目标函数取得极值的可行解称为**最优解**。
在最优解处目标函数的值称为**最优值**。
有些情况下可能不存在最优解。
通常有两种情况：
(1)根本没有可行解，即给定的约束条件之间是相互排斥的，可行区域为空集；
(2)目标函数没有极值，也就是说在n 维空间中的某个方向上，目标函数值可以无限增大，而仍满足约束条件，此时目标函数值无界。
例：
这个问题的解为 (x1,x2,x3,x4)=(0,3.5,4.5,1)最优值为16。
### 线性规划基本定理
约束条件(8.2)-(8.5)中n个约束以等号满足的可行解称为线性规划问题的基本可行解。
若n>m，则基本可行解中至少有n−m个分量为0，也就是说，基本可行解中最多有m个分量非零。
*线性规划基本定理：*如果线性规划问题有最优解，则必有一基本可行最优解。
上述定理的重要意义在于，它把一个最优化问题转化为一个组合问题，即在(8.2) -(8.5)式的m+n个约束条件中，确定最优解应满足其中哪n个约束条件的问题。
由此可知，只要对各种不同的组合进行测试，并比较每种情况下的目标函数值，直到找到最优解。
Dantzig于1948年提出了线性规划问题的单纯形算法。
单纯形算法的特点是：
（1）只对约束条件的若干组合进行测试，测试的每一步都使目标函数的值增加；
（2）一般经过不大于m或n次迭代就可求得最优解。
### 约束标准型线性规划问题的单纯形算法
当线性规划问题中没有不等式约束(8.2)和(8.4)式，而只有等式约束(8.3)和变量非负约束(8.5)时，称该线性规划问题具有标准形式。
先考察一类更特殊的标准形式线性规划问题。这一类线性规划问题中，每一个等式约束中，至少有一个变量的系数为正，且这个变量只在该约束中出现。
在每一约束方程中选择一个这样的变量，并以它作为变量求解该约束方程。
这样选出来的变量称为左端变量或基本变量，其总数为m个。剩下的n-m个变量称为右端变量或非基本变量。
这一类特殊的标准形式线性规划问题称为约束标准型线性规划问题。
虽然约束标准型线性规划问题非常特殊，但是对于理解线性规划问题的单纯形算法是非常重要的。
任意一个线性规划问题可以转换为约束标准型线性规划问题。
例：
任何约束标准型线性规划问题，只要将所有非基本变量都置为0，从约束方程式中解出满足约束的基本变量的值，可求得一个基本可行解。
单纯形算法的基本思想就是从一个基本可行解出发，进行一系列的基本可行解的变换。
每次变换将一个非基本变量与一个基本变量互调位置，且保持当前的线性规划问题是一个与原问题完全等价的标准线性规划问题。
基本可行解x=(7,0,0,12,0,10)。
### 单纯形法计算步骤如下：
#### **单纯形算法的第1步：**
选出使目标函数增加的非基本变量作为**入基变量**。
查看单纯形表的第1行（也称之为z行）中标有非基本变量的各列中的值。
选出使目标函数增加的非基本变量作为入基变量。
z行中的正系数非基本变量都满足要求。
在上面单纯形表的z行中只有1列为正，即非基本变量相应的列，其值为3。
选取非基本变量x3作为入基变量。
#### **单纯形算法的第2步：**
选取离基变量。
在单纯形表中考察由第1步选出的入基变量所相应的列。
在一个基本变量变为负值之前，入基变量可以增到多大。
如果入基变量所在的列与基本变量所在行交叉处的表元素为负数，那么该元素将不受任何限制，相应的基本变量只会越变越大。
如果入基变量所在列的所有元素都是负值，则目标函数无界，已经得到了问题的无界解。
如果选出的列中有一个或多个元素为正数，要弄清是哪个数限制了入基变量值的增加。
受限的增加量可以用入基变量所在列的元素（称为主元素）来除主元素所在行的“常数列”（最左边的列）中元素而得到。所得到数值越小说明受到限制越多。
应该选取受到限制最多的基本变量作为离基变量，才能保证将入基变量与离基变量互调位置后，仍满足约束条件。
上例中，惟一的一个值为正的z行元素是3，它所在列中有2个正元素，即4和3。
min{124,103}=4，应该选取x4为离基变量；
入基变量x3取值为3。
#### **单纯形算法的第3步:**
转轴变换。
转轴变换的目的是将入基变量与离基变量互调位置。
给入基变量一个增值，使之成为基本变量；
修改离基变量，让入基变量所在列中，离基变量所在行的元素值减为零，而使之成为非基本变量。
解离基变量所相应的方程，将入基变量x3用离基变量x4表示为x1−12x2+14x4=3
再将其代入其他基本变量和所在的行中消去x3，
x1+52x2+14x4+2x5=10x6−52x2−34x4+8x5=1
代入目标函数得到z=9+12x2−34x4−2x5
形成新单纯形表
#### **单纯形算法的第4步：**
转回并重复第1步，进一步改进目标函数值。
不断重复上述过程，直到z行的所有非基本变量系数都变成负值为止。这表明目标函数不可能再增加了。
在上面的单纯形表中，惟一的值为正的z行元素是非基本变量x2相应的列，其值为12。
因此，选取非基本变量x2作为入基变量。
它所在列中有惟一的正元素52，即基本变量x1相应行的元素。
因此，选取x1为离基变量。
再经步骤3的转轴变换得到新单纯形表。
新单纯形表z行的所有非基本变量系数都变成负值，求解过程结束。
整个问题的解可以从最后一张单纯形表的常数列中读出。
目标函数的最大值为11；
最优解为：x∗=(0,4,5,0,0,11)。
### 例子
上述方法严格按照步骤，虽然思路清晰，但过程冗长，不方便操作。下面用表格形式的单纯形方法来解上题。
在引进松弛变量，并将问题转化成标准形式后，建立初始单纯形表：
上述题目求解的是目标函数的min,下面再分享一道求解目标函数max的情况。两种情况极其类似，只是在判别式选择上的差异。
数理统计中的区间估计
最后更新于：2022-04-01 16:06:54
## 区间估计
用点估计θ^(X1,X2,…,Xn)来估计总体的未知参数 θ，一旦我们获得了样本观察值 (x1,x2,…,xn)，将它代入θ^(X1,X2,…,Xn)，即可得到θ的一个估计值。这很直观，也很便于使用。但是，点估计值只提供了θ的一个近似值，并没有反映这种近似的精确度。同时，由于θ本身是未知的，我们也无从知道这种估计的误差大小。因此，我们希望估计出一个真实参数所在的范围，并希望知道这个范围以多大的概率包含参数真值，这就是参数的区间估计问题。
### 定义
设θ为总体ξ的未知参数,ξ1,ξ2,…,ξn为ξ的一个子样，T1(ξ1,ξ2,…,ξn),T2(ξ1,ξ2,…,ξn) 为两个统计量。对于任意给定的α(0<α<1)，若T1,T2满足
P{T1≤θ≤T2}=1−α
——(4)
则称随机区间[T1,T2]为θ的置信水平为1−α的区间估计，α为显著性水平，T1,T2分别称为置信下限和置信上限.
注意：也称T2−T1为该区间估计的精度。
值得注意的是，置信区间(θ^1,θ^2)是一个随机区间，对于给定的样本(X1,X2,…,Xn)， 可能包含未知参数(θ^1,θ^2)，也可能不包含θ。但（4）表明，在重复取样下，将得到许多不同的区间θ^1(x1,x2,…,xn)、θ^2(x1,x2,…,xn)，根据贝努利大数定律，这些区间中大约有100(1−α) 的区间包含未知参数θ 。
置信度表示区间估计的可靠度，置信度1−α越接近于1越好。区间长度则表示估计的范围，即估计的精度，区间长度越短越好。当然，置信度和区间长度是相互矛盾的。在实际问题中，我们总是在保证可靠度的前提下，尽可能地提高精度。因此区间估计的问题，就是在给定α值的情况下，利用样本(X1,X2,…,Xn)去求两个估计量θ^1和θ^2 的问题。
### 置信区间的含义
以α=0.01为例，此时置信度为99。假设反复抽取样本1000次，则得到1000个随机区间[T1,T2]，在这1000个区间中，包含值的大约有990个，而不包含θ值的大约有10个。
### 构造区间估计的步骤
1.构造一个与θ有关的函数
{U不含其它未知参数已知U的分布
2.对给定的
α(0<α<1)
，求
a,b
使得
P{a≤U≤b}=1−α
3.解不等式
a≤U≤b⇔T1≤θ≤T2
,得到区间
[T1,T2]
### 正态总体均值与方差的区间估计
设ξ～N(a,σ2),ξ1,…,ξn为ξ的一子样
### 单个总体ξ～N(a,σ2)的情形
#### σ2=σ20已知时，求a的区间估计
因为ξ¯是a的最优无偏估计，因此在求a的区间估计时，自然从ξ¯出发来构造一个适合的函数。因为
ξ～N(a,σ20)⇒ξ¯～N(a,σ20n)
令U=ξ¯−aσ0/n√，则U～N(0,1)
对给定的α(0<α<1)，求uα，使得
P{|U|≤uα}=1−α——(∗)
临界值uα可由P{U≤uα}=1−α2，查N(0，1)分布表得到
（*）式变为
P{|ξ¯−aσ/n√|≤uα}=1−α
亦是
P{ξ¯−uασ0n√≤a≤ξ¯+uασ0n√}=1−α
因此a的置信水平为1−α的区间估计为
[ξ¯−uασ0n√,ξ¯+uασ0n√]
不同置信水平1−α下，a的区间估计为
##### 例题
设某种清漆的9个样品，其干燥时间(以小时计)分别为6.0，5.7，5.8，6.5，7.0，6.3，5.6，6.1，5.0 ,
设干燥时间总体服从正态分布N(a,0.62)，求a的置信
水平为0.95的置信区间。
解：
σ2=0.62已知，n=9,α=1−0.95=0.05
取U=ξ¯−aσ0/n√～N(0,1)
由P{|U|≤uα}=1−α可得P{U≤uα}=1−α2=0.975
查标准正态分布表得到uα=1.96，故所求的置信水平为0.95的置信区间为
[X¯¯¯−uασ0n√,X¯¯¯+uασ0n√]=[5.0608，6.392]
#### σ2未知时，求a的区间估计
取T=n−1−−−−−√ξ¯−aS～t(n−1)
对给定的α(0<α<1)，求一个t(n−1)(α)，使得
P{|T|≥tn−1(α)}=α
查t分布表可求得tn−1(α)
P{|n−1−−−−−√ξ¯−aS|≥tn−1(α)}=α
即
P{ξ¯−tn−1(α)Sn−1−−−−−√≤a≤ξ¯+tn−1(α)Sn−1−−−−−√}=α
得到a的置信度为1−α的置信区间为
[ξ¯−tn−1(α)Sn−1−−−−−√,ξ¯+tn−1(α)Sn−1−−−−−√]
##### 例题
设某种清漆的9个样品，其干燥时间(以小时计)分别为6.0，5.7，5.8，6.5，7.0，6.3，5.6，6.1，5.0 ,
设干燥时间总体服从正态分布N(a,σ2)，σ2>0未知，求a的置信
水平为0.95的置信区间。
解：
σ2>0未知，n=9,α=1−0.95=0.05
取T=n−1−−−−−√ξ¯−aS～t(n−1)=t(8)
由P{|T|≥t8(α)}=α=0.05
查t(8)分布表可以得到t8(α)=2.306，故所求的区间估计为
[ξ¯−tn−1(α)Sn−1−−−−−√,ξ¯+tn−1(α)Sn−1−−−−−√]
计算得ξ¯=6.0 S2=0.29
故所求a的区间估计为[5.558,6.442]
#### （总体均值未知时）σ2的区间估计
取χ2=nS2σ2～χ2(n−1)
给定的α(0<α<1)，求一个λ1，一个λ2，使得
P{λ1≤χ2≤λ2}=1−α
P{χ2n−1(1−α2)≤χ2≤χ2n−1(α2)}=1−α
P{χ2n−1(1−α2)≤nS2σ2≤χ2n−1(α2)}=1−α
即
P{nS2χ2n−1(1−α2)≤σ2≤nS2χ2n−1(α2)}=1−α
由此可得σ2的置信水平为1−α的区间估计为
[nS2χ2n−1(1−α2)，nS2χ2n−1(α2)]
##### 例题
设某种清漆的9个样品，其干燥时间(以小时计)分别为6.0，5.7，5.8，6.5，7.0，6.3，5.6，6.1，5.0 ,
设干燥时间总体服从正态分布N(a,σ2)，求σ2>0的置信
水平为0.95的置信区间。
解：
n=9,α=1−0.95=0.05，a未知
取χ2=nS2σ2～χ2(n−1)=χ2(8)
对给定的α(0<α<1)，查χ2(8)表得
χ28(0.975)=2.180,χ28(0.025)=17.535
计算得ξ¯=6.0 S2=0.29
故σ2的置信度为1−α的置信区间为
[nS2χ2n−1(1−α2)，nS2χ2n−1(α2)]
[9×0.2917.535，9×0.292.18]
[0.151,1.211]
因此σ2的置信区间为0.95的区间估计为[0.151,1.211]
### 二个正态总体的情形
#### 二个正态总体均值差a1−a2的区间估计
设ξ1,ξ2,…,ξn1与η1,η2,…,ηn1分别是来自正态总体N(a1,σ21)与N(a2,σ22)的子样，且这两个子样相互独立，ξ¯,η¯分别是这两个子样的均样，s21,s22分别是这两个子样的方差。
因为ξ¯,η¯分别为a1,a2的点估计，故取ξ¯−η¯为a1−a2的点估计。此时ξ¯−η¯服从正态分布，且
E(ξ¯−η¯)=a1−a2D(ξ¯−η¯)=D(ξ¯)+D(η¯)=σ21n1+σ22n2
对总体方差的不同情况可得a1−a2的不同置信区间。
##### 若σ21、σ22都已知
取U=ξ¯−η¯−(a1−a2)σ21n1+σ22n2√～N(0,1)
对于给定的α(0<α<1)，查正态分布表得uα
从而得到a1−a2置信水平为1−α的区间估计
(ξ¯−η¯−uασ21n1+σ22n2−−−−−−−−√,ξ¯−η¯+uασ21n1+σ22n2−−−−−−−−√)
##### 若σ21=σ22=σ2都未知
取T=n1n2(n1+n2−2)n1+n2−−−−−−−−−−√ξ¯−η¯−(a1−a2)n1S21+n2S22√～t(n1+n2−2)
对于给定的α(0<α<1)，由
P{|T|>t(n1+n2−2)(α)}=α
确认t(n1+n2−2)(α)，从而得到a1−a2置信区间为1−α的区间估计是
(ξ¯−η¯±t(n1+n2−2)(α)n1S21+n2S22−−−−−−−−−−−√n1+n2n1n2(n1+n2−2)−−−−−−−−−−−−−−−√)
##### 例题
为比较I, II 两种型号子弹的枪口速度,随机地取I 型子弹10 发, 得到枪口速度的平均值为x¯1=500(m/s)，标准差s1=1.10(m/s)，随机地取II 型子弹20 发, 得到枪口速度的平均值为x¯2=496(m/s)，标准差s2=1.20(m/s)，假设两总体都可认为近似地服从正态分布,且由生产过程可认为方差相等,求两总体均值差a1−a2的一个置信水平为0.95的区间.
解：
由假设两总体的方差相等, 但数值未知
故a1−a2置信度为1−α的置信区间是
(ξ¯−η¯±t(n1+n2−2)(α)n1S21+n2S22−−−−−−−−−−−√n1+n2n1n2(n1+n2−2)−−−−−−−−−−−−−−−√)
n1=10,n2=20,n1+n2−2=28
α−1−0.95=0.05,t0.05(28)=2.048
t(n1+n2−2)(α)n1S21+n2S22−−−−−−−−−−−√n1+n2n1n2(n1+n2−2)−−−−−−−−−−−−−−−√=2.048×10×1.102+20×1.202−−−−−−−−−−−−−−−−−−−√30200×28−−−−−−−−√=0.853
故a1−a2置信水平为1−α的区间估计是
(4−0.853,4+0.853)=(3.147,4.853)
#### 两个总体方差比的置信区间
设两正态总体N(a1,σ21)、N(a2,σ22)的参数都为未知的，子样容量分别为n1,n2，且两个子样相互独立，子样方差分别为S21,S22，求方差比σ21/σ22的置信区间
由n1S21σ21～χ2(n1−1),n2S22σ22～χ2(n2−1)构造
F=n1S21σ21/(n1−1)n2S22σ22/(n2−1)～F(n1−1,n2−1)
对于给定的α(0<α<1)，取λ1,λ2使满足
P{λ1≤F≤λ2}=1−α
令P{F<λ1}=α/2,P{F>λ2}=α/2
即取λ1=F(n1−1,n2−1)（1−α2）,λ2=F(n1−1,n2−1)（α2）
P{F(n1−1,n2−1)(1−α2)}≤(n2−1)n1S21(n1−1)n2S22．σ22σ21≤P{F(n1−1,n2−1)(α2)}=1−α
由此得σ21σ22的置信度为1−α的置信区间为
[(n2−1)n1S21(n1−1)n2S22．1F(n1−1,n2−1)(α2),(n2−1)n1S21(n1−1)n2S22．1F(n1−1,n2−1)(1−α2)]
注意如何查表得到F(n1−1,n2−1)(1−α2)
F(n1−1,n2−1)(1−α2)=1F(n2−1,n1−1)(α2)
##### 例子
研究机器A和机器B生产的钢管的内径，随机抽取机器A生产的管子18只，测得样本方差s21=0.34(mm2),抽取机器B生产的管子13只，测得样本方差s22=0.29(mm2）,设两样本相互独立，且由机器A和机器B生产的钢管的内径分别服从正态分布N(a1,σ21),N(a2,σ22)试求方差比σ21/σ22的一个置信水平为0.90 的置信区间.
解：
现在n1=18,s21=0.34;n2=13,s22=0.29
α=0.10F(n1−1,n2−1)(α2)=F(17,12)(0.05)=2.59F(n1−1,n2−1)(α2)=F(17,12)(0.95)=1F(12,17)(0.05)=12.38
于是得σ21σ22的置信度为0.90的置信区间为
(73.4464.09×2.59,73.4464.09×12.38)=(0.44,2.73)
### 练习
假设随机变数X～Ｎ(a,2.8)，现有X的10个观察值X1,…,X10，已知X¯¯¯=110∑10i=1xi=1500,
**1)求a的置信度为0.95置信区间**
解1）
由于σ2=2.82已知，故选U=X¯−aσ/n√～N(0,1)
由α=0.5⇒uα=1.96，a的置信度为0.95置信区间为
[X¯¯¯−uασn√,X¯¯¯+uασn√]=[1500−1.96×2.810−−√,1500+1.96×2.810−−√]=[1498.3,1501.7]
**2)要想使0.95的置信区间长度l小于1，观察值个数n最少应去多少？**
解2）
l=(X¯¯¯+uασn√)−(X¯¯¯−uασn√)=2uασn√
要使a的置信度为0.95置信区间的长度小于1,
即
2uασn√<1⇒2×1.96×2.8n<1
⇒n>(2×1.96×2.8)2=120.47
所以观察值个数n最少应取121
**3）如果样本容量n=100，那么区间(X¯¯¯−1,X¯¯¯+1)作为a的区间估计，其置信度是什么？**
解3）
置信区间若为(X¯¯¯−1,X¯¯¯+1)，则l=2
即有等式2uασn√=2⇒uα=100√2.8=3.57
P(|U|≤uα)=P(|U|≤3.57)=2Φ(3.57)−1=2×0.9998−1=0.9996
置信度为0.9996
数理统计中的点估计
最后更新于：2022-04-01 16:06:52
• 统计推断的基本问题有二：估计问题，和假设检验问题.
• 本章讨论总体参数的点估计和区间估计.理解这两种估计的思想，掌握求参数估计量的方法和评判估计量好坏的标准.
### 点估计
### 问题的提出
设灯泡寿命T～N(μ,σ2)，但参数μ和σ2未知. 现在要求通过对总体抽样得到的样本，构造两样本函数分别μ和σ2作出估计,称为估计量, 记为μ′和σ2′,代入观察值x=(x1,…,xn)，得相应估计值.在不致混淆时统称为**估计**.
借助于总体的一个样本，构造适当的样本函数来估计总体S未知参数的值的问题称为参数的点估计问题.
• 两种常用的构造估计量的方法: 矩估计法和极大似然估计法.
### 矩估计
#### 思想与方法
设总体k阶矩存在，
对于连续型总体X，它的m阶原点矩为
μk:=E(Xk)=∫+∞−∞xkdF(x,θ)
若X为离散型的，则
μk:=E(Xk)=∑i=1nxkF(x,θ)
这里θ为未知参数向量. 可见μk是θ的函数，改记为μk(θ) .
设测得10个灯泡寿命(失效时间)分别为
166，185，232，242，264，268，270，275，285，312
(小时）.
那么自然想到平均寿命为
(166+185+...+312)/10=249.9（小时）
即用样本均值的观测值x¯来估计总体的平均寿命(期望寿命) μ
即
μ^=X¯¯¯=1n(X1+X2+⋯+Xn)
对μk(θ),k阶样本原点矩为
μk^(θ)=Mk=1n(Xk1+Xk2+⋯+Xkn)
这就是矩估计的思想：
用样本的k阶矩作为总体k阶矩的估计量.如果未知参数有m个，则可建立m个方程
μ^k(θ1,…,θm)=Mk,k=1,…,m
(如总体μm存在). 从中解出θj=θj(X1,X2,…,Xn), 改记为θ^,并作为θj的估计量. 称这种估计量为**矩估计量**, 相应观察值称为**矩估计值**.
由上一篇文章讲得经验df函数性质可以知道
样本矩几乎处处收敛于总体矩，
• 样本矩的连续函数也几乎处处收敛于总体矩的相应的连续函数,它保证:几乎每次从容量足够大的样本观测值，都可得到相应总体参数的近似值.
#### 例题1
设总体X的二阶矩存在，求总体X的期望和方差的矩估计量.
解：
m=2，可得
（将μ^和σ^2当做未知量，将Xi当做已知量，解方程组）
解得
结论：不论总体有什麽样的分布，只要它的*期望*和*方差*存在，则它们的矩估计量都分别是其样本均值和样本的二阶中心矩.
为突出是矩估计量，也常加下标M，例如μ^M
#### 例题2
设总体X～U(0,θ), θ未知，(X1,…,Xn)是一个样本, 试求θ的矩估计量.
解:
直接由上例结果，令解得θ的矩估计量
θ^M=2X¯¯¯
#### 例题3
设总体 ，即 具有概率密度
这里a，b为未知参数，(X1,X2,…,Xn)为抽自X的简单随机样本
由于E(X)=a+b2, D(X)=(b−a)212
令
由此可解得a和b的矩估计为
a^=X¯¯¯−3√Snb^=X¯¯¯+3√Sn
其中S2n=1n∑ni=1(Xi−X¯¯¯)2
### 极大似然估计法
#### 思想和方法
假设在一个罐中放着许多黑球和白球，并假定已知它们的数目之比为 ，但不知哪种颜色的球多。如果我们有放回地从罐中抽取 个球，则其中的黑球数 服从二项分布：
P(X=k)=Ck3pkq3−k,k=0,1,2,3
其中p=罐中黑球数目罐中全部球的数目,q=1−p，由假设知道p可能取值为14或34.
现在根据样本中的黑球数，来估计未知参数 ，也就是说在14和34之间作一选择。对抽样的四种可能结果计算出相应的概率：
从表1中可见，如果样本中的黑球数为0，那么具有X=0的样本来自p=14的总体的可能性比来自p=34的总体的可能性大，这时应当估计p为14而不是 34。如果样本中黑球数为2，那么具有X=2的样本来自p=34的总体的可能性比来自p=14的总体的可能性大，这时应当估计p为34而不是14。从而可以选择估计量：
也就是说根据样本的具体情况来选择估计量 p^，使得出现该样本的可能性最大。
一般地，若总体X具有概率密度p(x,θ1,θ2,…,θk)，其中θ1,θ2,…,θk为未知参数，又设(x1,x2,…,xn)是样本的一组观察值，那么样本(X1,X2,…,Xn)落在点(x1,x2,…,xn)的邻域内的概率为∏ni=1p(xi;θ1,θ2,…,θk)dxi，它是θ1,θ2,…,θk的函数。
最大似然估计的直观想法是：既然在一次试验中得到了观察值(x1,x2,…,xn)，那么我们认为样本落入该观察值(x1,x2,…,xn) 的邻域内这一事件应具有最大的可能性，所以应选取使这一概率达到最大的参数值作为参数真值的估计。记
离散型时θ应使
L(x,θ):=L(x1,…,xn;θ)=∏i=1np(xi;θ)
最大;
连续型时θ应使
f(x1,…,xn;θ)dx1…dxn=∏i=1nf(xi;θ)dxi
也即, 使L(x,θ)=∏ni=1f(xi;θ)最大.
称L(x,θ)为样本的似然函数.
这样得到的估计值, 称为参数θ的极大似然估计值, 而相应的统计量称为参数θ的极大似然估计量.
求θ的最大似然估计就是求似然函数L(x;θ)的最大值点的问题。
如L(x;θ)关于θ可微, 这时也可以从方程
解出. (1.12)和(1.13)都称为**似然方程**.
由于在许多情况下，求lnL(x;θ)的最大值点比较简单，而且lnx是x的严格增函数，因此在lnL(x;θ)对θi（i=1,2,…,k）的偏导数存在的情况下， 可由(1.13)式求得.
解这一方程组，若lnL(x;θ)的驻点唯一，又能验证它是一个极大值点，则它必是lnL(x;θ)的最大值点，即为所求的最大似然估计。但若驻点不唯一，则需进一步判断哪一个为最大值点。还需指出的是，若 lnL(x;θ)对 θi（i=1,2,…,k）的偏导数不存在，则我们无法得到方程组（1.13），这时必须根据最大似然估计的定义直接求L(x,θ)的最大值点。
有时我们需要估计g(θ1,θ2,…,θk)，如果θ^1,θ^2,…,θ^k分别是θ1,θ2,…,θk 的最大似然估计，且g(θ1,θ2,…,θk)为连续函数，则g(θ^1,θ^2,…,θ^k) 是g(θ1,θ2,…,θk) 的最大似然估计。
#### 例题1
设X～N(μ,σ2)， x1,…,xn 为一个样本值求未知参数μ和σ2的极大似然估计量.
解：
似然函数为
它的对数为
解对数似然方程组（见1.13）：
可得
由于对数似然方程组有唯一解，且它一定是最大值点，于是 μ和σ2的最大似然估计为
#### 例题2
求事件发生的概率 的最大似然估计。
解：
若事件A发生的概率P(A)=p，定义随机变量
则X～B(1,p)，其概率分布为
P(X=xi)=pxi(1−p)1−xi,xi=0,1
设(X1,X2,…,Xn)为抽自X的样本，则似然函数为
由对数似然方程
解得
注意到∑ni=1xi≤n，容易验证d2lnLdp2在x¯处取得负值，于是x¯是lnL的最大值点，因而p的最大似然估计为p^=X¯¯¯
于是我们有结论：频率是概率的最大似然估计。
#### 例题3
设总体 X～U[a,b]，(X1,X2,…,Xn) 为抽自X的样本，求未知参数a,b的最大似然估计。
解：
由于X的密度函数为
因此似然函数为
显然，作为a,b的二元函数，L是不连续的。这时我们不能用方程组（1.13）来求最大似然估计，而必须从最大似然估计的定义出发来求L的最大值点。
为使L达到最大，b-a应尽量地小，但b又不能小于max{x1,x2,…,x3}，否则L(x1,x2,…,x3;a,b)=0 ；类似地，a 又不能大于min{x1,x2,…,x3}。因此a,b的最大似然估计为
## 估计的优良性准则
同一个未知参数，可以有几种不同的估计，这时就存在采用哪一种估计的问题。另一方面，对同一个参数，用矩估计法和最大似然估计法，即使得到同一个估计，也存在衡量该估计量优劣的问题。设θ为未知参数， θ^是θ的估计，直观上讲，θ^与θ越接近越好，为了度量θ^与θ的接近程度，我们可以采用|θ^−θ|作为衡量的标准，但由于θ^(X1,X2,…,Xn)依赖于样本，它本身是随机变量，而θ又是未知的，因此很难采用。下面我们从不同的角度，提出几种衡量估计优劣的标准。
### 一致性
#### **定义1：**
设θ^(X1,X2,…,Xn)是总体X分布的未知参数θ的估计量，若θ^依概率收敛于θ，即对任意的ε>0，
limn→∞P(|θ^−θ|<ε)=1
则称θ^是θ的一致估计。
满足一致性的估计量 θ^，当样本容量n 不断增大时， θ^观察值能越来越接近参数真值 。这很容易理解，当样本容量n越大时，信息越多，当然估计就越准确。
由大数定律知，样本均值X¯¯¯是总体均值 μ（即 E(X)）的一致估计。还有，样本修正方差S2是总体方差σ2（即D(X)）的一致估计。
#### 例题1
若总体X服从正态分布N(μ,σ2)， (X1,X2,…,Xn)是来自总体 X的容量为n的样本，EXi=μ ，DXi=σ2 ，i=1,2,…,n ，则由大数定律知，X¯¯¯依概率收敛于μ，即
limn→∞P(|X¯¯¯−μ|<ε)=1
也即未知参数μ的最大似然估计或矩估计μ^=X¯¯¯是μ的一致估计。
#### 例题2
若总体X服从泊松分布P(λ)，(X1,X2,…,Xn) 是从总体X中抽取的容量为n的样本，且EXi=λ ，DXi=λ ，i=1,2,…,n，则 X¯¯¯依概率收敛于 λ，故未知参数λ 的最大似然估计或矩估计 λ^=X¯¯¯是 λ的一致估计。
#### 例题3
若总体X服从0-1分布，P(X=1)=p,0<p<1， (X1,X2,…,Xn) 是从X中抽取的容量为n的样本EXi=p ，DXi=p(1−p) ，i=1,2,…,n则 X¯¯¯依概率收敛于 p，故未知参数p 的最大似然估计或矩估计 p^=X¯¯¯是 p的一致估计。
### 无偏性
设θ为总体分布的未知参数，θ^(X1,X2,…,Xn) 是θ的一个估计，它是一个统计量。对于不同的样本 (X1,X2,…,Xn)，θ^(X1,X2,…,Xn) 取不同的值。
#### 定义2
如果θ^(X1,X2,…,Xn)的均值等于未知参数θ ，即E[θ^(X1,X2,…,Xn)]=θ 对一切可能的θ成立 ————（3）
则称θ^(X1,X2,…,Xn)为θ 的无偏估计 。
无偏估计的意义是：用θ^(X1,X2,…,Xn)去估计未知参数 θ，有时候可能偏高，有时候可能偏低，但是平均说来等于未知参数 θ。
在（3）式中，对一切可能的θ ，是指在每个具体的参数估计问题中，参数θ取值范围内的一切可能的值。例如，若θ是正态总体N(μ,σ2)的均值μ，那么它的一切可能取值范围是 (−∞,+∞)。若θ是方差 σ2，则它的取值范围为(0,+∞)。我们之所以要求（3）对一切可能的θ 都成立，是因为在参数估计中，我们并不知道参数的真值。因此，当我们要求一个估计量具有无偏性时，自然要求它在参数的一切可能取值范围内处处都是无偏的。
#### 例题1
设(X1,X2,…,Xn) 是抽自均值为μ的总体的样本，考虑μ的如下估计量：
μ^1=X1μ^2=X1+X22μ^3=X1+X2+Xn−1+Xn4
假设n≥4
因为EXi=μ，容易验证Eμ^i=μ,i=1,2,3 ，所以μ都是 的的无偏估计，但是
μ^4=2X1μ^5=X1+X23
都不是μ的的无偏估计。
对于任一总体 X，由于EX¯¯¯=μ，所以X¯¯¯ 是μ的的无偏估计，但由于 ES2n=E[1n∑ni=1(Xi−X¯¯¯)2]=n−1nσ2，故S2n不是总体方差σ2的无偏估计，而修正的样本方差 是总体方差 S2n=1n−1∑ni=1(Xi−X¯¯¯)2的无偏估计。
若 θ^是 θ的估计， g(θ)为θ 的实函数，通常我们总是用g(θ^) 去估计g(θ) ，但是值得注意的是，即使 Eθ^=θ，也不一定有E(g(θ^))=g(θ) 。
#### 例题2
修正样本方差的标准差S不是总体标准差σ的无偏估计。
事实上，由于 σ2=E(S2)=DS2+[ES]2≥[ES]2，从而σ≥ES ，即 S不是σ的无偏估计。
若θ的估计θ^不是无偏的，但当n→∞ 时，Eθ^→θ ，则称θ^ 是θ的**渐近无偏估计**。显然，样本方差S2n是总体方差的一个渐近无偏估计。
无偏性对估计量而言是很基本的要求，它的直观意义是没有系统误差。由上例知，对于一个未知参数，它的无偏估计可以不止一个。那么，怎么来比较它们的好坏呢？我们很自然地想到，一个好的估计量应该方差比较小，只有这样才能得到比较稳定的估计值。
### 有效性
#### 定义3
设θ^1(X1,X2,…,Xn)和θ^2(X1,X2,…,Xn)均为参数θ的无偏估计，如果
Dθ^1<Dθ^2
则称**θ^1较θ^2有效**。当θ^(X1,X2,…,Xn)是所有无偏估计中方差最小时，称θ^(X1,X2,…,Xn)
为最小方差无偏估计。
#### 例题
设(X1,X2,…,Xn) 是来自总体X的容量为n的样本，证明总体均值μ （即 EX）的估计量μ^1=X¯¯¯比μ^2=∑ni=1aiXi有效，其中ai≥0,i=1,2,…,n且∑ni=1ai=1 。
证明
由于 Eμ^1=μ，Eμ^2=E(∑ni=1aiXi)=μ∑ni=1ai=μ ，所以μ^1,μ^2均是μ的无偏估计。
又
从而
所以X¯¯¯比∑ni=1aiXi有效。
由上例和一致性知，样本均值X¯¯¯是总体均值 μ（即EX）的一致最小方差无偏估计。同样还可以证明，样本修正方差S2是总体方差σ2 （即DX ）的一致最小方差无偏估计。
数理统计
最后更新于：2022-04-01 16:06:50
当研究并解决一个实际问题时， 我们会
遇到下面问题：
• 1. 这个随机现象可以用什么样的分布律
来刻划，这种分布律的选用合理吗？
• 2. 所选用的这一分布律的参数是多少？
如何估计和确定这些参数？
如何利用数据资料，作出尽可能精确可
靠的统计结论(统计推断):
1） **估计**——从局部观测资料的统计特征，推断总体的特征(分布与矩)；
2）**假设检验**——依据抽样数据资料，对总体的某种假设作检验，从而决定对此假定是拒绝抑或接受.
### 数理统计的基本概念
**总体**：研究对象全体; 也称母体, 记作S.
**样本**：总体中抽出作观测的个体;也称子样，记ω
**样本容量**：抽取的个体数目;也称样本大小.
### 例子
随机抽5支,得寿命数据(称为观察[测]值）：
725，520，683，992，742.(小时)
一般记为，x1 x2 x3 x4 x5.
又抽5支， x′1 x′2 x′3 x′4 x′5.
再抽5支， x′′1 x′′2 x′′3 x′′4 x′′5 .
…… ……
如此继续. 各组观察值彼此不同.
如此继续. 每组中的第一支灯的寿命，
也彼此不同. 这样，泛指所抽取的第一支荧光灯的寿命应是一个rv，记为
X1 . 同样第二支的寿命是rv X2 ，…
如此得一组rv : X1,X2,X3,X4,X5
称为大小为5的样本.
一般地则有大小（容量）为n 的样本,称x1,x2,...,xn为**样本观察值**[*现实*].
抽取的样本如能切实保证其随机性，那么应该彼此独立，且能反映总体的随机规律性，即所有样本彼此独立且与总体同分布.
这样的样本，我们称之为**简单样本**. 这种抽样方法，叫**简单抽样**.
注意，在有限总体中，各观察结果可能不独立.
### 样本的数字特征与分布
最简单又方便的样本函数g(X1,…,Xn)是Xi们的一次和二次的线性组合.
由于样本“平等”，线性组合中应有相等的权系数.
**一次时:**样本的算术平均值X¯¯¯;
**二次时:**中心化后的样本二阶中心矩S2n.
设X1,…,Xn为总体S的大小为n的样本， 分别称
X¯¯¯=1n∑i=1nXi S2=1n−1∑i=1n(Xi−X¯¯¯)2
为
**样本均值**
和
样本方差[（样本方差除以n-1的原因）](http://www.dutor.net/index.php/2009/10/sample-variance/)
，而依次称
Mk=1n∑i=1nXki S2n=1n∑i=1n(Xi−X¯¯¯)2
为
**样本的k阶矩**
和
**样本的二阶中心矩**
.
记号：**总体k阶矩**：
μk=EXk∫+∞−∞xkdFX(x)
**总体的k阶中心矩**
：
σk=∫+∞−∞(x−EX)kdFX(x)
μ=μ1,σ2=σ2
.
注意
1)M1=X¯¯¯,S2n没叫样本方差.
2) 比较总体的期望μ、方差σ2与矩μk:
1. 样本的均值、方差及k阶矩等都是rv，并且因n有限而总是存在的.
2. 总体的期望、方差及k阶矩等不一定存在.且即便存在，也是实数值， 而非rv.
3.代入观察值, 有相应的**样本矩的观察值**x,m以及s2 等.
性质 如果总体k阶矩存在，则样本的k阶矩的数学期望等于总体的k阶矩，而当n趋于无穷时，样本的k阶矩以概率收敛到总体的k阶矩，即
### 顺序统计量与经验df
仍从观察值出发设法求总体分布. 以五支荧光灯寿命数据725,520,683,992,742为例，构造
其df 函数(如后图)称为**经验df函数**.
设{xi}观察值重新依序排列为{x(n)}: x(1)≤x(2)≤⋯≤xn
令
称为由{xi}决定的**经验df**, 简记为F∗n(x).
将以从小到大为序重新排列的一个样本，称为**顺序统计量**，专记为x(1) x(2) … xn
下面一个非常重要的定理确立经验df 的重要地位. 此定理保证，几乎由每一组观察值得到的经验df，只要n足够大，都可作为总体df的近似. 定理中一致收敛性和几乎处处收敛性，给了我们充分的自由.从而由样本去找总体df，理论上有一个完满的解决.
limn→∞F∗n(x)=F(x)
### 抽样分布与统计量
#### 正态总体常用的样本函数
1.设总体S～N(μ,σ2). 则
样本均值X¯¯¯～N(μ,σ2n)，从而
Z:=X¯¯¯−μσ/n√～N(0,1)
2.K2n:=∑n1(Xi−μσ)2的分布χ2(n)
K2n是n个独立的标准正态变量的平方和,称n个独立的标准正态变量的平方和的分布为自由度为n的[χ2分布](http://baike.baidu.com/link?url=Nu_ktFPjY7pDAtSiJt5IXx6pOjijIZhxJp1RvQ1yFDskdSmu1gnhk6QLk9JRPZqXIorAfySMJqg2yQCo4Fo_mq).
3.(n−1)S2σ2～χ2(n−1)
样本均值与样本方差独立, 且
K2=(n−1)S2σ2=∑1n(Xi−X¯¯¯σ)2 ～ χ2(n−1)
在K2n=∑n1(Xi−μσ)2中用X¯¯¯易μ得K2.
4.T:=X¯−μS/n√～ t(n−1)
Z:=X¯−μσ/n√～N(0,1)中如σ未知，S2是σ2的无偏估计，自然用S代替Z中的σ引入T
如果Z ～Ｎ（０，１），Y～χ2(n)且独立，则称
t=ZY/N−−−−√～t(n)
即自由度n的[t分布](http://baike.baidu.com/link?url=scKS9Aozzu4_3ydPC18Kg4S5jrD4nkyvensgS2exsIZW-SgpuEXxDOw64SgbKV3UjdE3CNKH7bvVqfcMxD_jfq).
5.Fnm:=S21σ22S22σ21～F(n−1,m−1)
如果X～χ2(n),Y～χ2(n)，且两者相互独立，则称F=χ2(n)/nχ2(m)/m～Ｆ(n,m)
为自由度为n,m的[F分布](http://baike.baidu.com/view/1173064.htm)
#### 性质
• t 分布是对称的,且n→∞极限为正态(n≥30时近似的效果就很好) .
• t 分布只有k<n阶矩.
• κ2分布和F分布不对称，且x<0 时为0.
• κ2 分布的可加性：设U 与V 独立，且分别~κ2(n)和κ2(m)，则U＋V～κ2(n+m).
对给定的实数α∈(0,0.5), 使
P(X>y)=∫∞yfX(x)=α
成立的点y, 称为X 或其分布的上百分位α点. 特别对N(0,1)、t(n)、κ2(n)和F(n,m)分布, 分别记为
zα,tα(n),χ2α(n),Fα(n,m)
使
P(X>y)=∫∞yfX(x)=1−α
成立的点y, 称为X 或其分布的下百分位α点. 特别对N(0,1)、t(n)、κ2(n)和F(n,m)分布, 分别记为
z1−α,t1−α(n),χ21−α(n),F1−α(n,m)
百分位点的值，可由表查得.
#### 例题：
##### 例题1：
设X1,X2,…,Xn, 是来自总体X～N(0,σ2)的简单随机样本，求统计量
∑10i=1(−1)iXi∑20i=11X2i−−−−−−−−√
的分布。
解：
由题意可知Xk～N(0,σ2)可得
∑10i=1(−1)iXi～N(0,10σ2)
∑10i=1(−1)iXi /10−−√σ～N(0,1)
又因为∑20i=11(X2iσ)～χ2(10)
故由t分布定义可得
∑10i=1(−1)iXi∑20i=11X2i−−−−−−−−√ = ∑10i=1(−1)iXi10−−√σ(∑20i=11(X2i/10)σ)−1～t(10)
##### 例题2：
设X1,X2,…,Xn+1是正态总体的简单样本，前面容量为n的样本均值和样本二阶中心矩分别为X¯¯¯ 和S2n
试求下列样本函数的分布
1)(n−1)(X1−μ)2 / ∑ni=2(Xi−μ)2
2)Xn+1−X¯Snn−1n+1−−−√
解：
1)
(n−1)(X1−μ)2 / ∑ni=2(Xi−μ)2=(Xi−μ)2σ2∑ni=2(Xi−μσ)2n−1
分子服从χ2(1)，分母服从χ2(n−1)
所以整个式子服从F(1,n−1)
2)
Xn+1−X¯Snn−1n+1−−−√
分母部分变成：
S2n(n−1)σ2～χ2(n−1)
分子部分变成：
Xn+1−X¯σ～Ｎ(0,1)
因此原式变成：
Xn+1−X¯σS2n(n−1)σ2√ / n−1√
服从t(n−1)
无约束最优化方法
最后更新于：2022-04-01 16:06:47
## 无约束最优化方法
最优化(Optimization)是工程技术、经济管理、 科学研究、社会生活中经常遇到的问题, 如：结构设计，资源分配，生产计划，运输方案等等。
最优化: 在一定条件下，寻求使目标最大(小)的决策解决优化问题的手段
• 经验积累，主观判断
• 作试验，比优劣
• 建立数学模型，求解最优策略
### 优化问题的数学模型
无约束优化问题标准形式：
求n维决策变量x=[x1,x2,…,xn]T在可行域Ω中
使得目标函数f(x)→min,x⊆Ω⊆Rn
⇒
**可行解**只满足(2)，**最优解**满足(1)(2)
**无约束优化**只有(1)，**约束优化**需要(1)(2)
### 无约束优化的基本方法
### 求解无约束优化的基本思路：[下山法](http://blog.csdn.net/u013007900/article/details/43525733)
#### 基本思想：
和爬山法基本思想一样，在Rn中某一点，确定一个**搜索放向**然后沿着该方向的**移动步长**，得到使目标函数下降的最新点。
#### 迭代步骤：
Step 1 初始化：初始点x0，终止准则等
Step 2 迭代改进：方向dk，步长αk
xk+1=xk+(αd)k 且 f(xk+1)<f(xk)
Step 3 终止检验：得到近优解或k+1⇒k转2> 选择dk,αk使得f下降更快⇒利用不同的算法
### 最速下降法（梯度法）
泰勒级数：
其中Δx可正可负，但必须充分接近于0。
X沿D方向移动步长α后，变为X+αD。由泰勒展开式：
目标函数：
α确定的情况下即最小化：
向量的内积何时最小？当然是两向量方向相反时。所以X移动的方向应该和梯度的方向相反。
**下降方向：**∇fT(xk)dk<0
**最速下降方向：**dk=−∇f(xk)负梯度方向
**迭代改进形式：**xk+1=xk−∇f(xk)
算法特点：在初始状态改进较快，在最优解附近改进缓慢。
求函数f(x1,x2)=x21+4x22在极小点，以初始点X0=(1,1)T。
~~~
#include"matrix.h"//大家自己YY一个矩阵类吧
#include<iostream>
#include<iomanip>
#include<cmath>
#include<limits>
#include<cassert>
using namespace std;
const int SIZE=2;
const double ZETA=0.001;
inline double func(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
return x1*x1+4*x2*x2;
}
inline Matrix<double> gradient(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
Matrix<double> rect(SIZE,1);
rect.put(0,0,2*x1);
rect.put(1,0,8*x2);
return rect;
}
inline Matrix<double> Hesse(Matrix<double> &X){
Matrix<double> rect(SIZE,SIZE);
rect.put(0,0,2);
rect.put(0,1,0);
rect.put(1,0,0);
rect.put(1,1,8);
return rect;
}
int main(int agrc,char *argv[]){
Matrix<double> X(SIZE,1);
X.put(0,0,1);
X.put(1,0,1);
int iteration=0;
double value=func(X);
double newValue=numeric_limits<double>::max();
while(++iteration){
Matrix<double> G=gradient(X);
double factor=((G.getTranspose()*G).get(0,0))/((G.getTranspose()*Hesse(X)*G).get(0,0));
for(int i=0;i<G.getRows();++i)
for(int j=0;j<G.getColumns();++j)
G.put(i,j,G.get(i,j)*factor);
Matrix<double> newX=X-G;
//cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl;
newValue=func(newX);
if(fabs(newValue-value)<ZETA)
break;
else{
X=newX;
value=newValue;
}
//cout<<"本次迭代找到的极小值"<<value<<endl;
}
cout<<"迭代次数"<<iteration<<endl;
cout<<"极小值"<<value<<endl;
return 0;
}
~~~
> 这份算法使用的是利用Hesse矩阵和泰勒公式计算最优步长
若f(X)具有二阶连续偏导，由泰勒展开式可得：
H是f(X)的Hesse矩阵。
可得最优步长：
g是f(X)的梯度矩阵。
此时：
可见最速下降法中最优步长不仅与梯度有关，而且与Hesse矩阵有关。
### 牛顿迭代法
将f(xk+1)在xk点作泰勒展开至二阶项，用d替代dk，有
f(xk+1)=f(xk+d)=f(xk)+∇fT(xk)d+12dT∇2f(xk)d
求d使得f(xk+1)极小⇒右端对d的导数为0⇒∇f(xk)+∇2f(xk)d=0
可得到> **牛顿方程**∇2f(xk)d = −∇f(xk)
**牛顿方向**dk = −(∇2f(xk))−1∇f(xk)
**迭代格式**xk+1 = xk−(∇2f(xk))−1∇f(xk)
其中∇2f(xk)为[Hessen矩阵](http://baike.baidu.com/link?url=E2wdH1BhTKFos1JYd_NhHvkmGAyABsHJDIgSyygvzTuoAs41tvb4OGzzWThbUMwobdBAEPpS7nzia-GQUHHf7a)
可见x的搜索方向是−(∇2f(xk))−1∇f(xk)，函数值在此方向上下降，就需要它与梯度方向相反，即−∇f(xk)(∇2f(xk))−1∇f(xk)<0，所以要求在每一个迭代点上Hessian矩阵必须是正定的。
算法特点：
牛顿法是二次收敛的，并且收敛阶数是2。一般目标函数在最优点附近呈现为二次函数，于是可以想像最优点附近用牛顿迭代法收敛是比较快的。而在开始搜索的几步，我们用梯度下降法收敛是比较快的。将两个方法融合起来可以达到满意的效果。
收敛快是牛顿迭代法最大的优点，但也有致命的缺点：Hessian矩阵及其逆的求解计算量大，更何况在某个迭代点Xk处Hessian矩阵的逆可能根本就不存在（即Hessian
矩阵奇异），这样无法求得Xk+1。
求f(x1,x2)=(x1−2)4+(x1−2x2)2的极小点，初始点取X=(0,3)T
~~~
#include"matrix.h"
#include<iostream>
#include<iomanip>
#include<cmath>
#include<limits>
#include<cassert>
using namespace std;
const int SIZE=2;
const double ZETA=0.001;
inline double func(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
return pow(x1-2,4)+pow(x1-2*x2,2);
}
inline Matrix<double> gradient(Matrix<double> &X){
assert(SIZE==X.getRows());
double x1=X.get(0,0);
double x2=X.get(1,0);
Matrix<double> rect(SIZE,1);
rect.put(0,0,4*pow(x1-2,3)+2*(x1-2*x2));
rect.put(1,0,-4*(x1-2*x2));
return rect;
}
inline Matrix<double> Hesse(Matrix<double> &X){
Matrix<double> rect(SIZE,SIZE);
double x1=X.get(0,0);
double x2=X.get(1,0);
rect.put(0,0,12*pow(x1-2,2)+2);
rect.put(0,1,-4);
rect.put(1,0,-4);
rect.put(1,1,8);
return rect;
}
int main(int agrc,char *argv[]){
Matrix<double> X(SIZE,1);
X.put(0,0,0);
X.put(1,0,3);
int iteration=0;
double value=func(X);
double newValue=numeric_limits<double>::max();
while(++iteration){
Matrix<double> G=gradient(X);
Matrix<double> H=Hesse(X);
Matrix<double> newX=X-H.getInverse()*G;
cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl;
newValue=func(newX);
if(fabs(newValue-value)<ZETA)
break;
else{
X=newX;
value=newValue;
}
cout<<"本次迭代找到的极小值"<<value<<endl;
}
cout<<"迭代次数"<<iteration<<endl;
cout<<"极小值"<<value<<endl;
return 0;
}
~~~
### 拟牛顿法
Hessian矩阵在拟牛顿法中是不计算的，拟牛顿法是构造与Hessian矩阵相似的正定矩阵，这个构造方法，使用了目标函数的梯度（一阶导数）信息和两个点的“位移”（Xk−Xk−1）来实现。
有人会说，是不是用Hessian矩阵的近似矩阵来代替Hessian矩阵，会导致求解效果变差呢？事实上，效果反而通常会变好。
拟牛顿法与牛顿法的迭代过程一样，仅仅是各个Hessian矩阵的求解方法不一样。
在远离极小值点处，Hessian矩阵一般不能保证正定，使得目标函数值不降反升。而拟牛顿法可以使目标函数值沿下降方向走下去，并且到了最后，在极小值点附近，可使构造出来的矩阵与Hessian矩阵“很像”了，这样，拟牛顿法也会具有牛顿法的二阶收敛性。
对目标函数f(X)做二阶泰勒展开：
f(x)=f(xk+1)+(x−xk+1)T∇f(xk+1)+12(x−xk+1)TG(xk+1)(x−xk+1)
其中∇2f(xk) ⇔G(xk+1)
上式两边对x求导：
∇f(x)=∇f(xk+1)+G(xk+1)(x−xk+1)
当x=xk的时候，有
G−1(xk+1)[∇f(xk+1)−∇f(xk)]=xk+1−xk
这里我们用Hk表示在点xk处的Hessian矩阵的**逆矩阵**
Hk+1[∇fd(xk+1)−∇f(xk)]=xk+1−xk ......(5)
(5)式就是**拟牛顿方程**
下面提供两种拟牛顿法
### 非线性最小二乘(Least Square)拟合
问题：
给定(ti,yi),i=1,…n, 拟合一个函数y=f(t,x)，
其中x为待定的参数向量, f 对x非线性。
记误差ri(x)=yi−f(ti,x) r(x)=(r1(x),…,rn(x))T
优化模型：
minxR(x)=rT(x)r(x)
记r(x)的雅克比行列式为J(x)=(∂ri∂xj)n×m
∇R=2J(x)Tr(x)∇2R=2J(x)TJ(x)+2SS=∑ni=1ri(x)∇2ri(x)∇2ri(x)=(∂2ri∂xk∂xl)n×m
牛顿法要计算Hessian矩阵，其中S计算量大
若f 对x线性, 则化为线性最小二乘拟合, 此时
S=0
特定算法考虑如何忽略或近似矩阵S
#### Gauss-Newton算法(GN)：忽略矩阵S
∇R=2J(x)Tr(x)∇2R=2J(x)TJ(x)
牛顿方程：∇2f(xk)d=−∇f(xk)
f用R代替，下降方向dk满足
J(x)TJ(x)dk=−J(x)Tr(x)
GN算法收敛性特点:
S=∑ni=1ri(x) 收敛性依赖f对x的线性程度，即偏差r的大小（”小残差问题“）
#### LM算法: Levenbery (1944), Marquardt (1963)
G-N算法修正
J(x)TJ(x)dk=−J(x)Tr(x)
⇓
(J(x)TJ(x)+αkI)dk=−J(x)Tr(x)
其中αk>0为修正参数，防止JTJ出现病态。
LM算法特点
dk位于GN方向（αk很小）和负梯度方向(αk很大)之间
LM算法的另一种观点
min∥r(xk)+J(xk)(x−xk)∥2
s.t. ∥x−xk∥2≤hk (∗)置信域
其解为
(J(x)TJ(x)+αkI)dk=−J(x)Tr(x)
若GN方向满足约束(*)，αk=0,否则ak>0
置信域算法特点：
先确定最长步长hk，再确定方向dk
一般线性最小二乘法
最后更新于：2022-04-01 16:06:45
### 最小二乘法：
人们对由某一变量t 或多个变量t1,…,tn 构成的相关变量y感兴趣。如弹簧的形变与所用的力相关，一个企业的盈利与其营业额，投资收益和原始资本有关。为了得到这些变量同y之间的关系，便用不相关变量去构建y，使用如下函数模型
ym=f(t1,…,tq;x1,…,xp)
q个相关变量或p个附加的相关变量去拟和。
通常人们将一个可能的、对不相关变量t的构成都无困难的函数类型称作函数模型（如抛物线函数或指数函数）。参数x是为了使所选择的函数模型同观测值y相匹配。（如在测量弹簧形变时，必须将所用的力与弹簧的膨胀系数联系起来）。其目标是合适地选择参数，使函数模型最好的拟合观测值。一般情况下，观测值远多于所选择的参数。
其次的问题是怎样判断不同拟合的质量。高斯和勒让德的方法是，假设测量误差的平均值为0。令每一个测量误差对应一个变量并与其它测量误差不相关（随机无关）。人们假设，在测量误差中绝对不含系统误差，它们应该是纯偶然误差，围绕真值波动。除此之外，测量误差符合正态分布，这保证了偏差值在最后的结果y上忽略不计。
确定拟合的标准应该被重视，并小心选择，较大误差的测量值应被赋予较小的权。并建立如下规则：被选择的参数，应该使算出的函数曲线与观测值之差的平方和最小。用函数表示为：
minx⃗ ∑i=1n(ym−yi)2
用欧几里得度量表达为：
minx⃗ ∥y⃗ m(x⃗ )−y⃗ ∥2
最小化问题的精度，依赖于所选择的函数模型。
### 线性函数模型
典型的一类函数模型是线性函数模型。最简单的线性式是y=x0+x1t，写成行列式，为
minx0,x1∥∥∥∥∥⎛⎝⎜⎜⎜1⋮1t1⋮tn⎞⎠⎟⎟⎟(x0x1)−⎛⎝⎜⎜⎜y1⋮yn⎞⎠⎟⎟⎟∥∥∥∥∥2=minx∥Ax−b∥2
直接给出该式的参数解：
x1=∑ni=1tiyi−n⋅t¯y¯∑ni=1t2i−n⋅(t¯)2和x0=y¯−x1t¯
其中t¯=1n∑ni=1ti，为t值的算术平均值。也可解得如下形式：
x1=∑ni=1(ti−t¯)(yi−y¯)∑ni=1(ti−t¯)2
### 一般线性情况
若含有更多不相关模型变量t1,...,tq，可如组成线性函数的形式
y(t1,…,tq;x0,x1,…,xq)=x0+x1t1+⋯+xqtq
即线性方程组
x0+x1t11+⋯+xjt1j+⋯+xqt1q=y1x0+x1t21+⋯+xjt2j+⋯+xqt2q=y2⋮x0+x1ti1+⋯+xjtij+⋯+xqtiq=yi⋮x0+x1tn1+⋯+xjtnj+⋯+xqtnq=yn
通常人们将tij记作数据矩阵 A，参数xj记做参数矢量x，观测值yi记作b，则线性方程组又可写成：
⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜11⋮1⋮1t11t21ti1tn1⋯⋯⋯⋯t1j⋯t2j⋯tij⋯tnj⋯t1qt2qtiqtnq⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⋅⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜x0x1x2⋮xj⋮xq⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜y1y2⋮yi⋮yn⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟
即 Ax=b
上述方程运用最小二乘法导出为线性平差计算的形式为：
minx∥Ax−b∥2
### 最小二乘法的解
定理：线性方程组Ax=b的LS问题一定有解，且求解LS问题与求解线性方程组的**法方程组等价**
ATAx=ATb
证明：
rank(A)=rank(ATA)⇒rank(ATA)≤rank([ATA,ATb])=rank(aT[A,b])≤rank(AT)=rank(A)⇒rank(ATA)=rank([ATA,ATb])
从而方程组一定有解
f(x)=∥b−Ax∥2∥b−Ax∗∥2=minx∈Rn∥b−Ax∥2⎫⎭⎬⎪⎪⇒
gradf(x∗)=(∂f(x∗)∂x1,∂f(x∗)∂x2,…,∂f(x∗)∂xn)T=0⇒2(ATAx∗−ATb)=0
说明LS解必定是法方程组的解
ATAx∗=ATb⇒AT(Ax∗−b)=0⇒Ax∗−b⊥R(A)={y∈Rm|y=Ax,x∈Rn}⇒∥b−Ax∥2=∥(b−Ax∗)+A(x∗−x)∥2=∥b−Ax∗∥2+∥A(x∗−x)∥2≥∥b−Ax∗∥2
说明法方程组的解必定是LS解
推论：当秩(Am×n)=n时，ATA为对称正定矩阵，最小二乘法问题有唯一解
xLS=(ATA)−1ATb
其中(ATA)−1AT称为A的伪逆矩阵，标记为A†> 证明：
先将b拆成A的值域及其正交补两部分
b=b1+b2
b1=AA†b∈R(A)
b2=(I−AA†)b∈R(A)⊥
所以Ax−b1∈R(A)，可得
∥∥Ax−b∥∥2=∥∥Ax−b1+(−b2)∥∥2=∥∥Ax−b1∥∥2+∥∥b2∥∥2
故当且仅当
x
是
Ax=b1=AA†b
解时，
x
即为最小二乘解，即
x=A†b
。
又因为
N(A)=N(A†A)=R(I−A†A)={(I−A†A)h:h∈Cn}
故
Ax=AA†b
的通解为
x=A†b+(I−A†A)h:h∈Cn
因为
∥∥A†b∥∥2<∥∥A†b∥∥2+∥∥(I−A†A)h∥∥2=∥∥A†b+(I−A†A)h∥∥2,(I−A†A)h≠0
所以
A†b
又是二范数极小的最小二乘解。
### 例题：
求解超定方程组的最小二乘问题：
⎧⎩⎨⎪⎪⎪⎪2x1+4x2=13x1−5x2=3x1+2x2=62x1−x2=7⇒A=⎛⎝⎜⎜⎜⎜23124−521⎞⎠⎟⎟⎟⎟,b=⎛⎝⎜⎜⎜⎜1367⎞⎠⎟⎟⎟⎟
系数矩阵A是列满秩的，古可以求最小二乘法的解！
最小二乘法解为：
最小平方残差：
本例是用法方程组来求解最小二乘法问题的解！但是这样作也可能会引出好多问题！我们提倡用QR分解来求解。
层次分析法(Analytic Hierarchy Process)
最后更新于：2022-04-01 16:06:43
## 层次分析法（Analytic Hierarchy Process ，简称 AHP ）
AHP (Analytic Hierarchy Process)层次分析法是美国运筹学家Saaty教授于二十世纪80年代提出的一种实用的多方案或多目标的决策方法。其主要特征是，它合理地将定性与定量的决策结合起来，按照思维、心理的规律把决策过程层次化、数量化。
该方法自1982年被介绍到我国以来，以其定性与定量相结合地处理各种决策因素的特点，以及其系统灵活简洁的优点，迅速地在我国社会经济各个领域内，如能源系统分析、城市规划、经济管理、科研评价等，得到了广泛的重视和应用。
### 层次分析法的基本思路：先分解后综合
首先将所要分析的问题层次化，根据问题的性质和要达到的总目标，将问题分解成不同的组成因素，按照因素间的相互关系及隶属关系，将因素按不同层次聚集组合，形成一个多层分析结构模型，最终归结为最低层（方案、措施、指标等）相对于最高层（总目标）相对重要程度的权值或相对优劣次序的问题。
运用层次分析法建模，大体上可按下面四个步骤进行：
1. 建立递阶层次结构模型；
2. 构造出各层次中的所有判断矩阵；
3. 层次单排序及一致性检验；
4. 层次总排序及一致性检验。
### 实例
人们在日常生活中经常会碰到多目标决策问题，例如假期某人想要出去旅游，现有三个目的地（方案）：风光绮丽的杭州（P1）、迷人的北戴河（P2）和山水甲天下的桂林（P3）。假如选择的标准和依据（行动方案准则）有5个：景色，费用，饮食，居住和旅途。
则常规思维的方式如下：
常规思维过程=⎧⎩⎨⎪⎪确定这些准则在你心目中各占的比重多大就每一准则将三个地点进行对比将这两个层次的比较判断进行综合，作出选择∗∗∗
### (1)建立层次结构模型
### （2）构造判断矩阵
通过相互比较确定各准则对于目标的权重，即构造判断矩阵。在层次分析法中，为使矩阵中的各要素的重要性能够进行定量显示，引进了矩阵判断标度（1～9标度法） :
| 标度 | 含义 |
|-----|-----|
| 1 | 表示两个元素相比，具有同样的重要性 |
| 3 | 表示两个元素相比，前者比后者稍重要 |
| 5 | 表示两个元素相比，前者比后者明显重要 |
| 7 | 表示两个元素相比，前者比后者极其重要 |
| 9 | 表示两个元素相比，前者比后者强烈重要 |
| 2,4,6,8 | 表示上述相邻判断的中间值 |
| 倒数：若元素i和元素j的重要性之比为aij,那么元素j与元素i的重要性之比为aji=1aij |
|-----|
对于要比较的因子而言，你认为一样重要就是1:1，强烈重要就是9:1，也可以取中间数值6:1等，两两比较，把数值填入，并排列成判断矩阵（判断矩阵是对角线积是1的正反矩阵即可）]
设准则层包含5个准则，景色：C1，费用：C2，居住：C3，饮食：C4，旅途：C5。相对于目标层：选择旅游地，进行两两比较打分。
构造所有相对于不同准则的方案层判断矩阵
### （3）层次单排序
所谓层次单排序是指，对于上一层某因素而言，本层次各因素的重要性的排序。
具体计算是：对于判断矩阵B，计算满足BW=λmaxW的特征根与特征向量。
式中λmax为矩阵B的最大特征根，W为对应于λmax的正规化的特征向量，W的分量ωi即是相应元素单排序的权值。
#### 利用判断矩阵计算各因素C对目标层Z的权重（权系数）
1. 将A的每一列向量归一化得：ω˘ij=aij∑ni=1aij
1. 对ω˘ij按行求和得ω˘i=∑nj=1ω˘ij
1. 将ω˘j归一化ωi=ω˘i∑ni=1ω˘i，ω=（ω1,ω2,…,ωn）T，即为近似特征根（权向量）
1. 计算λ=1n∑ni=1(Aω)iωi，作为最大特征根的近似值。
得到排序结果：ω=(0.588,0.322,0.090)T,λmax=3.009
矩阵与向量的乘积计算
#### 判断矩阵的一致性检验
判断矩阵通常是不一致的，但是为了能用它的对应于特征根的特征向量作为被比较因素的权向量，其不一致程度应在容许的范围内.如何确定这个范围？
1. 一致性指标：CI=λ−nn−1,CI=0时A一致，CI越大，A的不一致性程度越严重
2. 随机一致性指标RI：
| n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 |
|-----|-----|-----|-----|-----|-----|-----|-----|-----|-----|-----|-----|
| RI | 0.00 | 0.00 | 0.58 | 0.90 | 1.12 | 1.24 | 1.32 | 1.41 | 1.45 | 1.49 | 1.51 |
3. 一致性比率（用于确定A的不一致性的容许范围）
CR=CIRI
当CR<0.1时，A的不一致性程度在容许范围内，此时可用A的特征向量作为权向量。
### 例题 解：
### 第一步：自上而下，先求判断矩阵A的最大特征根与特征向量。
A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1214131312117151547123351211351311⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
λmax=5.073对应于λmax的正规化的特征向量为：
ω(2)=(0.263,0.475,0.055,0.099,0.110)T
### 第二步：计算与准则层各准则相关的判断矩阵最大特征跟及权向量：
### 第三步，算出B2,B3,B4,B5的最大特征值分别为：
λmax(2)=3.002,λmax(3)=3,λmax(4)=3.009,λmax(5)=3
所对应的特征向量分别为：
### 第四步，一致性检验
A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1214131312117151547123351211351311⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥a14=3,a43=2⇒a13=a14a43=3×2=6
前述计算得到了最大特征根：λmax=5.073
CI=λmax−nn−1=5.073−55−1=0.01825
经过查表可以得到平均随机一致性指标RI，从而可检验矩阵一致性：
CR=CIRI=0.018251.12=0.016295<0.1
同理，对于第二层的景色、费用、居住、饮食、旅途五个判断矩阵的一致性检测均通过。
以W(3)k为列向量构成矩阵：
ps:这一次就不附赠代码啦~
龙格-库塔(Runge-Kutta)方法数学原理及实现
最后更新于：2022-04-01 16:06:40
## 龙格-库塔(Runge-Kutta)方法
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。
对于一阶精度的欧拉公式有：
yi+1=yi+hki
其中h为步长，则yi+1的表达式与y(xi+1)的Taylor展开式的前两项完全相同，即**局部截断误差**为O(h2)。
当用点xi处的斜率近似值k1与右端点xi+1处的斜率k2的算术平均值作为平均斜率k∗的近似值，那么就会得到二阶精度的改进欧拉公式：
yi+1=yi+h(k1+k2)
其中k1=f(xi,yi) ， k2=f(xi+h,yi+hk1)
依次类推，如果在区间[xi,xi+1]内多预估几个店上的斜率值k1,k2,…,km，并用他们的加权平均数作为平均斜率k∗的近似值，显然能够构造出具有很高精度的高阶计数公式。
上述两组公式在形式删过的共同点：都是用f(x,y)在某些点上值得线性组合得出y(xi+1)的近似值yi+1，且增加计算的次数，可以提高截断误差的阶，他们的误差估计可以用f(x,y)在xi处的Taylor展开来估计。
于是可考虑用函数f(x,y)在若干点上的函数值的线性组合老构造金斯公式，构造时要求近似公式在f(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使金斯公式达到所需要的阶数。既避免求高阶导数，又提高了计算方法精度的阶数。或者说，在[xi,xi+1]这一步内计算多个点的斜率值，若够将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格-库塔(Runge-Kutta)方法。
一般的龙格-库塔法的形式为
称为P阶龙格-库塔方法。
其中ai,bij,cj为待定参数，要求上式yi+1在点(xi,yi)处作Taylor展开，通过相同项的系数确定参数。
当然，经典的龙格-库塔方法是四阶的。也就是在[xi,xi+1]上用四个点处的斜率加权平均作为平均斜率k∗的近似值，构成一系列四阶龙格-库塔公式。具有四阶精度，即局部截断误差是O(h5)。
下面介绍最常用的一种四阶龙格-库塔方法。
设
yi+1=yi+c1K1+c2K2+c3K3+c4K4
这里K1,K2,K3,K4为四个不同点上的函数值，分别设其为
其中c1,c2,c3,c4,a2,a3,a4,b21,b31,b32,b41,b42,b43均为待定系数。
把K2,K3,K4分别在xi点占城h的幂级数，带入线性组合式中，将得到的公式与y(xi+1)在xi点上的泰勒展开式比较，使其两式右端知道h4的系数相等，经过较复杂的解方程过程便可得到关于ai,bij,cj的一组特解。
a2=a3=b21=b32=12
b31=b41=b42=0
a4=b43=1
c1=c4=16
c2=c3=13
带入之后得到
龙格-库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应正对问题的具体特点选择适合的算法。对于光滑性不太好的解，最好采用**低阶算法**而将步长h取小。
龙格-库塔法的C语言实现
~~~
#include "stdio.h"
#include "stdlib.h"
void RKT(t,y,n,h,k,z)
int n; /*微分方程组中方程的个数，也是未知函数的个数*/
int k; /*积分的步数(包括起始点这一步)*/
double t; /*积分的起始点t0*/
double h; /*积分的步长*/
double y[]; /*存放n个未知函数在起始点t处的函数值,返回时,其初值在二维数组z的第零列中*/
double z[]; /*二维数组,体积为n x k.返回k个积分点上的n个未知函数值*/
{
extern void Func(); /*声明要求解的微分方程组*/
int i,j,l;
double a[4],*b,*d;
b=malloc(n*sizeof(double)); /*分配存储空间*/
if(b == NULL)
{
printf("内存分配失败\n");
exit(1);
}
d=malloc(n*sizeof(double)); /*分配存储空间*/
if(d == NULL)
{
printf("内存分配失败\n");
exit(1);
}
/*后面应用RK4公式中用到的系数*/
a[0]=h/2.0;
a[1]=h/2.0;
a[2]=h;
a[3]=h;
for(i=0; i<=n-1; i++)
z[i*k]=y[i]; /*将初值赋给数组z的相应位置*/
for(l=1; l<=k-1; l++)
{
Func(y,d);
for (i=0; i<=n-1; i++)
b[i]=y[i];
for (j=0; j<=2; j++)
{
for (i=0; i<=n-1; i++)
{
y[i]=z[i*k+l-1]+a[j]*d[i];
b[i]=b[i]+a[j+1]*d[i]/3.0;
}
Func(y,d);
}
for(i=0; i<=n-1; i++)
y[i]=b[i]+h*d[i]/6.0;
for(i=0; i<=n-1; i++)
z[i*k+l]=y[i];
t=t+h;
}
free(b); /*释放存储空间*/
free(d); /*释放存储空间*/
return;
}
main()
{
int i,j;
double t,h,y[3],z[3][11];
y[0]=-1.0;
y[1]=0.0;
y[2]=1.0;
t=0.0;
h=0.01;
RKT(t,y,3,h,11,z);
printf("\n");
for (i=0; i<=10; i++) /*打印输出结果*/
{
t=i*h;
printf("t=%5.2f\t ",t);
for (j=0; j<=2; j++)
printf("y(%d)=%e ",j,z[j][i]);
printf("\n");
}
}
void Func(y,d)
double y[],d[];
{
d[0]=y[1]; /*y0'=y1*/
d[1]=-y[0]; /*y1'=y0*/
d[2]=-y[2]; /*y2'=y2*/
return;
}
~~~
ps:如果有时间的话，可能会回过头来加一分解方程的推到吧…
分支限界(Branch and Bound )算法
最后更新于：2022-04-01 16:06:38
## 分支定界 (branch and bound) 算法
是一种在问题的解空间树上搜索问题的解的方法。但与回溯算法不同，分支定界算法采用广度优先或最小耗费优先的方法搜索解空间树，并且，在分支定界算法中，每一个活结点只有一次机会成为扩展结点。
除少数整数规划可以用线性规划的单纯形法直接求解外，一般整数规划必需寻找新的求解方法。这里我们介绍全整数规划的分枝定界方法，它也可以用于求解混合整数规划问题和0-1规划问题。下面我们从一个例子出发来讨论它的思路和步骤。
例1 解下面的整数规划问题
解:
（1）首先我们注意到问题(4.3)的可行解集为图4-2的阴影部分内的整数格子点组成的集合，暂时不考虑整数限制条件(5)。解相应的线性规划(1)~(4)，即(4.3)的松弛问题：
用线性规划的方式得到如下的图案
得最优解
由(1)和(5)可知ILP的最优解
，且必为整数，从而可以断言
。这一过程称为定界，即给出ILP0问题目标函数最优解z*的下界和上界。
（2）其次我们注意到线性规划(4.4)的解
具有小数，但这两个变量在(4.3)中都必须是整数，那就是说必须把小数部分**划掉**。我们注意到，对!
而言(对!)
也是如此)，最终的最优解不会在3和4之间取值，亦即必然有：
这种表达式实际上是将!
在3和4间的小数部分划掉了，把可行域![]
分成了!
和!
，显然这种分法把原来线性规划的解A点**(2.25,3.75)**排除出去了。但没有排除任何整数可行解。这一过程称为**分枝**，即用两个矛盾的约束条件(4.5)分别代入原问题(4.3)形成两个子问题ILP1和ILP2：
解ILP1的松弛问题ILP1得到
解ILP2的松弛问题ILP2得到
（3）修改上下界：从LP1和LP2的解我们知道有
（4）再分枝：
下面我们就是要划掉!
LP1的解中!
的小数部分，增加约束!
对ILP1进一步的分枝，即
求解LP3得：
对于LP4，不难看出，无可行解。
（5）再修改界，此时我们又有!
（6）再分枝，继续对ILP3进行分枝(由于!)
是小数，增限制条件!又得到
求解LP5得到：
求解LP6得到：
至此，所有的子问题都已探明，求解结束。我们得到了ILP0(即原问题)的最优解：
### 分枝定界法的一般步骤如下：
第一步，先不考虑原问题的整数限制，求解相应的松弛问题，若求得最优解，检查它是否符合整数约束条件；如符合整数约束条件，即转下一步。
第二步，定界。在各分枝问题中，找出目标函数值最大者作为整数规划最优值!
的上界记为!
，从已符合整数条件的分枝中，找出目标函数值最大者作为下界，记为
。即 
第三步，分枝。根据对变量重要性的了解，在最优解中选择一个不符合整数条件的
，令
，(
不为整数)则用下列两个约束条件：
其中
表示不超过
的最大整数，分别加入问题形成两个子问题。
第四步，应用对目标函数估界的方法，或对某一分枝重要性的了解，确定出首先要解的某一分枝的后继问题，并解此问题。若所获得的最优解符合整数条件，则就是原问题的解，若不符合整数条件，再回到第二步，并参照第四步终止后继问题。
在上述过程中，要不断应用分枝、定界、估界来进行判断。当我们求解某子问题的松弛问题时，只要出现下列情况之一，该问题就已探明：
1. 松弛问题没有可行解，则原问题也没有可行解；
2. 松弛问题的最优解恰好全取整数，则该最优解也是其对应的子问题的最优解；
3. 松弛问题的最大值小于现有的下界z，则无论其最优解是否取整数值，都将对应的子问题剪枝；
已探明的子问题就不再用分枝了，如果所有的子问题都已探明，则原整数规划的最优解就一定可以求出，或可以判定它无解。
JAVA实现分支定界代码：
~~~
package sy2;
import sy1.*;
/*
* 假设这里所解的整数规划问题的目标函数取的是max
*/
public class FenZhiDingJie {
double A[][]; //原矩阵的系数矩阵
String D[]; //原矩阵的符号矩阵
double b[]; //原矩阵的常数矩阵
public int count=0;
int index=-1;
public double C[];//目标函数的初始系数向量
public double Uz;//目标函数值下界
public double Lz;//目标函数值上界
public double z;//现在的目标函数值
public double X[];//定义最优解
public double zX[];//定义最优解
public double yX[];//定义最优解
public double Z;//定义整数线性规划的最优值
public double Ix[];//定义整数线性规划的最优解
public double As[][];
public double bs[];
public String Ds[];
public double Xs[];
public int lc;
int M,N;
//初始化分支定解法
FenZhiDingJie(double[][] a,double[] B,String[] d,double[] c){
M=B.length;
N=c.length;
A=a;
b=B;
D=d;
C=c;
X=new double[N];
zX=new double[N];
yX=new double[N];
Ix=new double[N];
}
//分支定解过程
public void FZDJ(double[][] a,double[] B,String[] d,double[] c,double x[]){
boolean flag1=true;
boolean flag2=true;
//利用两阶段法解出该整数线性规划的最优值的上界（第一次使用两阶段法）
while(count<10){
if(count==0){
System.out.println("\n第 "+count+" 次迭代");
TwoStepMethod tm=new TwoStepMethod(A,b,D,C);
if(!tm.flag){
break;
}
X=tm.X;
x=X;
Lz=0;
Uz=-tm.z;
System.out.print("该整数线性规划问题对应的松弛线性规划问题的最优解是：x={");
for(int i=0;i<N-1;i++){
System.out.print(tm.X[i]+",");
}
System.out.print(tm.X[N-1]+"}\n");
System.out.println("该整数线性规划问题对应的松弛线性规划问题的最优值是："+Uz);
index=AllInteger(X,index);
if(index==-1){
System.out.println("该整数线性规划问题的最优解即为它的松弛线性规划问题的最优解，如上所示！");
}
//第一次分支
else{
count++;//跳入分支定解环节
}// else*/
}// if
else{//接下来实现继续的分支定解法
//System.out.println(Lz);
double A1[][]=new double[B.length+1][N];
double yA1[][]=new double[B.length+1][N];
double ytA[][]=new double[B.length+1][N];
double tA[][]=tempA(a,B.length,index);
for(int i=0;i<B.length+1;i++){
for(int j=0;j<N;j++){
A1[i][j]=tA[i][j];
yA1[i][j]=tA[i][j];
ytA[i][j]=tA[i][j];
//System.out.print(ytA[i][j]+" ");
}
//System.out.println();
}
double yb1[]=new double[B.length+1];
double b1[]=new double[B.length+1];
double ytb[]=tempb(B,(int)x[index]+1);
double ztb[]=tempb(B,(int)x[index]);
for(int i=0;i<B.length+1;i++){
b1[i]=ztb[i];
yb1[i]=ytb[i];
//System.out.println(b1[i]+" "+yb1[i]);
}
//左分支
String D1[]=new String[B.length+1];
String ztD[]=LeftD(d);
for(int i=0;i<B.length+1;i++){
D1[i]=ztD[i];
}
String yD1[]=new String[B.length+1];
String ytD[]=RightD(d);
for(int i=0;i<B.length+1;i++){
yD1[i]=ytD[i];
//System.out.print(yD1[i]+" ");
}
TwoStepMethod ztm=new TwoStepMethod(tA,ztb,ztD,c);//左分支利用两阶段法求解
TwoStepMethod ytm=new TwoStepMethod(ytA,ytb,ytD,c);//左分支利用两阶段法求解
if(count==1){
As=yA1;
bs=yb1;
Ds=yD1;
Xs=ytm.X;
}
if(ztm.flag){
double zz=ztm.z;
zX=ztm.X;
index=AllInteger(zX,index);
if(index!=-1){
//System.out.println(DingJie(z,Uz));
if(DingJie(zz,Uz)){
//左分支
System.out.print("第"+count+"次迭代得到一个松弛线性规划的最优解：\nzx={");
for(int i=0;i<N-1;i++){
System.out.print(zX[i]+",");
}// for i
System.out.print(zX[N-1]+"}\n");
System.out.println("对应的最优值是："+-zz);
count++;
if(-zz<Lz){
flag1=false;
lc=count;
}
}// if
else{
System.out.println("不存在可行解");
break;
}//else
}else{
System.out.print("第 "+count+" 次迭代得到一个可行解:\nzx={");
for(int i=0;i<N-1;i++){
System.out.print(zX[i]+",");
}//for i
System.out.print(zX[N-1]+"}\n");
System.out.println("对应的最优值是："+-zz);
Lz=-zz;
Z=Lz;//记录可行解对应的目标函数值
Ix=zX;//记录可行解
count++;
//lc=count;
}//else
}// if
else{
flag1=false;
System.out.println(count+++"次迭代无可行解");
}
if(ytm.flag){
double yz=ytm.z;
yX=ytm.X;
index=AllInteger(yX,index);
if(index!=-1){
//System.out.println(DingJie(z,Uz));
if(DingJie(yz,Uz)){
if(DingJie(yz,Uz)){
//右分支
System.out.print("第"+count+"次迭代得到一个松弛线性规划的最优解：\nyx={");
for(int i=0;i<N-1;i++){
System.out.print(yX[i]+",");
}
System.out.print(yX[N-1]+"}\n");
System.out.println("对应的最优值是："+-yz);
//System.out.println(2);
if(-yz<Lz){
flag2=false;
break;
}
//System.out.println(2);
count++;
}//if
else{
System.out.println("不存在可行解");
break;
}
}
}else{
System.out.println("第 "+count+" 次迭代得到一个可行解:\nyx={");
for(int i=0;i<N-1;i++){
System.out.print(yX[i]+",");
}
System.out.print(yX[N-1]+"}\n");
System.out.println("对应的最优值是："+-yz);
Uz=-yz;
Z=Uz;//记录可行解对应的目标函数值
Ix=yX;//记录可行解
}
}// if
else{
flag2=false;
System.out.println(count+++"次迭代无可行解");
}
if(flag1){
index=AllInteger(zX,index);
FZDJ(A1,b1,D1,c,zX);
count++;
}else{
break;
}
index=AllInteger(Xs,index);
FZDJ(As,bs,Ds,c,Xs);
break;
}
}//else
}//while
//确定可行解存在性
public boolean DingJie(double z,double Uz){
boolean flag=true;
//首先保证算出的最优解对应的目标函数值位于整数线性规划问题的目标函数值的上下界中
if(z<=Lz || z>=Uz){
flag=true;//表示该分支继续往下分不可能找到一个可行解
}// if
else{
flag=false;
}
return flag;
}
//判断所给向量的元素是否全为整数如果不是返回对应第一个非整数的下标
public int AllInteger(double X[],int index){
boolean flag=true;
for(int i=0;i<X.length;i++){
if(X[i]-(int)X[i]!=0){//表示最优解中存在分数
flag=false;//表示此解不是该整数规划的最优解
index=i;//记录下最优解中非整数解的位置
break;
}// if
}// for
if(flag){
index=-1;
}
return index;
}
//左分支各参数处理过程
public String[] LeftD(String d[]){
String D[]=new String[d.length+1];
//对D做相应的转换
for(int i=0;i<D.length;i++){
if(i<d.length){
D[i]=d[i];
}
else{
D[i]="<";
}
}
return D;
}
//右分支各参数处理过程
public String[] RightD(String d[]){
String D[]=new String[d.length+1];
//对D做相应的转换
for(int i=0;i<D.length;i++){
if(i<d.length){
D[i]=d[i];
}
else{
D[i]=">";
}
}
return D;
}
//每次分支扩展参数过程
public double[][] tempA(double a[][],int m,int index){
double A[][]=new double[m+1][N];
//对原约束条件添加一个限制
for(int i=0;i<m+1;i++){
if(i<m){
for(int j=0;j<N;j++){
A[i][j]=a[i][j];
}// for
}// if
else{
for(int j=0;j<N;j++){
if(j==index){
A[i][j]=1;
}//if
else
{
A[i][j]=0;
}// else
}//for j
}// else
}//for i
return A;
}
public double[] tempb(double B[],int Lx){
double b[]=new double[B.length+1];
//对b做相应的变换
for(int i=0;i<B.length+1;i++){
if(i<B.length){
b[i]=B[i];
}// if
else{
b[i]=Lx;
}// else
}// for i
return b;
}
}
~~~
~~~
package sy2;
import java.util.Scanner;
import sy1.TwoStepMethod;
/*
* 代码测试程序
*/
public class Text {
public static void main(String args[]){
int M=2,N=2;//M是约束条件个数，N是未知数个数
double A[][]=new double[M][N];
double b[]=new double[M];
String D[]=new String[M];//原符号向量
double C[]=new double[N];//原目标函数的系数向量
Scanner S=new Scanner(System.in);
for(int i=0;i<M;i++){
for(int j=0;j<N;j++){
A[i][j]=S.nextDouble();
//System.out.print(" ");
}
//System.out.println();
}
for(int i=0;i<M;i++){
b[i]=S.nextDouble();
//System.out.print(" ");
}
for(int i=0;i<M;i++){
D[i]=S.next();
//System.out.print(" ");
}
for(int i=0;i<N;i++){
C[i]=S.nextDouble();
//System.out.print(" ");
}
System.out.println("请输出初始矩阵的方程个数"+M);
System.out.println("请输出初始未知数的个数"+N);
System.out.println("请输出初始矩阵的系数矩阵：");
for(int i=0;i<M;i++){
for(int j=0;j<N;j++){
System.out.print(A[i][j]+" ");
}
System.out.println();
}
System.out.println("请输出初始矩阵的常数项：");
for(int i=0;i<M;i++){
System.out.print(b[i]+" ");
}
System.out.println("\n"+"请输出初始矩阵的符号项：");
for(int i=0;i<M;i++){
System.out.print(D[i]+" ");
}
System.out.println("\n"+"请输出目标函数的系数向量：");
for(int i=0;i<N;i++){
System.out.print(C[i]+" ");
}/*
TwoStepMethod tm=new TwoStepMethod(A,b,D,C);
if(tm.flag){
System.out.print("\n请输出通过两阶段法得到的最优解：X={");
for(int i=0;i<N-1;i++){
System.out.print(tm.X[i]+",");
}
System.out.print(tm.X[N-1]+"}\n");
System.out.println("请输出通过两阶段法得到的最优值是："+tm.z);
}*/
System.out.println("\n\n分支定界法过程：");
double A1[][]=new double[M][N];
for(int i=0;i<N;i++){
for(int j=0;j<M;j++){
A1[i][j]=A[i][j];
}
}
double b1[]=new double[M];
for(int i=0;i<M;i++){
b1[i]=b[i];
}
String D1[]=new String[M];
for(int i=0;i<M;i++){
D1[i]=D[i];
}
double x[]=new double[N];
FenZhiDingJie fzdj=new FenZhiDingJie(A, b, D, C);
fzdj.FZDJ(A1, b1, D1, C,x);
System.out.print("该整型线性规划问题的最优解为X={");
for(int i=0;i<N-1;i++){
System.out.print(fzdj.Ix[i]+",");
}
System.out.print(fzdj.Ix[N-1]+"}");
System.out.println("\n该整型线性规划问题的最优值为：Z="+fzdj.Z);
}
}
~~~
模拟退火算法分析
最后更新于：2022-04-01 16:06:36
## 一. 爬山算法 ( Hill Climbing )
介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。
爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如图1所示：假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。
图1
这题是poj2420，原作者认为是模拟退火，但是并不具备模拟退火的概率事件特点，因此我认为是爬山算法，并借此加入
~~~
<span style="font-size:18px;">#include <iostream>
#include <cmath>
using namespace std;
struct point
{
double x,y;
}p[105];
int dir[8][2] = {-1,-1,-1,0,-1,1,0,-1,0,1,1,-1,1,0,1,1};
double getdis(point a, point b)
{
return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}
double allDis(int n , point f)
{
double sum = 0;
for(int i = 0 ; i < n ; i++)
sum += getdis(p[i],f);
return sum;
}
point fermat(int n)
{
double step = 0;
for (int i = 0 ; i < n ; i++)
step += fabs(p[i].x) + fabs(p[i].y);
point f;
f.x = 0;
f.y = 0;
for (int i = 0 ; i < n ; i++)
f.x += p[i].x , f.y +=p[i].y;
f.x /= n;
f.y /= n;
point t;
while(step > 1e-10)
{
for (int i = 0 ; i < 8 ; i++)
{
t.x = f.x + dir[i][0]*step;
t.y = f.y + dir[i][1]*step;
if(allDis(n,t) < allDis(n,f))
f = t;
}
step *=0.7; //步长改动
}
return f;
}
int main(void)
{
int n;
while (cin >> n)
{
for(int i=0; i<n; i++)
cin >> p[i].x >> p[i].y;
double ans = allDis(n, fermat(n));
int t = ans*10;
if (t%10 < 5)
cout << t/10 << endl;
else
cout << t/10+1 << endl;
}
return 0;
}</span>
~~~
## 二. 模拟退火(SA,Simulated Annealing)思想
爬山法是完完全全的贪心法，每次都鼠目寸光的选择一个当前最优解，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以图1为例，模拟退火算法在搜索到局部最优解A后，会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点，于是就跳出了局部最大值A。
模拟退火算法描述：
若J( Y(i+1) )>= J( Y(i) ) (即移动后得到更优解)，则总是接受该移动
若J( Y(i+1) )< J( Y(i) ) (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能趋向稳定）
这里的“一定的概率”的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。
根据热力学的原理，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为：
P(dE) = exp( dE/(kT) )
其中k是一个常数，exp表示自然指数，且dE<0。这条公式说白了就是：温度越高，出现一次能量差为dE的降温的概率就越大；温度越低，则出现降温的概率就越小。又由于dE总是小于0（否则就不叫退火了），因此dE/kT < 0 ，所以P(dE)的函数取值范围是(0,1) 。
随着温度T的降低，P(dE)会逐渐降低。
我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。
关于爬山算法与模拟退火，有一个有趣的比喻：
爬山算法：兔子朝着比现在高的地方跳去。它找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是爬山算法，它不能保证局部最优值就是全局最优值。
模拟退火：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向高处，也可能踏入平地。但是，它渐渐清醒了并朝最高方向跳去。这就是模拟退火。
下面给出模拟退火的伪代码表示。
## 三. 模拟退火算法伪代码
~~~
/*
* J(y)：在状态y时的评价函数值
* Y(i)：表示当前状态
* Y(i+1)：表示新的状态
* r： 用于控制降温的快慢
* T： 系统的温度，系统初始应该要处于一个高温的状态
* T_min ：温度的下限，若温度T达到T_min，则停止搜索
*/
while( T > T_min )
{
dE = J( Y(i+1) ) - J( Y(i) ) ;
if ( dE >=0 ) //表达移动后得到更优解，则总是接受移动
Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
else
{
// 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也
if ( exp( dE/T ) > random( 0 , 1 ) )
Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
}
T = r * T ; //降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
/*
* 若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。若r过小，则搜索的过程会很快，但最终可能会达到一个局部最优值
*/
i ++ ;
}
~~~
hdu 5017 Ellipsoid
这是我第一次接触模拟退火，是2014年西安网络赛的题目
当时想着用计算几何和解析几何做，后来学长说用模拟退火可以做
~~~
#include <cstdio>
#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
const double EPS = 1e-9;
const double INF = 1e18;
const double dx[8] = {1.0, 1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0};
const double dy[8] = {0.0, 1.0, 1.0, 1.0, 0.0, -1.0, -1.0, -1.0};
double a, b, c, d, e, f;
double dis(double x, double y, double z){
return sqrt(x*x + y*y + z*z);
}
double getZ(double x, double y){
double A = c, B = d*y + e*x, C = a*x*x + b*y*y + f*x*y - 1.0;
double delta = B*B - 4.0*A*C;
if(delta < 0.0) return INF;
double z1 = (-B + sqrt(delta)) / (2.0 * A);
double z2 = (-B - sqrt(delta)) / (2.0 * A);
return z1*z1 < z2*z2 ? z1 : z2;
}
void work(){
double x = 0.0, y = 0.0, z = getZ(x, y);
double step = 0.8;
while(step > EPS){
for(int i=0; i<8; i++){
double nx = x + dx[i]*step;
double ny = y + dy[i]*step;
double nz = getZ(nx, ny);
if(nz >= INF) continue;
if(dis(nx, ny, nz) - dis(x, y, z) < 0.0){
x = nx, y = ny, z = nz;
}
}
step *= 0.99;
}
printf("%.7f\n", dis(x, y, z));
}
int main(){
while(scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e, &f) == 6){
work();
}
return 0;
}
~~~
## 四. 使用模拟退火算法解决旅行商问题
旅行商问题 ( TSP , Traveling Salesman Problem ) ：有N个城市，要求从其中某个问题出发，唯一遍历所有城市，再回到出发的城市，求最短的路线。
旅行商问题属于所谓的NP完全问题，精确的解决TSP只能通过穷举所有的路径组合，其时间复杂度是O(N!) 。
使用模拟退火算法可以比较快的求出TSP的一条近似最优路径。（使用遗传算法也是可以的，我将在下一篇文章中介绍）模拟退火解决TSP的思路：
1. 产生一条新的遍历路径P(i+1)，计算路径P(i+1)的长度L( P(i+1) )
2. 若L(P(i+1)) < L(P(i))，则接受P(i+1)为新的路径，否则以模拟退火的那个概率接受P(i+1) ，然后降温
3. 重复步骤1，2直到满足退出条件
产生新的遍历路径的方法有很多，下面列举其中3种：
1. 随机选择2个节点，交换路径中的这2个节点的顺序。
2. 随机选择2个节点，将路径中这2个节点间的节点顺序逆转。
3. 随机选择3个节点m，n，k，然后将节点m与n间的节点移位到节点k后面。
## 五. 算法评价
模拟退火算法是一种随机算法，并不一定能找到全局的最优解，可以比较快的找到问题的近似最优解。 如果参数设置得当，模拟退火算法搜索效率比穷举法要高。
遗传算法（GA）的C语言实现
最后更新于：2022-04-01 16:06:34
## 问题：
在下面的程序中将要运用遗传算法对一个多项式求最小值
要求在(-8,8)间寻找使表达式达到最小的x，误差为0.001
## 问题分析：
编码：采用常规码，即二进制码编码。构造简单，交叉、变异的实现非常容易，同时解的表达也很简洁、直观。可以每0.001取一个点，这样理论误差讲小于0.0005，可以满足题目中的误差要求。此事总的求解空间为：
N = (8 - (-8)) * 1000 = 160000
可以用n = 14位二进制来表示。
群体规模m：
群体规模m可以选择 n ~ 2n 的一个确定的数，这里选择 m = 20
初始种群的选取：
在这里初始种群将在值域范围内随机选取
终止规则：
①最优解在连续的20次循环中改变量小于0.01，此事认为这个最优解为满足题目要求的最优解，求解成功，退出程序
②总的循环次数大于1200次时，循环也将结束，这种情况按照求解失败处理
交叉规则：
采用最常用的双亲双子法
选择：
在进行交叉、变异后，种群中的个体个数达到2m个，将这2m个染色体按其适应度进行排序，保留最优的m个淘汰其他的，使种群在整体上得到进化
~~~
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define SUM 20 //总共的染色体数量
#define MAXloop 1200 //最大循环次数
#define error 0.01 //若两次最优值之差小于此数则认为结果没有改变
#define crossp 0.7 //交叉概率
#define mp 0.04 //变异概率
//用于求解函数y=x^6-10x^5-26x^4+344x^3+193x^2-1846x-1680在（-8，8）之间的最小值
struct gen //定义染色体结构
{
int info; //染色体结构，用一整型数的后14位作为染色体编码
float suitability; //次染色体所对应的适应度函数值，在本题中为表达式的值
};
struct gen gen_group[SUM];//定义一个含有20个染色体的组
struct gen gen_new[SUM];
struct gen gen_result; //记录最优的染色体
int result_unchange_time; //记录在error前提下最优值为改变的循环次数
struct log //形成链表，记录每次循环所产生的最优的适应度
{
float suitability;
struct log *next;
}llog,*head,*end;
int log_num; //链表长度
/**************函数声明******************/
void initiate(); //初始化函数，主要负责产生初始化种群
void evaluation(int flag); //评估种群中各染色体的适应度，并据此进行排序
void cross(); //交叉函数
void selection(); //选择函数
int record(); //记录每次循环产生的最优解并判断是否终止循环
void mutation(); //变异函数
void showresult(int); //显示结果
//-----------------------以上函数由主函数直接调用
int randsign(float p); //按照概率p产生随机数0、1，其值为1的概率为p
int randbit(int i,int j); //产生一个在i，j两个数之间的随机整数
int randnum(); //随机产生一个由14个基因组成的染色体
int convertionD2B(float x);//对现实解空间的可能解x进行二进制编码（染色体形式）
float convertionB2D(int x); //将二进制编码x转化为现实解空间的值
int createmask(int a); //用于交叉操作
void main()
{
int i,flag;
flag=0;
initiate(); //产生初始化种群
evaluation( 0 ); //对初始化种群进行评估、排序
for( i = 0 ; i < MAXloop ; i++ )
{
cross(); //进行交叉操作
evaluation( 1 ); //对子种群进行评估、排序
selection(); //对父子种群中选择最优的NUM个作为新的父种群
if( record() == 1 ) //满足终止规则1，则flag=1并停止循环
{
flag = 1;
break;
}
mutation(); //变异操作
}
showresult( flag ); //按照flag显示寻优结果
}
void initiate()
{
int i , stime;
long ltime;
ltime=time(NULL);
stime=(unsigned)ltime/2;
srand(stime);
for( i = 0 ; i < SUM ; i++ )
{
gen_group[i].info = randnum(); //调用randnum()函数建立初始种群
}
gen_result.suitability=1000;
result_unchange_time=0;
head=end=(struct log *)malloc(sizeof(llog));//初始化链表
if(head==NULL)
{
printf("\n内存不够！\n");
exit(0);
}
end->next = NULL;
log_num = 1;
}
void evaluation(int flag)
{
int i,j;
struct gen *genp;
int gentinfo;
float gentsuitability;
float x;
if( flag == 0 ) // flag=0的时候对父种群进行操作
genp = gen_group;
else genp = gen_new;
for(i = 0 ; i < SUM ; i++)//计算各染色体对应的表达式值
{
x = convertionB2D( genp[i].info );
genp[i].suitability = x*(x*(x*(x*(x*(x-10)-26)+344)+193)-1846)-1680; //提取公因式比原式更快
}
for(i = 0 ; i < SUM - 1 ; i++)//按表达式的值进行排序，
{
for(j = i + 1 ; j < SUM ; j++)
{
if( genp[i].suitability > genp[j].suitability )
{
gentinfo = genp[i].info;
genp[i].info = genp[j].info;
genp[j].info = gentinfo;
gentsuitability = genp[i].suitability;
genp[i].suitability = genp[j].suitability;
genp[j].suitability = gentsuitability;
}
}
}
}
void cross()
{
int i , j , k ;
int mask1 , mask2;
int a[SUM];
for(i = 0 ; i < SUM ; i++) a[i] = 0;
k = 0;
for(i = 0 ; i < SUM ; i++)
{
if( a[i] == 0)
{
for( ; ; )//随机找到一组未进行过交叉的染色体与a[i]交叉
{
j = randbit(i + 1 , SUM - 1);
if( a[j] == 0) break;
}
if(randsign(crossp) == 1) //按照crossp的概率对选择的染色体进行交叉操作
{
mask1 = createmask(randbit(0 , 14)); //由ranbit选择交叉位
mask2 = ~mask1; //形成一个类似 111000 000111之类的二进制码编码
gen_new[k].info = (gen_group[i].info) & mask1 + (gen_group[j].info) & mask2;
gen_new[k+1].info=(gen_group[i].info) & mask2 + (gen_group[j].info) & mask1;
k = k + 2;
}
else //不进行交叉
{
gen_new[k].info=gen_group[i].info;
gen_new[k+1].info=gen_group[j].info;
k=k+2;
}
a[i] = a[j] = 1;
}
}
}
void selection()
{
int i , j , k;
j = 0;
i = SUM/2-1;
if(gen_group[i].suitability < gen_new[i].suitability)
{
for(j = 1 ; j < SUM / 2 ; j++)
{
if(gen_group[i+j].suitability > gen_new[i-j].suitability)
break;
}
}
else
if(gen_group[i].suitability>gen_new[i].suitability)
{
for(j=-1;j>-SUM/2;j--)
{
if(gen_group[i+j].suitability<=gen_new[i-j].suitability)
break;
}
}
for(k=j;k<SUM/2+1;k++)
{
gen_group[i+k].info = gen_new[i-k].info;
gen_group[i+k].suitability = gen_new[i-k].suitability;
}
}
int record() //记录最优解和判断是否满足条件
{
float x;
struct log *r;
r=(struct log *)malloc(sizeof(llog));
if(r==NULL)
{
printf("\n内存不够！\n");
exit(0);
}
r->next = NULL;
end->suitability = gen_group[0].suitability;
end->next = r;
end = r;
log_num++;
x = gen_result.suitability - gen_group[0].suitability;
if(x < 0)x = -x;
if(x < error)
{
result_unchange_time++;
if(result_unchange_time >= 20)return 1;
}
else
{
gen_result.info = gen_group[0].info;
gen_result.suitability = gen_group[0].suitability;
result_unchange_time=0;
}
return 0;
}
void mutation()
{
int i , j , m;
float x;
float gmp;
int gentinfo;
float gentsuitability;
gmp = 1 - pow(1 - mp , 11);//在基因变异概率为mp时整条染色体的变异概率
for(i = 0 ; i < SUM ; i++)
{
if(randsign(gmp) == 1)
{
j = randbit(0 , 14);
m = 1 << j;
gen_group[i].info = gen_group[i].info^m;
x = convertionB2D(gen_group[i].info);
gen_group[i].suitability = x*(x*(x*(x*(x*(x-10)-26)+344)+193)-1846)-1680;
}
}
for(i = 0 ; i < SUM - 1 ; i++)
{
for(j = i + 1 ; j < SUM ; j++)
{
if(gen_group[i].suitability > gen_group[j].suitability)
{
gentinfo = gen_group[i].info;
gen_group[i].info = gen_group[j].info;
gen_group[j].info = gentinfo;
gentsuitability = gen_group[i].suitability;
gen_group[i].suitability = gen_group[j].suitability;
gen_group[j].suitability = gentsuitability;
}
}
}
/*
*为了提高执行速度，在进行变异操作的时候并没有直接确定需要进行变异的位
*而是先以cmp概率确定将要发生变异的染色体，再从染色体中随进选择一个基因进行变异
*由于进行选择和变异后父代种群的次序已被打乱，因此，在变异前后对种群进行一次排序
*/
}
void showresult(int flag)//显示搜索结果并释放内存
{
int i , j;
struct log *logprint,*logfree;
FILE *logf;
if(flag == 0)
printf("已到最大搜索次数，搜索失败！");
else
{
printf("当取值%f时表达式达到最小值为%f\n",convertionB2D(gen_result.info),gen_result.suitability);
printf("收敛过程记录于文件log.txt");
if((logf = fopen("log.txt" , "w+")) == NULL)
{
printf("Cannot create/open file");
exit(1);
}
logprint=head;
for(i = 0 ; i < log_num ; i = i + 5)//对收敛过程进行显示
{
for(j = 0 ; (j < 5) & ((i + j) < log_num-1) ; j++)
{
fprintf(logf , "%20f" , logprint->suitability);
logprint=logprint->next;
}
fprintf(logf,"\n\n");
}
}
for(i = 0 ; i< log_num ; i++)//释放内存
{
logfree=head;
head=head->next;
free(logfree);
fclose(logf);
}
getchar();
}
int randsign(float p)//按概率p返回1
{
if(rand() > (p * 32768))
return 0;
else return 1;
}
int randbit(int i, int j)//产生在i与j之间的一个随机数
{
int a , l;
l = j - i + 1;
a = i + rand() * l / 32768;
return a;
}
int randnum()
{
int x;
x = rand() / 2;
return x;
}
float convertionB2D(int x)
{
float y;
y = x;
y = (y - 8192) / 1000;
return y;
}
int convertionD2B(float x)
{
int g;
g = (x * 1000) + 8192;
return g;
}
int createmask(int a)
{
int mask;
mask=(1 << (a + 1)) - 1;
return mask;
}
~~~
本例中所实现的是最基本的遗传算法，在实际应用中，往往要事先建立问题的数学模型，选择适合的适应度函数，并且根据问题的特点确定求解空间及染色体形式。本例事先用mathmatica绘出函数的曲线图，确认函数的最小值位于±8之间，然后在此区间内求解。
交叉、变异、选择的实现方法常常由选定的适应度函数以及染色体形式决定，可供选择的方式有很多。读者可以自行查阅资料
前言
最后更新于：2022-04-01 16:06:31
> 原文出处：[计算机数学模型](http://blog.csdn.net/column/details/mathematicalmodelin.html)
作者：[u013007900](http://blog.csdn.net/u013007900)
**本系列文章经作者授权在看云整理发布，未经作者允许，请勿转载！**
# 计算机数学模型
> 里面收集部分工业界常用算法与数学证明，以及部分数学建模问题的求解。