数学建模方法总结(matlab)

1. 前言

在这篇文章中，我会把我所接触的数学建模的知识的代码分享给大家，有的是我自己写的，更多的是网上借鉴并修改为可执行可用的代码

需要说明的是，我不太了解其中的数学实现的具体方法和原理，我只分享源码和作用条件以及使用的经验说明（更详细见注释），网上有不少文章对某些方法讲得非常透彻，因此我没必要赘述

各位只需要忘记代码或使用条件时从我这里 copy 就好了

2. 灰色预测模型

这是用一种在数学上对数据进行累加和累减以取巧性地削弱预测数据与原始数据的关联性

最好用于短期预测，只能用于递增序列

2.1. GM（1，1）

我也不懂为什么叫 GM(1,1)，但是用得最多的就是它，又称一阶灰色方程

clc,clear;close all
y =  [0.6 1.1 1.6]
yuceshu = 3;
% 本程序主要用来计算根据灰色理论建立的模型的预测值
% 应用的数学模型是 GM(1,1)
% 原始数据的处理方法是一次累加
% 渡辺曜改进了这个代码,通常GM预测第一个数是不够优秀的
% y 是应该输入的数组，yuceshu是你想预测的后几个数据 我们一般就预测两个
% 示例 huiseyuce([1,2,3],2)
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:nyy(i)=yy(i-1)+y(i);
end
B=ones(n-1,2);
for i=1:(n-1)B(i,1)=-(yy(i)+yy(i+1))/2;B(i,2)=1;
end
BT=B';
for j=1:n-1YN(j)=y(j+1);
end
YN=YN';
A=inv(BT*B)*BT*YN;
a=A(1);
u=A(2);
t=u/a;
i=1:n+yuceshu;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;
yys(1)=y(1);
for j=n+yuceshu:-1:2ys(j)=yys(j)-yys(j-1);
end
x=1:n;
% 把下面这里的2改成1 就可以看到GM拟合的第一个数据
xs=2:n+yuceshu;
yn=ys(2:n+yuceshu);
plot(x,y,'^r',xs,yn,'*-b');
det=0;sum1=0;
sumpe=0;
for i=1:nsumpe=sumpe+y(i);
end
pe=sumpe/n;
for i=1:nsum1=sum1+(y(i)-pe).^2;
end
s1=sqrt(sum1/n);
sumce=0;
for i=2:nsumce=sumce+(y(i)-yn(i));
end
ce=sumce/(n-1);
sum2=0;
for i=2:nsum2=sum2+(y(i)-yn(i)-ce).^2;
end
s2=sqrt(sum2/(n-1));
c=(s2)/(s1);
disp(['后验差比值为：',num2str(c)]);
if c<0.35disp('系统预测精度好')
else if c<0.5disp('系统预测精度合格')else if c<0.65disp('系统预测精度勉强')elsedisp('系统预测精度不合格')endend
enddisp(['下个拟合值为 ',num2str(ys(n+1))]);
for i=2:yuceshudisp(['再下个拟合值为',num2str(ys(n+i))]);
end

2.2. 分数阶 GM(1,1)

分数阶灰色方程的预测效果优于原始灰色模型，我只找到了 python 版本的可运行代码

首先，将 GM(1,1)打包成一个类

GM.py

##渡辺笔记
##传统的灰色模型
import numpy as np
import xlrd
import xlwt
import datetime
class Grey_model(object):def __init__(self,input_value):##初始化过程，先建立好变量空间，即累加序列，背景值序列，B和Y矩阵，拟合序列，预测值序列等。self.input_value=input_valueself.accumulation_value=np.zeros(len(input_value))self.background_value=np.zeros(len(input_value)-1)self.y_matrix_value=np.mat(np.zeros((len(input_value)-1,1)))self.b_matrix_value=np.mat(np.ones((len(input_value)-1,2)))self.accumulation()##计算累加序列def accumulation(self):for i in range(len(self.input_value)):self.accumulation_value[i]=np.sum(self.input_value[0:i+1])##计算Z值def background_values(self):for i in range(len(self.accumulation_value)-1):self.background_value[i]=(self.accumulation_value[i]+self.accumulation_value[i+1])/2##计算B矩阵def b_matrix(self):self.background_values()for i in range (self.b_matrix_value.shape[0]):self.b_matrix_value[i,0]=-self.background_value[i]##计算Y矩阵def y_matrix(self):for i in range(self.y_matrix_value.shape[0]):self.y_matrix_value[i]=self.accumulation_value[i+1]-self.accumulation_value[i]##计算参数矩阵U：U=(B^T*Y)^-1*B^T*Ydef u_matrix(self):self.y_matrix()self.b_matrix()self.u_matrix_values=(self.b_matrix_value.T*self.b_matrix_value)**-1 * self.b_matrix_value.T * self.y_matrix_value##下面把矩阵格式转化为数组格式，再转化为列表格式self.u_matrix_array=np.array(self.u_matrix_values.T)[0]self.u_matrix_value=self.u_matrix_array.tolist()##计算预测值def predict(self,number_of_forecast):self.u_matrix()self.predicted_accumulation_value=[]##使用了float(%.2f% 来调整输出值的小数的个数self.predicted_value=[float(%.2f%(self.input_value[0]))]for i in range(len(self.input_value)+int(number_of_forecast)):self.predicted_accumulation_value.append(((self.accumulation_value[0]-self.u_matrix_value[1]/self.u_matrix_value[0])*np.exp(-self.u_matrix_value[0]*(i)))+self.u_matrix_value[1]/self.u_matrix_value[0])##上式中e^-at,由于python中从0开始算，故t减去1for i in range(len(self.predicted_accumulation_value)-1):self.predicted_value.append(float(%.2f%(self.predicted_accumulation_value[i+1]-self.predicted_accumulation_value[i])))##计算误差def test_error(self):MAPE_list=[]for i in range(len(self.input_value)):MAPE_list.append(abs((self.predicted_value[i]-self.input_value[i])/self.input_value[i]))self.MAPE=str(100*(np.mean(MAPE_list[1:])))[:5]+ '%'## 打开xls，只能用这种格式了
work_book = xlrd.open_workbook('first_shuju.xls')
## 按sheet表名称获取sheet对象，名称分大小写
sheet_1 = work_book.sheet_by_name('Sheet1')
col_0_value = sheet_1.col_values(6,2,15)##输入变量，order表示分数阶的阶数
input_value = col_0_value
if __name__ == '__main__':input_value = col_0_valuenumber_of_forecast=3A=Grey_model(input_value)A.predict(number_of_forecast)A.test_error()print('预测值')print(A.predicted_value)print('\nMAPE值')print(A.MAPE)## xlwt_data = A.predicted_valuedef save_excel(target_list, output_file_name):将数据写入xls文件if not output_file_name.endswith('.xls'):output_file_name += '.xls'workbook = xlwt.Workbook(encoding='utf-8')ws = workbook.add_sheet(sheet1)rows = len(target_list)lines = len(target_list[0])for i in range(rows):for j in range(lines):ws.write(i, j, target_list[i][j])workbook.save(output_file_name)## output_file_name = 'putong_GM.xls'
## xlwt_data = tuple(xlwt_data)
## data = []
## data.append(xlwt_data)
## print(data)
## save_excel(data, output_file_name)

在这个类的基础上，运用分数阶进行优化

FGM.py

##渡辺笔记
##分数阶灰色模型
from GM import Grey_model
import numpy as np
from scipy.special import gamma
import xlrd
import xlwt##更新累加序列（式1）
def updata_fractional_accumulation(input_value,order):accumulation_value=np.zeros(len(input_value))for i in range(len(input_value)):for j in range(i+1):tmp=gamma(i-j+order-1+1)/(gamma(i-j+1)*gamma(order-1+1))accumulation_value[i]+=tmp*input_value[j]return accumulation_value##更新还原方法（参考式2和3）
def updata_predicted_results(predicted_accumulation_value,order):if order!=1:predicted_one_accumulation_value=updata_fractional_accumulation(predicted_accumulation_value,1-order)else:predicted_one_accumulation_value=predicted_accumulation_valuepredicted_value=[float(%.2f%(predicted_one_accumulation_value[0]))]for i in range(len(predicted_accumulation_value)-1):predicted_value.append(float(%.2f%(predicted_one_accumulation_value[i+1]-predicted_one_accumulation_value[i])))return predicted_value## 打开xls，只能用这种格式了
work_book = xlrd.open_workbook('first_shuju.xls')
## 按sheet表名称获取sheet对象，名称分大小写
sheet_1 = work_book.sheet_by_name('Sheet1')
col_0_value = sheet_1.col_values(6,2,19)print(col_0_value)##输入变量，order表示分数阶的阶数
input_value = col_0_value
number_of_forecast=10
order=0.1
A=Grey_model(input_value)##引入GM模型的计算模块
A.accumulation_value=updata_fractional_accumulation(input_value,order)##更新累加值
A.predict(number_of_forecast)
A.predicted_value= updata_predicted_results(A.predicted_accumulation_value,order)##更新预测值的还原方式
A.test_error()##误差用MAPE值来衡量
print('预测值')
print(A.predicted_value)
print('\nMAPE值')
print(A.MAPE)def save_excel(target_list, output_file_name):将数据写入xls文件if not output_file_name.endswith('.xls'):output_file_name += '.xls'workbook = xlwt.Workbook(encoding='utf-8')ws = workbook.add_sheet(sheet1)rows = len(target_list)lines = len(target_list[0])for i in range(rows):for j in range(lines):ws.write(i, j, target_list[i][j])workbook.save(output_file_name)xlwt_data = A.predicted_value
output_file_name = 'fenjieshu_putong_GM.xls'
xlwt_data = tuple(xlwt_data)
data = []
data.append(xlwt_data)
print(data)
save_excel(data, output_file_name)

3. 最短路径

3.1. Floyd

感觉大多都能用 floyd 解决

% 这是floyd算法，以下注释为渡边所加
% floyd不能解决多条最短路径和负权回路的问题
% 解稠密图效果最佳，边权可正可负，是双源解法
% 对于稠密图，效率要高于执行|V|次Dijkstra算法，也要高于执行|V|次SPFA算法。
% 缺点：时间复杂度比较高，不适合计算大量数据，注意
% 输入示例
% W=[0 50 inf 40 25 10; 50 0 15 20 inf 25; inf 15 0 10 20 inf; 40 20 10 0 10 25;25 inf 20 10 0 55; 10 25 inf 25 55 0];
% [D, P, dis, path] = floyd(W, 1, 4)
function [D, P, dis, path] = floyd(W, start, stop) %start为指定起始结点，stop为指定终止结点
D = W; %最短距离矩阵赋初值
n = length(D); %n为结点个数
P = zeros(n,n); %路由矩阵赋初值for i = 1:nfor j = 1:nP(i,j) = j;end
endfor k = 1:nfor i = 1:nfor j = 1:nif D(i,k) + D(k,j) < D(i,j)D(i,j) = D(i,k) + D(k,j);   %核心代码P(i,j) = P(i,k);endendend
endif nargin ~= 3errordlg('参数个数输入有误！', 'Warning!')
elsedis = D(start, stop); %指定两结点间的最短距离m(1) = start;i = 1;while P(m(i),stop) ~= stopk = i + 1;m(k) = P(m(i),stop);i = i + 1;endm(i+1) = stop;path = m; %指定两结点之间的最短路径
end

3.2. Dijkstra

Dijkstra 运算速度快，但是用不了负权图

% 这是 dijkstra 算法，以下注释为渡边所加
% 不能用于负权图，但是由于时间复杂度低，适合在大数据中运算
% 是单源解法
% 输入示例
% W=[0 50 inf 40 25 10; 50 0 15 20 inf 25; inf 15 0 10 20 inf; 40 20 10 0 10 25;25 inf 20 10 0 55; 10 25 inf 25 55 0];
% [dis, path] = dijkstra(W, 1, 4)
function [dis,path] = dijkstra( W,st,e )
%   W  权值矩阵   st 搜索的起点   e 搜索的终点
n=length(W);%节点数
D = W(st,:);
visit= ones(1:n); visit(st)=0;
parent = zeros(1,n);%记录每个节点的上一个节点path =[];for i=1:n-1temp = [];%从起点出发，找最短距离的下一个点，每次不会重复原来的轨迹，设置visit判断节点是否访问for j=1:nif visit(j)temp =[temp D(j)];elsetemp =[temp inf];endend[value,index] = min(temp);visit(index) = 0;%更新 如果经过index节点，从起点到每个节点的路径长度更小，则更新，记录前趋节点，方便后面回溯循迹for k=1:nif D(k)>D(index)+W(index,k)D(k) = D(index)+W(index,k);parent(k) = index;endendenddistance = D(e);%最短距离
%回溯法  从尾部往前寻找搜索路径
t = e;
while t~=st && t>0path =[t,path];p=parent(t);t=p;
end
path =[st,path];%最短路径
dis = distance;
end

3.3. matlab 自带最短路径

clc,clear,close all;
% 程序基于但不基于 Dijkstra，解决了负权问题
s = [1 1 1 2 2 3 3 4 5 5 6 7];
t = [2 4 8 3 7 4 6 5 6 8 7 8];
weights = [10 10 1 10 1 -10 1 1 12 12 12 12];
% names = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H'};
% G = digraph(s,t,weights,names);
G = digraph(s,t,weights);
p = plot(G,'Layout','force','EdgeLabel',G.Edges.Weight);[path1,d] = shortestpath(G,6,8);
highlight(p,path1,'EdgeColor','r')

4. 层次分析法

感觉这个方法就像我们土木工程所谓的“拍脑袋”，先拍脑袋得出一个比较矩阵，待反复测试符合规范后，就可以拿去比较了

% 本代码只需要在提示输入处，输入构成的比较矩阵 就会输出各指标的权重值。
% 例子： 选拔干部考虑5个条件：品德x1，才能x2，资历x3，年龄x4，群众关系x5。某决策人用成对比较法，得到成对比较阵如下：
% [1,2,7,5,5;
% 1/2,1,4,3,3;
% 1/7,1/4,1,1/2,1/3;
% 1/5,1/3,2,1,1;
% 1/5,1/3,3,1,1]
% 其中的x1=1;x2=2;x3=7;x4=5;x5=5; 怎么来的呢？
% 其实这些x的值就是由决策人决定了，对应值也就是决策人认为的重要性了
clc;
clear;
% 判断矩阵A
A= [1 3 3 31/3 1 3 51/5 1/3 1 31/5 1/5 1/3 1];
[m,n]=size(A);                     %获取指标个数
RI=[0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.51];
R=rank(A);                         %求判断矩阵的秩
[V,D]=eig(A);                      %求判断矩阵的特征值和特征向量，V特征值，D特征向量；
tz=max(D);
B=max(tz);                         %最大特征值
[row, col]=find(D==B);             %最大特征值所在位置
C=V(:,col);                        %对应特征向量
CI=(B-n)/(n-1);                    %计算一致性检验指标CI
CR=CI/RI(1,n);
if CR<0.10disp('CI=');disp(CI);disp('CR=');disp(CR);disp('对比矩阵A通过一致性检验，各向量权重向量Q为：');Q=zeros(n,1);for i=1:nQ(i,1)=C(i,1)/sum(C(:,1)); %特征向量标准化endQ'                              %输出权重向量
elsedisp('对比矩阵A未通过一致性检验，需对对比矩阵A重新构造');
end
sc = Q';

5. 聚类分析

这类方法对于我的专业毫无实际用处，我并不需要知道管道应该分成几类，因为规范早就规定好了的

聚类分析方法丰富，要根据实际情况谨慎选取

5.1. k-mean

% k_mean 算法，渡边笔记
% 对处理大型数据集，该算法保持可收缩性，高效性；
% 当簇接近正态分布时，效果最好；
% 若簇中含有异常点，将导致均值偏离严重；（即对噪声、孤立点敏感）；
% 不适用于发现非凸形状的簇，或大小差别大的簇
% 中k是事先给定的，这个k值难以估计，很多时候并不知道数据集应该分成多少个类别最合适
clear; clc;data = [11 2 02 2 24 3 3]; % 输入数据
k = 3; % 分类数X = mapminmax(data',0,1)'; % 按列最小最大规范化到[0,1]
T = kmeans(X,k); % 直接调用kmeans函数
for i = 1:kkong = find(T == i);fprintf('第 %d 类 :',i);disp(data(kong)')
end

5.2. 层次聚类

% 层次聚类，渡边笔记
% 看图说话，分类清晰明了
data =[17.0    0.31    0.42    0.74    1.16    1.81    3.79
18.0    0.1 0.32    0.31    0.69    0.76    2.68
19.0    0.18    0.35    0.47    0.77    1.57    3.93
20.0    0.23    0.15    0.5 1.04    0.99    4.11
21.0    0.13    0.15    0.2 1.42    1.53    2.51]';X=mapminmax(data,0,1)'; % 按行最小最大规范化到[0,1]
Y = pdist(X); % 计算矩阵X中样本两两之间的距离，但得到的Y是个行向量
D = squareform(Y); % 将行向量的Y转换成方阵，方便观察两点距离(方阵的对角线元素都是0)
Z = linkage(Y); % 产生层次聚类树，默认采用最近距离作为类间距离的计算公式
dendrogram(Z); % 图示层次聚类树
title('层次聚类法');ylabel('标准距离');

5.3. 直接聚类

% 直接聚类，渡笔记
% 不管机理，不管距离，比较愚蠢但是有用且方便
clear; clc;
data =[36.6    4.46    36.72   6.37    32.39   4.63    15.43
36.80   12.46   47.21   18.66   9.22    1.69    10.76
67.88   2.76    39.1    4.2 36.92   1.87    15.15
48.9    2.85    36.85   7.23    38.29   3.51    11.27
60.5    2.23    27.25   6.8 45  8.73    9.99
];
X=mapminmax(data,0,1); % 按列最小最大规范化到[0,1]T1=clusterdata(X,0.2); % 如果0<cutoff<2，则当不一致系数大于cutoff时，分到不同类(簇)中
T2=clusterdata(X,3); % 如果cutoff是一个≥2的整数，则形成的不同类别数为cutoff
k1 = max(T1);
k2 = max(T2);
for i = 1:k1kong = find(T1 == i);fprintf('第 %d 类 :',i);disp(data(kong)')
end
disp('--------------------------------');
% for i = 1:k2
%     kong = find(T2 == i);
%     fprintf('第 %d 类 :',i);
%     disp(data(kong)')
% end

6. 最小生成树

好像是个示例，太久远忘了干什么用的了

zhuixiao_shengchengshu.m

clc,clear,close all;
% 使用 Kruskal 的算法来计算 无向图的最小生成树
s = [1 1 1 2 5 3 6 4 7 8 8 8];
t = [2 3 4 5 3 6 4 7 2 6 7 5];% 使用加权边创建并绘制一个立方体图
weights = [100 10 10 10 10 20 10 30 50 10 70 10];
G = graph(s,t,weights);
% figure;
p = plot(G,'EdgeLabel',G.Edges.Weight);% 计算并在图上方绘制图的最小生成树。T 包含的节点与 G 相同，但包含的边仅为后者的子集
% figure;
[T,pred] = minspantree(G);
% p = plot(G,'EdgeLabel',G.Edges.Weight);
highlight(p,T) % 高亮显示

7. 误差检验

7.1. 一般检验

看着用吧，我觉得没什么意义，反正计算机计算出的误差都挺小的，最多能在论文中凑字数

% 渡辺笔记
% 各类检验，除决定系数是1最好，都是0最好
YReal = [1 2 3 4 5];
YPred = [1 2 3 4 5.1];
% 平均绝对百分比误差（MAPE）
mape = mean(abs((YReal - YPred)./YReal));
disp(mape);
% 均方根误差（RMSE）
rmse = sqrt(mean((YPred-YReal).^2));
disp(rmse);
% 残差平方和（SSE）
sse = sum((YReal - YPred).^2);
disp(sse);
% 均方误差（MSE）
mse = mean(sum((YReal - YPred).^2));
disp(mse);
% 平均绝对误差（MAE）
mae = mean(abs(YReal - YPred));
disp(mae);
% 决定系数（R2-R-Square）
r2 = 1 - (sum((YPred - YReal).^2) / sum((YReal - mean(YReal)).^2));
disp(r2)

7.2. T 检验

% 渡边笔记 t检验
% 待检验的数据应该近服从正态分布
clc,clear;close all;
x = [1 2 2 3 3 3 3 4 4];
y = [1 2 2 3 3 3 4 4 4];
ALPHA = 0.05;
[H,P,CI,STATS]=ttest2(x,y, ALPHA);
% 其中，x，y均为行向量（维度必须相同），各表示一组数据
% ALPHA为可选参数，表示设置一个值作为t检验执行的显著性水平
% 在不设置ALPHA的情况下默认ALPHA为0.05，即计算x和y在5%的显著性水平下是否来自同一分布(假设是否被接受）
% H = 0，P > 0.05，表明零假设被接受，即x，y在统计上可看做来自同一分布的数据
% H = 1，P < 0.05，表明零假设被拒绝，即x，y在统计上认为是来自不同分布的数据，即有区分度
disp(['H = ',num2str(H)]);
disp(['P = ',num2str(P)]);
if H == 0disp('x，y在统计上可看做来自同一分布的数据')
elsedisp('x，y在统计上可看做来自不同分布的数据')
end

7.3. 回归检验

clc,clear;close all;
% 回归检验；渡边笔记
x = [0 1 2 3 4];y = [1 1 4 10 16];x = x';
y = y';
% x,y 需要转置才可以使用X = [ones(size(y)) x]; % 定义回归方程[B,BINT,R,RINT,STATS] = regress(y,X); % 回归
% [B,BINT,R,RINT,STATS] = REGRESS(Y,X,alpha)
% B 为回归系数向量
% BINT 回归系数的估计区间
% R 残差
% RINT 置信区间
% STATS 用于检验回归模型的统计量。有4个数值：判定系数r2，F统计量观测值，检验的p的值，误差方差的估计
% r2越接近1，回归方程越显著；F>F1−α(k,n−k−1)时拒绝H0，F越大，回归方程越显著；p<α时拒绝H0
% ALPHA 显著性水平（缺少时默认0.05）
y2 = B(2).* x + B(1);   %预测值
plot(x', y' ,'+', x, y2);     %回归效果图
disp('回归分析统计量为');
disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 这么多原始数据 解释，最优为 1
disp([ 'F : ',num2str(STATS(2))]);
disp([ '检验 p : ',num2str(STATS(3))]); % 日常小于 0.05 最好
disp([ '误差方差 : ',num2str(STATS(4))]); % 越小越好

7.4. 交叉检验

% 交叉检验，渡边笔记
% 提高数据的利用率
% 有效的避免过学习以及欠学习状态的发生,最后得到的结果也比较具有说服性
% 常用的 K 值有 3,6,10 等
clear,clc;close all;% 导入数据
data = reshape(1:144,12,12);
% 每行的是这个样本属性的数据，最后一个数据是标签
[data_r, data_c] = size(data);
% 将数据样本随机分割为 n 部分，做 n 折交叉检验
n = 5;
indices = crossvalind('Kfold', data_r, n);
for i = 1 : n% 获取第i份测试数据的索引逻辑值test = (indices == i);% 取反，获取第i份训练数据的索引逻辑值train = ~test;% 测试用的数据test_data = data(test, 1 : data_c - 1);test_label = data(test, data_c);% 训练用的数据train_data = data(train, 1 : data_c - 1);train_label = data(train, data_c);% 下面放入检验数据的代码end

8. 差分

% 数据差分，渡边笔记
% 差分就是后一个减去前一个，使得数据平稳化
% n 阶就是做 n 次差分，具体直接看例子就好
clc,clear;close all;
Y = [0 5 15 30 50 75 105];
disp(['未差分时为 : ',num2str(Y)])
Y_1= diff(Y, 1);
disp(['一阶差分结果为 : ',num2str(Y_1)])
Y_2 = diff(Y, 2);
disp(['二阶差分结果为 : ',num2str(Y_2)])
Y_3 = diff(Y, 3);
disp(['三阶差分结果为 : ',num2str(Y_3)])

9. 立法白噪音

% 立法数白噪声检验，渡边笔记
% 如果易知原始数据不是平稳的，则不能做随机性检验
% 接下来要求差分,目的: 变成平稳的数据
% p 如果比 0.05 小就不是白噪声序列，可以使用时间序列
% 某种现象的指标数值按照时间顺序排列而成的数值序列
% 因为时间序列是某个指标数值长期变化的数值表现
% 所以时间序列数值变化背后必然蕴含着数值变换的规律性
% 这些规律性就是时间序列分析的切入点
% 如果原数据平稳，但是没通过，可以直接差分
clc,clear;close all;x = []; % 时间，一般做题就是顺序时间排列
yanchi=[6 12 18];   % 通常是做6 12 18 24步延迟，这个数据的选择上限请根据报错来调整
y = [970279
1259308
1127571
1163959
1169540
1076938
991350
953275
951508
904434
889381
864015
836236
]';;
[h1] = adftest(y);    %检验是否平稳
if h1 == 1disp('数据是平稳的');y_1 = y;
elsedisp('数据是不平稳的');i = 1;while 1y_1 = diff(y,i);   % 在这里对数据进行 n 阶差分[h1,p1,adf,ljz] = adftest(y_1);    %检验是否平稳if h1 == 1disp(['差分后是平稳的，做了 ',num2str(i),' 阶差分']);subplot(2,2,3)plot(y_1);  % 一阶差分数据时序图title('一阶差分数据时序图')subplot(2,2,4)autocorr(y_1); % 一阶差分数据的自相关系数图title('一阶差分数据自相关系数图');breakendi = i + 1;end
end
% 随时间的变化值
subplot(2,2,1)
plot(y); % 原始数据时序图
title('原始数据时序图')
subplot(2,2,2)
autocorr(y); % 原始数据的自相关系数图
title('原始数据自相关系图像')[H,pValue,Qstat,CriticalValue]=lbqtest(y_1,'lags',yanchi); % H.结果，pValue.p值, Qstat.卡方统计量
fprintf('%15s%15s%15s','延迟阶数','卡方统计量','p值');
fprintf('\n');
for i=1:length(yanchi)    % i = 1,时候为6，i = 2时候为12fprintf('%18f%19f%19f',yanchi(i),Qstat(i),pValue(i));fprintf('\n');
endif sum(find(pValue > 0.05))disp('但是没通过立法白噪声检验');i = 1;while 1y_1 = diff(y,i);   % 在这里对数据进行 n 阶差分[h1,p1,adf,ljz] = adftest(y_1);    %检验是否平稳if h1 == 1disp(['差分后是平稳的，做了 ',num2str(i),' 阶差分']);subplot(2,2,3)plot(y_1);  % 一阶差分数据时序图title('一阶差分数据时序图')subplot(2,2,4)autocorr(y_1); % 一阶差分数据的自相关系数图title('一阶差分数据自相关系数图');breakendi = i + 1;end
end% 再来一次[H,pValue,Qstat,CriticalValue]=lbqtest(y_1,'lags',yanchi); % H.结果，pValue.p值, Qstat.卡方统计量
fprintf('%15s%15s%15s','延迟阶数','卡方统计量','p值');
fprintf('\n');
for i=1:length(yanchi)    % i = 1,时候为6，i = 2时候为12fprintf('%18f%19f%19f',yanchi(i),Qstat(i),pValue(i));fprintf('\n');
end
if sum(find(pValue < 0.05))disp('数据通过立法白噪声检验');
elsedisp('数据没通过立法白噪声检验');
end

10. 熵权法

% 渡边笔记 熵权法
clc;clear;close all;
% 熵权法的思想是：变量数值变化越大，变异程度越大，则其权重应该更大；反之权重则越小x = [
2.41    52.59   0   9.78    1.17
1.42    53.21   0   6.31    1.63
4.71    35.16   1   9.17    3.02
14.69   15.16   2.13    10.35   7.97
0.94    72.99   0   7.39    0.61
1.43    72.62   0   8.16    0.51
2.21    67.5    0   9.84    0.85
3.79    51.21   0   12.95   1.43
1.23    85.09   3.97    4.08    0.13
1.71    82.07   2.88    4.97    0.33
3.63    66.9    3.18    8.57    0.71
5.72    49.77   3.44    10.52   1.83
1.49    79.51   6.53    2.58    0.27
1.66    81.44   5.18    2.74    0.36
2.41    76.32   5.88    4.13    0.54
4.42    59.65   7.64    8.38    1.02
3.27    88.42   3.36    2.85    0.14
11.27   70.05   5.77    6.07    0.19
13.18   62.45   5.66    7.85    0.74
15.83   56.28   2.92    9.97    1.14
11.59   80.23   1.04    3.64    0.2
26.67   55.7    2.02    8.13    0.38
28.51   51.07   2.12    9.66    1.46
3.69    87.26   0   3.12    0.18
3.27    84.43   0   5.43    0.31
3.98    79.99   0   6.62    0.57
1.59    86.5    0   6.14    0.14
4.31    82.26   0   4.71    0.2
4.6 72.79   0   8.27    0.52
4.99    81.93   0   7.52    0.16
4.66    75.09   0   10.24   0.33
5.08    61.02   1.57    15.7    0.53
12.49   83.06   0   1.2 1.06
4.67    92.77   0   0.33    0.58
5.8 90.32   0   0.91    0.8
97.76   0   0   0   2.14
94.75   0   0   1.42    2.83
93.76   0   0   1.18    3.24
3.48    81.43   7.45    1.33    0.14
4.2 80  5.3 2.21    0.18
8.83    71.28   5.34    2.9 0.43
5.39    79.6    6.87    2.64    0.31
7.67    74.74   5.91    3.4 0.66
19.65   55.4    4.87    6.14    1.2
2.63    90.74   3.18    1.42    0.14
2.8 89.7    2.85    1.96    0.14
4.07    85.12   3.43    3.52    0.25
5.7 83.4    0   4.48    0.1
4.03    81.35   0   6.18    0.19
4.11    73.45   0   9.71    0.45
2.78    89.53   0   4.23    0.2
3.92    83.2    0   7.59    0.32
5.21    71.37   3.09    10.29   0.72
18.98   76.81   0   1.05    0.31
19.79   73.56   0   0.88    0.42
19.86   70.07   0   1.72    0.74
16.61   67.57   3.77    3.15    1.16
6.91    82.18   4.19    0   0.1
2.93    83.06   1.93    5.14    0.32
8.47    78.11   4.04    4.02    0.31
12.29   70.48   3.89    4.32    0.69
3.98    84.81   4.76    1.97    0.18
7.67    78.13   4.22    4.57    0.35
14.04   66.89   4.41    6.27    0.47
14.62   59.29   5.28    8.35    0.77
1.97    85.16   4.87    3.27    0.23
2.16    86.83   3.82    2.25    0.15
4.81    74.9    5.05    5.97    0.5
7.44    57.98   6.75    10.73   1.04
2.04    86.01   4.79    2.95    0.13
3.49    79.79   5.67    4.28    0.15
6.47    68.02   6.71    5.74    0.2
7.94    59.12   7.14    5.93    1.42
];[m,n]=size(x);
lamda=ones(1,n); % 人为修权，1代表不修改计算后的指标权重
for i=1:nx(:,i)=(x(:,i)-min(x(:,i)))/(max(x(:,i))-min(x(:,i)))+1; % 对原始数据进行非负数化、归一化处理，值介于1-2之间
end
for i=1:mfor j=1:np(i,j)=x(i,j)/sum(x(:,j));end
end
k=1/log(m);
for i=1:mfor j=1:nif p(i,j)~=0e(i,j)=p(i,j)*log(p(i,j));elsee(i,j)=0;endend
end
for j=1:nE(j)=-k*sum(e(:,j));
end
d=1-E;
for j=1:nw(j)=d(j)/sum(d);% 指标权重计算
end
for j=1:nw(j)=w(j)*lamda(j)/sum(w.*lamda);% 修改指标权重
end
for i=1:mscore(i,1)=sum(x(i,:).*w);% 计算综合分数
end
disp('各指标权重为：')
disp(w)
disp('各项综合分数为：')
disp(score)
Out = mean (score,1)

11. 数据修复

% 判断缺失值和异常值并修复，顺便光滑噪音，渡边笔记
clc,clear;close all;
x = 0:0.06:10;
y = sin(x)+0.2*rand(size(x));
y(22:34) = NaN;
y(89:95) = 50;
testdata = [x' y'];subplot(2,2,1);
plot(testdata(:,1),testdata(:,2));
title('原始数据');%% 判断数据中是否存在缺失值
if sum(isnan(testdata(:)))disp('数据存在缺失值');
elsedisp('数据不存在缺失值');
end%% 判断数据中是否存在异常值
% 1.mean 三倍标准差法 2.median 离群值法 3.quartiles 非正态的离群值法
% 4.grubbs 正态的离群值法 5.gesd 多离群值相互掩盖的离群值法
choice_1 = 5;
yichangzhi_fa = char('mean', 'median', 'quartiles', 'grubbs','gesd');
yi_chang = isoutlier(y,strtrim(yichangzhi_fa(choice_1,:)));
if sum(yi_chang)disp('数据存在异常值');
elsedisp('数据不存在异常值');
end%% 对异常值赋空值
F = find(yi_chang == 1);
y(F) = NaN; % 令数据点缺失
testdata = [x' y'];subplot(2,2,2);
plot(testdata(:,1),testdata(:,2));
title('去除差异值');%% 对数据进行补全
% 数据补全方法选择
% 1.线性插值 linear 2.分段三次样条插值 spline 3.保形分段三次样条插值 pchip
% 4.移动滑窗插补 movmean
chazhi_fa = char('linear', 'spline', 'pchip', 'movmean');
choice_2 = 2;
if choice_2 ~= 4testdata_1 = fillmissing(testdata,strtrim(chazhi_fa(choice_2,:))); % strtrim 是为了去除字符串组的空格
elsetestdata_1 = fillmissing(testdata,'movmean',10); % 窗口长度为 10 的移动均值
endsubplot(2,2,3);
plot(testdata_1(:,1),testdata_1(:,2));
title('数据补全结果');%% 进行数据平滑处理
% 滤波器选择 1.Savitzky-golay 2.rlowess 3.rloess
choice_3 = 2;
lvboqi = char('Savitzky-golay', 'rlowess', 'pchip', 'rloess');
% 通过求 n 元素移动窗口的中位数，来对数据进行平滑处理
windows = 8;
testdata_2 = smoothdata(testdata_1(:,2),strtrim(lvboqi(choice_3,:)),windows) ;subplot(2,2,4);
plot(x,testdata_2)
title('数据平滑结果');

12. 相关性分析

clc,clear,close all
% 相关性分析，渡边笔记
% 列为指标，行为数据data = [1 22 93 4 ];
% Pearson相关系数
r1 = corr(data,'type','Pearson');
disp(r1);
% Kendall相关系数
r2 = corr(data,'type','Kendall');
disp(r2);
% Spearman相关系数
r3 = corr(data,'type','Spearman');
disp(r3);

13. 马氏链

13.1. 一般马氏链

% 马氏链模型，渡边笔记
% 系统在每个时期所处的状态是随机的
% 从一时期到下时期的状态按一定概率转移
% 下时期状态只取决于本时期状态和转移概率。即已知现在，将来与过去无关（无后效性）
% 过去不影响未来的判断，是马氏链的本质
% P_0 是初始分布，P_n 是转移矩阵，则在 n 未来的概率为
% P_0 * ( P_n ^ n )
clc,clear;close all;% 第一次买盐，概率为 P_0
% 已知 买盐倾向 P_n
% P_0 = [0.2 0.4 0.4];
% P_n = [0.8 0.1 0.1
%     0.5 0.1 0.4
%     0.5 0.3 0.2];
% P = P_0 * P_n ^ 3;
% disp(P)% 一般马氏链模型
a=[160 120 120180 90 30180 30 90]; % 输入统计矩阵
P_0 = [0.4 0.3 0.4]; % 初始概率% 求状态转移矩阵
[m,~] = size(a);
ni=sum(a,2);  %计算矩阵f的行和
for i = 1:mfor j = 1:mZT(i,j) = a(i,j)./ni(i);   %求状态转移的频率end
end
disp('转移矩阵为 : ');
disp(ZT)
fprintf('%c', 8); % 删掉换行符
% 求 n 期概率
n = 2;
P = P_0 * ( ZT ^n );
disp(['第 ',num2str(n),' 期概率为 ',num2str(P)]);% 概率平衡时，有 P = P * P_n 恒成立
X = ZT';
X(1,:) = [];
for i = 2:mX(i-1,i) = X(i-1,i) - 1;
end
X(m,:) = ones(1,m);
Y = zeros(m-1,1);
Y(m,:) = ones(1);
x = X\Y; x = x';
disp(['n 趋向于无穷的 X 的解为 : ',num2str(x)]);

13.2. 带利润的马氏链

% 带利润的马氏链模型，渡边笔记
clc,clear;close all;a=[5 54 6 ]; % 输入统计矩阵
R = [9 33 -7];  % 输入利润矩阵% 对统计矩阵，求状态转移矩阵
[m,~] = size(a);
ni=sum(a,2);  %计算矩阵f的行和
for i = 1:mfor j = 1:mZT(i,j) = a(i,j)./ni(i);   %求状态转移的频率end
end
disp('转移矩阵为 : ');
disp(ZT)
fprintf('%c', 8); % 删掉换行符% 求 n 期利润
n = 3;% n = 1 时
for i = 1:mfor j = 1:mv(i,j) = R(i,j) * ZT(i,j);endSum = sum(v,2);
end
% n = n 时for i = 1:mfor j = 1:mV(i,j) = Sum(j).^(n-1) * ZT(i,j);endV_Sum = sum(V,2);V_Sum = V_Sum + Sum;if n > 1 % 边界条件disp(['处于状态 ',num2str(i),' 的 ',num2str(n),' 期利润为 ',num2str(V_Sum(i,:))])elsedisp(['处于状态 ',num2str(i),' 的 ',num2str(n),' 期利润为 ',num2str(Sum(i,:))])end
end% 带利润的马氏链一般不存在平衡
% 因为钱是赚不完的

14. 蒙特卡罗

% 渡边笔记
% Monte_Carlo 是一种方法的概述，不是什么特定的算法
% 主要思想是利用大量随机数来实现真实值的拟合
%% 例1 计算 0-3 上 x^2 的积分，正解应该是 9
N=500; x=unifrnd(0,3,N,1); % 设定随机数及其终止上限
M = mean(3*x.^2);disp(['M = ',num2str(M)]); % 将随机数代入原式求平均值，特别的，蒙特卡洛需要式乘以积分区间,下同
%% 例2 求相交面积
%给定曲线y =2 – x2 和曲线y3 = x2，曲线的交点为：P1(–1，1)、P2(1，1),函数图像如下
x_1 = -1.5:.1:1.5; y_1 = 2 - x_1.^2;
x_2 = -2:.1:2; y_2 = x_2.^2.^(1/3);
subplot(1,2,1); plot(x_1,y_1,'r',x_2,y_2,'b') %曲线围成平面有限区域，用蒙特卡罗方法计算区域面积。
P=rand(1000,2); %随机产生1000个点
%定义x y 的范围
x=2*P(:,1)-1; % 这一步是巧妙设置在(-1,2)
y=2*P(:,2);
II=find(y<=2-x.^2&y.^3>=x.^2); %找出在函数范围的数
%计算索引的长度
M=length(II);
%计算面积
S = 4*M/1000;disp(['S = ',num2str(S)]);
subplot(1,2,2); plot(x(II),y(II),'g.');
close all % 删除这个就可以显示图像
%% 例3 求圆周率
n = 1000;  %总的实验次数
m = 0;  %落在圆中点的次数
%循环实验
for i = 1:nx = 2 * rand - 1;y = 2 * rand - 1;if (x^2 + y^2 <= 1)m = m + 1;end
end
PI=4*m/n; disp(['PI = ',num2str(PI)]);

15. 模糊评价

先打包成一个 function

fuzzy_zhpj.m

function[B]=fuzzy_zhpj(model,A,R) %模糊综合评判
B=[];
[m,s1]=size(A);
[s2,n]=size(R);
if(s1~=s2)disp('A的列不等于R的行');
elseif(model==1)                 %主因素决定型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=0;if(A(i,k)<R(k,j))x=A(i,k);elsex=R(k,j);endif(B(i,j)<x)B(i,j)=x;endendendendelseif(model==2)               %主因素突出型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=A(i,k)*R(k,j);if(B(i,j)<x)B(i,j)=x;endendendendelseif(model==3)              %加权平均型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)B(i,j)=B(i,j)+A(i,k)*R(k,j);endendendelseif(model==4)             %取小上界和型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=0;x=min(A(i,k),R(k,j));B(i,j)=B(i,j)+x;endB(i,j)=min(B(i,j),1);endendelseif(model==5)            %均衡平均型C=[];C=sum(R);for(j=1:n)for(i=1:s2)R(i,j)=R(i,j)/C(j);endendfor(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=0;x=min(A(i,k),R(k,j));B(i,j)=B(i,j)+x;endendendelsedisp('模型赋值不当');end
end
end

使用方法

clc,clear;
A1=[0.1 0.9];
% A1=[0.4 0.35 0.15 0.1];
R=[
1   4
2 2];sco = fuzzy_zhpj(1,A1,R)
% fuzzy_zhpj(2,A1,R)
% fuzzy_zhpj(3,A1,R)
% fuzzy_zhpj(1,A2,R)

16. 模拟退火

16.1. 运用数学方法

clc,clear
% 模拟退火(求最小)
temperature = 100; % 初始温度
final_temperature = 1; % 最低温度
time_temperature = 1; % 用于计算步数
time_tuihuo = 10; % 退火步长
cooling_rate = 0.95; % 降温步数
% 初始化随机数生成器，以使结果具备可重复性
% rng(0,'twister');
% 生成范围a-b的随机数
a = -10; b = 10;
% rang_math = (b-a).*rand(1000,1) + a;
rang_math = (b-a).*rand + a;
f = @(x)(...x.^2 ...    % 定义函数，求最小值);
current_old = f(rang_math);while final_temperature < temperaturerang_math = (b-a).*rand + a; % 随机数current_new =  f(rang_math); % 生成新解diff = current_new - current_old;% Metropolis 准则if(diff<0)||(rand < exp(-diff/(temperature))) % 接受新解的条件route = rang_math;current_old = current_new;time_temperature = time_temperature + 1;endif time_temperature == time_tuihuotemperature = temperature * cooling_rate;time_temperature = 0;endendplot([a:b],f([a:b])); % 原函数图像
hold on;
plot(route,current_old,'ko', 'linewidth', 2); % 作解的图像

16.2. 运用 matlab 自带工具

######## 2.1. 一元示例

% 渡边自己手打的
% 模拟退火解一元函数
clc,clear;close all;% 设置区间
x = -10:1:10;f = @(x)(...
x.^4 ...    % 定义函数，求最小值
);% f = @(x)f_xy(x(1),x(2));
x0 = rand(1,1);   % 退火开始点
lb = [];
ub = []; % 退火实施范围，可以不设置
% 实时观测
options = saoptimset('MaxIter',20,... % 迭代次数
'StallIterLim',300,... % 最高温度
'TolFun',1e-100,... % 最低温度
'AnnealingFcn',@annealingfast,'InitialTemperature',...
100,'TemperatureFcn',@temperatureexp,'ReannealInterval',500,'PlotFcns',...
{@saplotbestx, @saplotbestf, @saplotx, @saplotf,@saplottemperature});[jie,fval] = simulannealbnd(f,x0,lb,ub,options);
fprintf('最优解为 x = %.10f, y = %.10f\n', jie, fval );
figure;
fun = x.^2 + 2*x; % 再定义一次
plot(x,fun);
hold on;
plot(jie,fval,'ko', 'linewidth', 2);

######## 2.2. 二元示例

% 渡边自己手打的
% 模拟退火解二元函数
clc,clear;close all;% 设置区间
xx = -10:1:10;
yy = -10:1:10;
[x, y]=meshgrid(xx, yy);f_xy = @(x,y)(...
x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20 ...    % 定义函数，求最小值
);
fun = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20; % 再定义一次
f = @(x)f_xy(x(1),x(2));
x0 = rand(1,2);   % 退火开始点
lb = [];
ub = []; % 退火实施范围，可以不设置
% 实时观测
options = saoptimset('MaxIter',20,... % 迭代次数
'StallIterLim',300,... % 最高温度
'TolFun',1e-100,... % 最低温度
'AnnealingFcn',@annealingfast,'InitialTemperature',...
100,'TemperatureFcn',@temperatureexp,'ReannealInterval',500,'PlotFcns',...
{@saplotbestx, @saplotbestf, @saplotx, @saplotf,@saplottemperature});[jie,fval] = simulannealbnd(f,x0,lb,ub,options);
fprintf('最优解为 x = %.10f, y = %.10f\n', jie(1),jie(2) );
fprintf('最优值为 z = %.10f\n',fval);
figure;
surf(x,y,fun);
hold on;
plot3(jie(1),jie(2),fval,'ko', 'linewidth', 3);

17. 神经网络

17.1. BP

% 渡辺笔记，BP 神经网络
clc,clear,close all;
%输入数据矩阵,每行为对应的属性
p = [1 1 12 2 23 3 3];
%目标（输出）数据矩阵
t = [6 6 6];%对训练集中的输入数据矩阵和目标数据矩阵进行归一化处理
[pn, inputStr] = mapminmax(p);
[tn, outputStr] = mapminmax(t);%建立BP神经网络
net = newff(pn, tn, [3 7 2], {'purelin', 'logsig', 'purelin'});
%每10轮回显示一次结果
net.trainParam.show = 10;
%最大训练次数
net.trainParam.epochs = 5000;
%网络的学习速率
net.trainParam.lr = 0.05;
%训练网络所要达到的目标误差
net.trainParam.goal = 0.65 * 10^(-3);
%网络误差如果连续6次迭代都没变化，则matlab会默认终止训练。为了让程序继续运行，用以下命令取消这条设置
net.divideFcn = '';
%开始训练网络
net = train(net, pn, tn);
%使用训练好的网络，基于训练集的数据对BP网络进行仿真得到网络输出结果
%(因为输入样本（训练集）容量较少，否则一般必须用新鲜数据进行仿真测试)
answer = sim(net, pn);
%反归一化
answer1 = mapminmax('reverse', answer, outputStr);

17.2. 灰色

clc,clear,close all;
% 渡边笔记 灰色神经网络filename_1 = 'first_shuju.xlsx';
sheet_1 = '';
range_3 = 'G3:G18';
for j = 1:16years(j) = j+3;
end
[num_3, ~,~] = xlsread(filename_1,sheet_1,range_3);
[m,~] = size(num_3);
IN=1:m;
for i = 1:msr(i) = num_3(i);
end
OUT=sr;
[X,minx,maxx,T,mint,maxt]=premnmx(IN,OUT);
q=50;
q1=0;
q0=70;
while(q1<q)q=q0;[M,N]=size(X);[L,N]=size(T);net=newff(minmax(X),[q,L],{'tansig','purelin'},'trainlm');net.trainParam.lr=0.05;net.trainParam.epochs=10000;net.trainParam.goal=1e-6;[net,tr]=train(net,X,T);Y=sim(net,X);Y=postmnmx(Y,mint,maxt);%灰色关联分析，调整网络隐层节点p=0.5;e=0.5;%此两个系数的设定是根据一些论文，已经实验的尝试得出的an=repmat(net.b{1},1,N);op=tansig(net.iw{1,1}*X+an);op1=op';T0=T';T1=repmat(T0,1,q);DIF=abs(T1-op1);MIN=min(min(DIF));MAX=max(max(DIF));Si=(MIN+p*MAX)./(DIF+p*MAX);ri=sum(Si)/N;D=find(ri>=e);[q0,q1]=size(D);q0=q1;
end
q0;
ri;
D;
q=q1;%进行测试
PRD=1:m;
PRD=PRD';
P=tramnmx(PRD,minx,maxx);
TNEW=sim(net,P');
TNEW=postmnmx(TNEW,mint,maxt);YY=OUT;
YC=TNEW;
figure
plot(years,YY,'b+-',years,YC,'r')
legend('灰色神经网络','生态用水量');
xlabel('年份');
RES0=YC-YY;
res0=RES0./YY;
figure
bar(years,res0);
xlabel('训练次数');
ylabel('训练误差');

17.3. 拟合

拟合神经网络是什么我忘记了，反正能用就行了

shenjin.m

function [net,tr] = shenjin(choice,P,T)
net=newff(minmax(P),[20,1],... % 隐单元的神经元数 20，输出 1，可调
{'tansig','purelin'});
if(choice==1)%  采用 L-M 优化算法 TRAINLMnet.trainFcn='trainlm';%  设置训练参数net.trainParam.epochs = 500;net.trainParam.goal = 1e-6;net=init(net);%  重新初始化
elseif(choice==2)%  采用贝叶斯正则化算法 TRAINBRnet.trainFcn='trainbr';%  设置训练参数net.trainParam.epochs = 500;randn('seed',192736547);net = init(net);%  重新初始化
elseif(choice==3)%  采用动量梯度下降算法 TRAINGDM%  当前输入层权值和阈值inputWeights=net.IW{1,1};inputbias=net.b{1};%  当前网络层权值和阈值layerWeights=net.LW{2,1};layerbias=net.b{2};%  设置训练参数net.trainParam.show = 50;net.trainParam.lr = 0.05;net.trainParam.mc = 0.9;net.trainParam.epochs = 1000;net.trainParam.goal = 1e-3;
end
%  调用相应算法训练 BP 网络
[net,tr]=train(net,P,T);
end

调用方法

% L-M 优化算法（trainlm）和贝叶斯正则化算法（trainbr），用以训练 BP 网络
% 使其能够拟合某一附加有白噪声的正弦样本数据，渡边笔记
% 神经网络的初始和常数设置，请在shenjin.m修改
clc,clear;close all;% 0.全选，如果数据不多时可以操作
% 1.L-M 优化算法 TRAINLM
% 2.贝叶斯正则化算法 TRAINBR
% 3.动量梯度下降算法 TRAINGDM
% trainbfg---BFGS算法（拟牛顿反向传播算法）训练函数；
% trainbr---贝叶斯归一化法训练函数；
% traincgb---Powell-Beale共轭梯度反向传播算法训练函数；
% traincgp---Polak-Ribiere变梯度反向传播算法训练函数；
% traingd---梯度下降反向传播算法训练函数；
% traingda---自适应调整学习率的梯度下降反向传播算法训练函数；
% traingdm---附加动量因子的梯度下降反向传播算法训练函数；
% traingdx---自适应调整学习率并附加动量因子的梯度下降反向传播算法训练函数；
% trainlm---Levenberg-Marquardt反向传播算法训练函数；
% trainoss---OSS（one step secant）反向传播算法训练函数；
% trainrp---RPROP（弹性BP算法）反向传播算法训练函数；
% trainscg---SCG（scaled conjugate gradient）反向传播算法训练函数；
% trainb---以权值/阈值的学习规则采用批处理的方式进行训练的函数；
% trainc---以学习函数依次对输入样本进行训练的函数；
% trainr---以学习函数随机对输入样本进行训练的函数。
% 当调用train函数时，上述训练函数被用于训练网络
L = 3; %共有三种方法,我觉得够用了
choice = 0;% P 为输入矢量
P = [1 2 3 4 5];
% T 为目标矢量
T = [1 2 3 4 5];%  绘制样本数据点
% 3种全选，比较谁更合适
% choice = 0
if choice == 0for i = 1:L% 神经网络[net,tr] = shenjin(i,P,T);%  对 BP 网络进行仿真A = sim(net,P);title(['神经网络 方法',num2str(i),' 计算结果']);%  计算仿真误差E = T - A;MSE(i) = mse(E); % mse 函数用于计算均方误差% 均方误差（mean-square error, MSE）是反映估计量与被估计量之间差异程度的一种度量endfor i = 1:Ldisp(['神经网络 方法',num2str(i),' 的_均方误差_为 : ',num2str(MSE(i))])end[~,Min] = min(MSE);disp(['神经网络 方法',num2str(Min),' 比较优'])
elsesubplot(1,2,1)plot(P,T,'+');title('原始数据点');% 神经网络[net,tr] = shenjin(i,P,T);% 对 BP 网络进行仿真A = sim(net,P);title(['神经网络 方法',num2str(choice),' 计算结果']);%  计算仿真误差E = T - A;MSE = mse(E); % mse 函数用于计算均方误差% 均方误差（mean-square error, MSE）是反映估计量与被估计量之间差异程度的一种度量disp(['神经网络 方法',num2str(choice),' 的_均方误差_为 : ',num2str(MSE)])
end% 这个用于分析是否可以很好的进行线性回归
% [m,b,r] = postreg(A, P);
% 神经网络的简单使用方法
% P = [1 2 3];
% [Y] = sim(net,P);
% 输入即可输出，维度应与元数据保存一致

18. 贪心算法

% 渡边笔记
% 贪心算法是一种思想，即寻找当前最优解，至使结果最优解
% 贪心算法必须经过证明才可以使用
% 其证明围绕着：整个问题的最优解一定由在贪心策略中存在的子问题的最优解得来的
% 一旦贪心算法得以证明，它就是最快最优的算法
% 实际上，贪心算法适用的情况很少
% 判断一个问题是否适合用贪心法求解，目前还没有一个通用的方法，在竞赛中，需要凭个人的经验来判断
% 贪心算法毕竟是贪心算法，只能缓一时，而不是长久之计，
% 问题的模型、参数对贪心算法求解结果具有决定性作用，这在某种程度上是不能接受的
% 于是聪明的人类就发明了各种智能算法（也叫启发式算法）
% 但在我看来所谓的智能算法本质上就是贪心算法和随机化算法结合
% 例如传统遗传算法用的选择策略就是典型的贪心选择，正是这些贪心算法和随机算法的结合，我们才看到今天各种各样的智能算法
%% 例1 设有n个正整数，将它们连接成一排，组成一个最大的多位整数。
% n=3时，3个整数13，312，343，连成的最大整数为34331213。
% n=4时，4个整数7，13，4，246，连成的最大整数为7424613。
% 此题很容易想到使用贪心法，比如把整数按从大到小的顺序连接起来，例子符合，但测试结果不全对
% 反例：12，121应该组成12121而非12112
% 正确的标准是：先把整数转换成字符串，然后比较a+b和b+a，如果a+b>=b+a，就把a排在b的前面，反之则把a排在b的后面
%% 例2 最短路径问题
clear,clc;close all;
n=4; % 设置点数
x = zeros(1,n); % 产生一个行向量存储坐标
y = zeros(1,n);
luxian = 1:1:n; % 生成1到n的矩阵，存储路线
paixu = 1:1:n;% 这一步是作出随机点的图
for i = 1 : nx(i) = randn * n ;                  %生成正态分布的随机坐标点y(i) = randn * n ;hold ontext(x(i)+0.3,y(i)+0.3,num2str(i))     %用text做好标记，
end% 这一步是计算点与点之间的距离
d = zeros(n) ;  % 生成一个n*n的矩阵用来存储点与点之间的距离，可删去
for i = 1 : nfor j = 1 : nd(i,j) = sqrt( ( x(i) - x(j) ) ^ 2 + ( y(i) - y(j) ) ^ 2);end
endluxian(1) = 1;                             %设置起点，路线的起点是1
num = 1;
for a = 1:(n-2)                         %去除起点和终点需要n-2次判断
paixu(:,1)=[];                       %删除上一个最优点
dis = zeros(1,(n-a));                %用来存剩下各个点的距离for b = 1:(n-a)                     %用来获取剩下各个点的距离dis(b) = d (num, paixu(b));end
num1 = find( dis == min(dis) );     %根据距离最小得到最优点位置
t = paixu(1);                      %通过中间变量互换位置paixu(1) = paixu(num1);       %最优点位置换到第一个paixu(num1) = t;num = paixu(1);                    %获取下次进行操作的数luxian(a+1) = paixu(1);              %将最优点存入luxian数组（空出开头）
end
luxian(n) = paixu(2);  %补上最后一个点plot(x(luxian),y(luxian),'--o') ; %标出点用虚线相连得到路径
grid on;title('贪心算法计算最优路径');

19. 小波分析

%%初始化程序
clear,clc
t1=clock;%% 载入噪声信号数据，数据为.mat格式，并且和程序放置在同一个文件夹下
x = 0:0.06:10;
YSJ  = sin(x)+0.2*rand(size(x));plot(YSJ);
title('原始数据');%% 数据预处理，数据可能是存储在矩阵或者是EXCEL中的二维数据，衔接为一维的，如果数据是一维数据，此步骤也不会影响数据
[c,l]=size(YSJ);
Y=[];
for i=1:cY=[Y,YSJ(i,:)];
end
[c1,l1]=size(Y);
X=[1:l1];%% 绘制噪声信号图像
figure(1);
plot(X,Y);
xlabel('横坐标');
ylabel('纵坐标');
title('原始信号');%% 硬阈值处理
lev=3;
xd=wden(Y,'heursure','h','one',lev,'db4');%硬阈值去噪处理后的信号序列
figure(2)
plot(X,xd)
xlabel('横坐标');
ylabel('纵坐标');
title('硬阈值去噪处理')
set(gcf,'Color',[1 1 1])%% 软阈值处理
lev=3;
xs=wden(Y,'heursure','s','one',lev,'db4');%软阈值去噪处理后的信号序列
figure(3)
plot(X,xs)
xlabel('横坐标');
ylabel('纵坐标');
title('软阈值去噪处理')
set(gcf,'Color',[1 1 1])
%% 固定阈值后的去噪处理
lev=3;
xz=wden(Y,'sqtwolog','s','sln',lev,'db4');%固定阈值去噪处理后的信号序列
figure(4)
plot(X,xz);
xlabel('横坐标');
ylabel('纵坐标');
title('固定阈值后的去噪处理')
set(gcf,'Color',[1 1 1])
%% 计算信噪比SNR
Psig=sum(Y*Y')/l1;
Pnoi1=sum((Y-xd)*(Y-xd)')/l1;
Pnoi2=sum((Y-xs)*(Y-xs)')/l1;
Pnoi3=sum((Y-xz)*(Y-xz)')/l1;
SNR1=10*log10(Psig/Pnoi1);
SNR2=10*log10(Psig/Pnoi2);
SNR3=10*log10(Psig/Pnoi3);
%% 计算均方根误差RMSE
RMSE1=sqrt(Pnoi1);
RMSE2=sqrt(Pnoi2);
RMSE3=sqrt(Pnoi3);
%% 输出结果
disp('-------------三种阈值设定方式的降噪处理结果---------------');
disp(['硬阈值去噪处理的SNR=',num2str(SNR1),'，RMSE=',num2str(RMSE1)]);
disp(['软阈值去噪处理的SNR=',num2str(SNR2),'，RMSE=',num2str(RMSE2)]);
disp(['固定阈值后的去噪处理SNR=',num2str(SNR3),'，RMSE=',num2str(RMSE3)]);
t2=clock;
tim=etime(t2,t1);
disp(['------------------运行耗时',num2str(tim),'秒-------------------'])

20. 遗传算法

我是非常不建议使用这个的，太复杂太麻烦了，每次换算遗传因子都要想半天，而且代码冗余，又长又臭

20.1. TSP 示例

cross.m

%交叉操作函数  cross.m
function [A,B]=cross(A,B)
L=length(A);
if L<10W=L;
elseif ((L/10)-floor(L/10))>=rand&&L>10W=ceil(L/10)+8;
elseW=floor(L/10)+8;
end
%%W为需要交叉的位数
p=(L-W+1);%随机产生一个交叉位置
%fprintf('p=%d ',p);%交叉位置
for i=1:Wx=find(A==B(1,p+i-1));y=find(B==A(1,p+i-1));[A(1,p+i-1),B(1,p+i-1)]=exchange(A(1,p+i-1),B(1,p+i-1));[A(1,x),B(1,y)]=exchange(A(1,x),B(1,y));
endend

exchange.m

%对调函数 exchange.mfunction [x,y]=exchange(x,y)
temp=x;
x=y;
y=temp;end

fit.m

%适应度函数fit.m，每次迭代都要计算每个染色体在本种群内部的优先级别，类似归一化参数。越大约好！
function fitness=fit(len,m,maxlen,minlen)
fitness=len;
for i=1:length(len)fitness(i,1)=(1-(len(i,1)-minlen)/(maxlen-minlen+0.0001)).^m;
end

Mutation.m

%变异函数 Mutation.mfunction a=Mutation(A)
index1=0;index2=0;
nnper=randperm(size(A,2));
index1=nnper(1);
index2=nnper(2);
%fprintf('index1=%d ',index1);
%fprintf('index2=%d ',index2);
temp=0;
temp=A(index1);
A(index1)=A(index2);
A(index2)=temp;
a=A;end

myLength.m

%染色体的路程代价函数  mylength.m
function len=myLength(D,p)%p是一个排列
[N,NN]=size(D);
len=D(p(1,N),p(1,1));
for i=1:(N-1)len=len+D(p(1,i),p(1,i+1));
endend

plot_route.m

%连点画图函数 plot_route.mfunction plot_route(a,R)
scatter(a(:,1),a(:,2),'rx');
hold on;
plot([a(R(1),1),a(R(length(R)),1)],[a(R(1),2),a(R(length(R)),2)]);
hold on;
for i=2:length(R)x0=a(R(i-1),1);y0=a(R(i-1),2);x1=a(R(i),1);y1=a(R(i),2);xx=[x0,x1];yy=[y0,y1];plot(xx,yy);hold on;
endend

调用方法

% 渡边笔记
% 多种群遗传算法解决 TSP 问题
% TSP问题即旅行商问题
% 经典的TSP可以描述为
% 一个商品推销员要去若干个城市推销商品，该推销员从一个城市出发，需要经过所有城市后，回到出发地
% 应如何选择行进路线，以使总的行程最短。从图论的角度来看，该问题实质是在一个带权完全无向图中，找一个权值最小的哈密尔顿回路clear,clc;close all;
% 输入参数
M=100;               %种群的个数
ITER=200;            %迭代次数
m=2;                 %适应值归一化淘汰加速指数
Pc=0.8;              %交叉概率
Pmutation=0.05;      %变异概率
% 城市坐标
pos=[0 01 1-1 -12 3] ;%生成城市之间距离矩阵
[N,~] = size(pos); %%城市的个数
D=zeros(N,N);
for i=1:Nfor j=i+1:Ndis=(pos(i,1)-pos(j,1)).^2+(pos(i,2)-pos(j,2)).^2;D(i,j)=dis^(0.5);D(j,i)=D(i,j);end
end%生成初始群体
popm=zeros(M,N);
for i=1:Mpopm(i,:)=randperm(N);%随机排列，比如[2 4 5 6 1 3]
end
%随机选择一个种群
R=popm(1,:);%初始化种群及其适应函数
fitness = zeros(M,1);
len = zeros(M,1);
for i=1:M   % 计算每个染色体对应的总长度len(i,1)=myLength(D,popm(i,:));
end
maxlen=max(len);%最大回路
minlen=min(len);%最小回路
fitness=fit(len,m,maxlen,minlen);
rr=find(len==minlen);%找到最小值的下标，赋值为rr
R=popm(rr(1,1),:);%提取该染色体，赋值为Rfitness=fitness/sum(fitness);
distance_min=zeros(ITER+1,1);  %%各次迭代的最小的种群的路径总长
nn=M;
iter=0;
while iter<=ITER%选择操作p=fitness./sum(fitness);q=cumsum(p);%累加for i=1:(M-1)len_1(i,1)=myLength(D,popm(i,:));r=rand;tmp=find(r<=q);popm_sel(i,:)=popm(tmp(1),:);end[fmax,indmax]=max(fitness);%求当代最佳个体popm_sel(M,:)=popm(indmax,:);%交叉操作nnper=randperm(M);
%    A=popm_sel(nnper(1),:);
%    B=popm_sel(nnper(2),:);for i=1:M*Pc*0.5A=popm_sel(nnper(i),:);B=popm_sel(nnper(i+1),:);[A,B]=cross(A,B);
%      popm_sel(nnper(1),:)=A;
%      popm_sel(nnper(2),:)=B;popm_sel(nnper(i),:)=A;popm_sel(nnper(i+1),:)=B;end
%变异操作for i=1:Mpick=rand;while pick==0pick=rand;endif pick<=Pmutationpopm_sel(i,:)=Mutation(popm_sel(i,:));endend%%求适应度函数NN=size(popm_sel,1);len=zeros(NN,1);for i=1:NNlen(i,1)=myLength(D,popm_sel(i,:));endmaxlen=max(len);minlen=min(len);distance_min(iter+1,1)=minlen;fitness=fit(len,m,maxlen,minlen);rr=find(len==minlen);R=popm_sel(rr(1,1),:);
%    for i=1:N
%        fprintf('%d ',R(i));
%   end
%    fprintf('\n');popm=[];popm=popm_sel;iter=iter+1;%pause(1);
end%画出所有城市坐标
subplot(1,3,1);
scatter(pos(:,1),pos(:,2),'rx');fprintf('最短路程是 ');
for i=1:Nfprintf('%d ',R(i));%把R顺序打印出来
end
fprintf('\n最短距离是 %d\n',minlen);
title('点图');
subplot(1,3,2); plot_route(pos,R);  % 路程图
title('路程图');
subplot(1,3,3); plot(distance_min); % 最短距离随迭代次数的变化
title('最短距离的变化');

21. 蚁群算法

21.1. TSP 示例

yiyu%% 蚁群算法求解TSP问题%% 清空环境变量
clear
clc%% 导入数据
load citys_data.mat%% 计算城市间相互距离
n = size(citys,1);
D = zeros(n,n);
for i = 1:nfor j = 1:nif i ~= jD(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2));elseD(i,j) = 1e-4;endend
end%% 初始化参数
m = 50;                              % 蚂蚁数量
alpha = 1;                           % 信息素重要程度因子
beta = 5;                            % 启发函数重要程度因子
rho = 0.1;                           % 信息素挥发因子
Q = 1;                               % 常系数
Eta = 1./D;                          % 启发函数
Tau = ones(n,n);                     % 信息素矩阵
Table = zeros(m,n);                  % 路径记录表
iter = 1;                            % 迭代次数初值
iter_max = 200;                      % 最大迭代次数
Route_best = zeros(iter_max,n);      % 各代最佳路径
Length_best = zeros(iter_max,1);     % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);      % 各代路径的平均长度%% 迭代寻找最佳路径
while iter <= iter_max% 随机产生各个蚂蚁的起点城市start = zeros(m,1);for i = 1:mtemp = randperm(n);start(i) = temp(1);endTable(:,1) = start;% 构建解空间citys_index = 1:n;% 逐个蚂蚁路径选择for i = 1:m% 逐个城市路径选择for j = 2:ntabu = Table(i,1:(j - 1));           % 已访问的城市集合(禁忌表)allow_index = ~ismember(citys_index,tabu);allow = citys_index(allow_index);  % 待访问的城市集合P = allow;% 计算城市间转移概率for k = 1:length(allow)P(k) = Tau(tabu(end),allow(k))^alpha * Eta(tabu(end),allow(k))^beta;endP = P/sum(P);% 轮盘赌法选择下一个访问城市Pc = cumsum(P);target_index = find(Pc >= rand);target = allow(target_index(1));Table(i,j) = target;endend% 计算各个蚂蚁的路径距离Length = zeros(m,1);for i = 1:mRoute = Table(i,:);for j = 1:(n - 1)Length(i) = Length(i) + D(Route(j),Route(j + 1));endLength(i) = Length(i) + D(Route(n),Route(1));end% 计算最短路径距离及平均距离if iter == 1[min_Length,min_index] = min(Length);Length_best(iter) = min_Length;Length_ave(iter) = mean(Length);Route_best(iter,:) = Table(min_index,:);else[min_Length,min_index] = min(Length);Length_best(iter) = min(Length_best(iter - 1),min_Length);Length_ave(iter) = mean(Length);if Length_best(iter) == min_LengthRoute_best(iter,:) = Table(min_index,:);elseRoute_best(iter,:) = Route_best((iter-1),:);endend% 更新信息素Delta_Tau = zeros(n,n);% 逐个蚂蚁计算for i = 1:m% 逐个城市计算for j = 1:(n - 1)Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i);endDelta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i);endTau = (1-rho) * Tau + Delta_Tau;% 迭代次数加1，清空路径记录表iter = iter + 1;Table = zeros(m,n);
end%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);%% 绘图
figure(1)
plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],...[citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-');
grid on
for i = 1:size(citys,1)text(citys(i,1),citys(i,2),['   ' num2str(i)]);
end
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),'       起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),'       终点');
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['蚁群算法优化路径(最短距离:' num2str(Shortest_Length) ')'])
figure(2)
plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r:')
legend('最短距离','平均距离')
xlabel('迭代次数')
ylabel('距离')
title('各代最短距离与平均距离对比')

citys_data.mat

1304 2312
3639    1315
4177    2244
3712    1399
3488    1535
3326    1556
3238    1229
4196    1004
4312    790
4386    570
3007    1970
2562    1756
2788    1491
2381    1676
1332    695
3715    1678
3918    2179
4061    2370
3780    2212
3676    2578
4029    2838
4263    2931
3429    1908
3507    2367
3394    2643
3439    3201
2935    3240
3140    3550
2545    2357
2778    2826
2370    2975

22. 元胞自动机

22.1. 森林火灾示例

clc,clear;
n=200;
Pltg = 5e-6;
Pgrw = 1e-2;
NW = [n 1:n-1];
SE = [2:n 1];
veg = zeros(n);
imh = image( cat(3,(veg == 1),(veg == 2),zeros(n)) );for i = 1:30num = (veg(NW,:)==1) + (veg(:,NW)==1) + (veg(:,SE)==1) + (veg(SE,:)==1);veg = 2*( (veg==2) | (veg==0 & rand(n)<Pgrw) ) - ...((veg==2) & (num>0|rand(n)<Pltg ) );set(imh,'cdata',cat(3,(veg==1),(veg==2),zeros(n)));drawnow
end

22.2. 交通流示例

gaplength.m

function gaps = gaplength(x,L)
ncar = length(x);
gaps = inf*ones(1,ncar);
if ncar>1gaps = x([2:end 1]) -x;gaps(gaps<0) = gaps(gaps<0) +L;
end

ns.m

function [flux,vmean] = ns(rho,p,L,tmax)ncar = ceil(L*rho);
vmax = 5;
x = sort(randperm(L,ncar));
v = vmax*ones(1,ncar);
flux = 0;vmean = 0;h = plotcirc(L,x,0.5);for t = 1:tmaxv = min(v+1,vmax);gaps = gaplength(x,L);v = min(v,gaps - 1);vdrops = (rand(1,ncar)<p);v = max(v - vdrops,0);x = x + v;passed = x>L;x(passed) = x(passed)-L;if t>tmax/2flux = flux +sum(v)/L;vmean = vmean+mean(v);endplotcirc(L,x,0.5,h);
endflux = flux/(tmax/2);
vmean = vmean/(tmax/2);

plotcirc.m

function h=plotcirc(L,x,dt,h)
W = 0.05;R=1;
ncar = length(x);theta = 0:2*pi/L:2*pi;
xc = cos(theta);yc = sin(theta);
xin = (R-W/2)*xc;yin = (R-W/2)*yc;
xot = (R+W/2)*xc;yot = (R+W/2)*yc;xi = [xin(x);xin(x+1);xot(x+1);xot(x)];
yi = [yin(x);yin(x+1);yot(x+1);yot(x)];if nargin == 3color = randperm(ncar);h = fill(xi,yi,color);hold onplot([xin;xot],[yin;yot],'k','linewidth',1.5);plot([xin;xot]',[yin;yot]','k','linewidth',1.5);
elsefor i=1:ncarset(h(i),'xdata',xi(:,i),'ydata',yi(:,i));end
end
pause(dt)

调用方法

clc;
clear;
[flux,vmean] = ns(0.2,1e-9,25,30);

23. 主成分分析

% 渡边笔记 主成分分析
% 当研究的问题涉及多个变量，并且变量间相关性明显，即包含的信息有所重叠时
% 可以考虑用主成分分析的方法，这样更容易抓住事物的主要矛盾，使问题简化
clc,clear;close all;
% 与主成分分析相关的 MATLAB 函数有pcacov、princomp函数% pcacov函数用来根据 协方差矩阵  或 相关系数矩阵 进行主成分分析
% [COEFF,latent,explained]=pcacov（V）
% V 是总体或样本的协方差矩阵或相关系数矩阵，对于p维总体，V是pxp的矩阵
% 输出参数 COEFF 是p个主成分的系数矩阵，它是pxp的矩阵，它的第i列是第i个主成分的系数向量
% 输出参数 latent 是p个主成分的方差构成的向量，即V的p个特征值的大小(从大到小)构成的向量
% 输出参数 explained 是p个主成分的贡献率向量，已经转化为百分比
V = [1     0.79    0.36    0.76    0.25    0.510.79  1       0.31    0.55    0.17    0.350.36  0.31    1       0.35    0.64    0.580.76  0.55    0.35    1       0.16    0.380.25  0.17    0.64    0.16    1       0.630.51  0.35    0.58    0.38    0.63    1 ];
[a,~] = size(V);
[COEFF_1,latent,explained] = pcacov(V);
result_1(1,:) = {'特征值', '差值', '贡献率', '累积贡献率'};
result_1(2: a+1 ,1) = num2cell(latent);
% diff函数式用于求导数和差分的.
result_1(2: a ,2) = num2cell(-diff(latent));
% cumsum函数通常用于计算一个数组各行的累加值。
result_1(2: a+1 ,3:4) = num2cell([explained, cumsum(explained)]);
% disp(result_1);% princomp函数用来根据样本观测值矩阵进行主成分分析，其调用格式如下：
% [COEFF,SCORE,latent,tsquare] = princomp(X,'econ')
% 根据样本观测值矩阵X进行主成分分析。输入参数X是n行p列的矩阵，每一行对应一个观测(样品)，每一列对应一个变量
% 输出参数COEFF是p个主成分析的系数矩阵，是pxp的矩阵，它的第i列对应第i个主成分的系数向量
% 输出参数SCORE是n个样品的p个主成分得分矩阵，它是n行p列的矩阵
% 每一行对应一个观测，每一列对应一个主成分，第i行第j列元素表示第i个样品的第j个主成分得分
% SCORE与X是一一对应的关系，是X在新坐标系中的数据，可以通过X*系数矩阵得到
% 返回样本协方差矩阵的特征值向量latent，它是由p个特征值构成的列向量，其中特征值按降序排列。
% 返回一个包含n个元素的列向量tsquare，它的第i个元素的第i个观测对应的霍特林T^2统计量
% 表述了第i个观测与数据集(样本观测矩阵)的中心之间的距离，可用来寻找远离中心的极端数据
% 通过设置参数‘econ’参数，使得当 n <= p 时
% 只返回latent中的前n-1个元素(去掉不必要的0元素)及COEFF和SCORE矩阵中相应的列
X = [2000 1000 3000 900 951900  990  3000 700 662000   900  3400 800 56];
% 行为样本，列为特征量
XZ = zscore(X); % 数据标准化
[COEFF_2,SCORE,latent,tsquare] = pca(XZ);
[m, n] = size(X); % 求X的行数和列数
result_2 = cell(m, 4); % 定义一个n+1行，4列的元胞数组
result_2(1,:) = {'特征值', '差值', '贡献率', '累积贡献率'};
explained = 100*latent/sum(latent);% 计算贡献率
%result_2中第1列的第2行到最后一行存放的数据(latent)特征值
result_2(2:end,1) = num2cell(latent);
%result1中第2列的第2行到倒数第2行存放的数据(latent的方差，特征值的方差)
result_2(2:end-1,2) = num2cell(-diff(latent));
%result1中第3列和第4列的第2行到最后一行分别存放主成分的贡献率和累积贡献率
result_2(2:end,3:4) = num2cell([explained, cumsum(explained)]);
disp(result_2);
[b_1,b_2] = size(X);
[b_3,b_4] = size(SCORE);
disp(SCORE);

24. 主成分回归

这是我论文的代码，找不到原生的主成分回归了

clc,clear,close all;
% 渡辺论文代码 m_1
% 文件的路径
filename_1 = 'first_shuju.xlsx';
sheet_1 = '';
% Excle 的 工作表 名称
range_1 = 'B43:Q58';
range_2 = 'B1:T1';
range_3 = 'G3:G18';
Year_0 = 2004;
years(1) = Year_0;
for j = 2:13years(j) = years(j-1)+1;
end[num_1, ~,~] = xlsread(filename_1,sheet_1,range_1);
[~, num_2,~] = xlsread(filename_1,sheet_1,range_2);
[num_3, ~,~] = xlsread(filename_1,sheet_1,range_3);
data = num_1(1:13,:);X = num_1(1:13,:);
XZ = zscore(X); % 数据标准化
[COEFF_2,SCORE,latent,tsquare] = pca(XZ);
[m, n] = size(X); % 求X的行数和列数
result_2 = cell(m, 4); % 定义一个n+1行，4列的元胞数组
result_2(1,:) = {'特征值', '差值', '贡献率', '累积贡献率'};
explained = 100*latent/sum(latent);% 计算贡献率
%result_2中第1列的第2行到最后一行存放的数据(latent)特征值
result_2(2:end,1) = num2cell(latent);
%result1中第2列的第2行到倒数第2行存放的数据(latent的方差，特征值的方差)
result_2(2:end-1,2) = num2cell(-diff(latent));
%result1中第3列和第4列的第2行到最后一行分别存放主成分的贡献率和累积贡献率
result_2(2:end,3:4) = num2cell([explained, cumsum(explained)]);
disp(result_2);
disp(COEFF_2);
xlswrite('result.xlsx',result_2);
xlswrite('SCORE.xlsx',SCORE);y = num_3(1:13);
% x,y 需要转置才可以使用
num_1_2 = num_1(1:13,:);
ppp = 0;
for i=1:13for ii = 1:6pp(i,ii) = num_1_2(i,:) * COEFF_2(:,ii);end
end
SCORE = pp;
X = [ones(size(y)) SCORE(:,1) SCORE(:,2) SCORE(:,3) SCORE(:,4) ...SCORE(:,5) SCORE(:,6)]; % 定义回归方程[B,BINT,R,RINT,STATS] = regress(y,X); % 回归
% [B,BINT,R,RINT,STATS] = REGRESS(Y,X,alpha)
% b 为回归系数向量
% BINT 回归系数的估计区间
% R 残差
% RINT 置信区间
% STATS 用于检验回归模型的统计量。有4个数值：判定系数r2，F统计量观测值，检验的p的值，误差方差的估计
% r2越接近1，回归方程越显著；F>F1−α(k,n−k−1)时拒绝H0，F越大，回归方程越显著；p<α时拒绝H0
% ALPHA 显著性水平（缺少时默认0.05）
y2 = B(2).* SCORE(:,1)' + B(3).*SCORE(:,2)' + B(4).*SCORE(:,3)' + B(5).*SCORE(:,4)' ...+ B(6).*SCORE(:,5)' + B(7).*SCORE(:,6)' + B(1)';   %预测值
plot([2004:2019],num_3','ro-',years,y2,'b*');
legend('生态用水量', '主成分回归');
xlabel('时间/年');
ylabel('城市生态用水量/亿立方米');
disp('回归分析统计量为');
disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 这么多原始数据 解释，最优为 1
disp([ 'F : ',num2str(STATS(2))]);
disp([ '检验 p : ',num2str(STATS(3))]); % 日常小于 0.05 最好
disp([ '误差方差 : ',num2str(STATS(4))]); % 越小越好
xlswrite('zhuchengf_huigui.xlsx',B);
hold on;num_1_2 = num_1(14:16,:);
ppp = 0;
for n = 1:3for ii = 1:6pp_1(n,ii) = num_1_2(n,:) * COEFF_2(:,ii);end
end
SCORE = pp_1;
y3 = B(2).* SCORE(:,1)' + B(3).*SCORE(:,2)' + B(4).*SCORE(:,3)' + B(5).*SCORE(:,4)' ...+ B(6).*SCORE(:,5)' + B(7).*SCORE(:,6)' + B(1)';   %预测值
plot([2017:2019],y3,'m*');
line([2016.5,2016.5],[0,17],'linestyle','--');
axis([2004 2019 0 17]);
legend('城市生态用水量-原始数据', '主成分回归在数据集的结果','主成分回归在测试集的结果');
xlswrite('save_4.xlsx',y3);

25. 分治算法

25.1. 快速排序

QuickSort.m

function a=QuickSort(a,leftIndex,rightIndex)
% a是待排序序列
%leftIndex和rightIndex是开始的左右两个边界
%%----------------------------------------------------------%%
% if leftIndex>rightIndex
%     break;
% end
if leftIndex<rightIndexi=leftIndex;j=rightIndex;temp=a(i);%选定的这个数挖掉，相当于挖坑while i<jwhile (i<j)&&(a(j)>=temp)%从右往左，找到第一个小于设定的数，j=j-1;enda(i)=a(j);%将找到的第一个小于设定的数填坑到最开始挖的坑里面去while (i<j)&&(a(i)<=temp)%从做到由，找到第一个大于选定的数i=i+1;enda(j)=a(i);%将找到的第一个大于选定的数填入上一步挖的坑里面去
%     if i==j
%         a(j)=temp;
%     endenda(j)=temp;%最后，i=j，将选定的数再填到上一步挖的坑里面去a=QuickSort(a,leftIndex,j-1);%对左边序列进行递归a=QuickSort(a,i+1,rightIndex);%对右边序列进行递归
end
end

调用方法

% 分治算法，渡边笔记
% 分治算法是一种思路，把大问题归纳为小问题来解决
% 使用分治法，问题应该满足
% 该问题的规模缩小到一定的程度就可以容易地解决
% 该问题可以分解为若干个规模较小的相同问题，即该问题具有最优子结构性质
% 利用该问题分解出的子问题的解可以合并为该问题的解
% 该问题所分解出的各个子问题是相互独立的，即子问题之间不包含公共的子子问题
% 第一条特征是绝大多数问题都可以满足的，因为问题的计算复杂性一般是随着问题规模的增加而增加
% 第二条特征是应用分治法的前提它也是大多数问题可以满足的，此特征反映了递归思想的应用
% 第三条特征是关键，能否利用分治法完全取决于问题是否具有第三条特征，如果具备了第一条和第二条特征，而不具备第三条特征，则可以考虑用贪心法或动态规划法。
% 第四条特征涉及到分治法的效率，如果各子问题是不独立的则分治法要做许多不必要的工作，
% 重复地解公共的子问题，此时虽然可用分治法，但一般用动态规划法较好
clc,clear;close all;
a=[72 57 88 60 42 83 73 48 85]; % 输入排序数leftIndex=1;
[~, rightIndex]=size(a);
disp(['未排序的序列为：',num2str(a)])
a=QuickSort(a,leftIndex,rightIndex);
disp(['未排序的序列为：',num2str(a)])

25.2. 最近点对

cloest.m

%分治求最近点对
function [d,x1,x2] = cloest(S,low,high)
P=[];
if(high - low == 1)[d,x1,x2] = deal(pdist2(S(low,:),S(high,:)),S(low,:),S(high,:));
elseif(high - low == 2)d1 = pdist2(S(low,:),S(low + 1,:));d2 = pdist2(S(low + 1,:),S(high,:));d3 = pdist2(S(low,:),S(high,:));if((d1 < d2) && (d1 < d3))[d,x1,x2] = deal(d1,S(low,:),S(low + 1,:));elseif(d2 < d3)[d,x1,x2] = deal(d2,S(low+1,:),S(high,:));else[d,x1,x2] = deal(d3,S(low,:),S(high,:));end
elsemid = floor((low+high)/2);[d1,x11,x12] = cloest(S,low,mid);[d2,x21,x22] = cloest(S,mid+1,high);if(d1 <= d2)[d,x1,x2] = deal(d1,x11,x12);else[d,x1,x2] = deal(d2,x21,x22);endindex = 0;for i = mid:-1:lowif(S(mid,1) - S(i,1) >= d)break;endindex = index + 1;P(index,:) = S(i,:);endfor i = mid+1:highif(S(i,1) - S(mid,1) >= d)break;endindex = index + 1;P(index,:) = S(i,:);endsortrows(P,2);for i = 1:indexfor j = i+1:indexif(P(j,2) - P(i,2) >= d)break;elsed3 = pdist2(P(i,:),P(j,:));if(d3 < d)[d,x1,x2] = deal(d3,P(i,:),P(j,:));endendendend
end

调用方法

% 分治法求最近点，渡边笔记
clear;clc
n = 20;
% 随机生成20个点
A=rand(n,2)*10;
% 将20个点按横坐标升序排列
A = sortrows(A,1);
% 分治法求随机点的最近点对
[mindist1,y1,y2] = cloest(A,1,n);
% 使用红色的点标记随机点
x=A(:,1);
y=A(:,2);
for i = 1:nplot(x,y,'r.'); hold ontext(A(i,1),A(i,2),num2str(i));hold on
end
%使用绿色十字标记分治法的最近点对
plot(y1(1),y1(2),'go', 'linewidth', 2);hold on
plot(y2(1),y2(2),'go', 'linewidth', 2);hold on

26. 杂项

26.1. 有向图

clc,clear,close all;
s = [1 1 1 2 2 3 3 4 5 5 6 7];
t = [2 4 8 3 7 4 6 5 6 8 7 8];
weights = [10 10 1 10 1 10 1 1 12 12 12 12];
% names = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H'};
% G = digraph(s,t,weights,names);
G = digraph(s,t,weights);
plot(G,'Layout','force','EdgeLabel',G.Edges.Weight)

26.2. 数据拟合

懒得修改了，反正方法也就那样

clc,clear
% 产生数据
lamuda_0 = xlsread('shuju.xlsx','','B3:U23');
lamuda_test = xlsread('shuju.xlsx','','A3:A23');
lamuda_test = lamuda_test';
x = xlsread('shuju.xlsx','','V3:V23');
% 拟合
for j = 1:21lamuda = 355 - lamuda_test(j);p = fittype('x * exp( S*lamuda )');f = fit(x,lamuda_0(j,:)',p);
%     plot(f,x,lamuda_0(j,:)');S(j) = f.S;
end

26.3. 标准化

% 标准化处理，渡边笔记
% 归一化的依据非常简单，不同变量往往量纲不同
% 归一化可以消除量纲对最终结果的影响，使不同变量具有可比性
% 如两个人体重差10KG，身高差0.02M
% 在衡量两个人的差别时体重的差距会把身高的差距完全掩盖，归一化之后就不会有这样的问题
% 数据归一化应该针对对于的 y，而不是针对每条数据，针对每条数据是完全没有意义的，因为只是等比例缩放，对之后的分类没有任何作用
clc,clear;close all;
x = [1 2 31 2 3];
% Min-max 极大极小标准化
% 当有新数据加入时，可能导致xmax和xmin的变化，需要重新计算
X_min_max = mapminmax(x, 0, 1); % 默认标准化区间为 0-1
disp('极大极小标准化为 :')
disp(X_min_max);% z-score 标准化
% 其结果没有任何意义，只能用于比较
% X_z_score = zscore(x);
% disp('z-score 标准化 :')
% disp(X_z_score);

26.4. excel 数据转置

clc,clear
% 省会城市 顺序表
% 北京 天津 石家庄 太原 呼和浩特 沈阳 长春 哈尔滨 上海 南京
% 杭州 合肥 福州 南昌 济南 郑州 武汉 长沙 广州 南宁
% 海口 重庆 成都 贵阳 昆明 拉萨 西安 兰州 西宁 银川
% 乌鲁木齐
filename_lujin_1 = 'C:\Users\99791\Desktop\生态用水\数据\国家统计局\分省\';
% 文件的前面的路径名称
filename_name_1 = '分省_城市绿地面积.xls';
% 文件的名字
filename_1 = [filename_lujin_1,filename_name_1];
sheet_1 = '';
% Excle 的 工作表 名称
range_1 = 'B5:Q35';
% 读取的范围
[num_1, tex_1, raw_1] = xlsread(filename_1,sheet_1,range_1);
range_1_1 = 'A2';
[num_2, tex_2, raw_2] = xlsread(filename_1,sheet_1,range_1_1);
num_1 = num_1';
dubian_bianliang_1 = num_1;
% 同上
filename_lujin_2 = 'C:\Users\99791\Desktop\生态用水\数据\';
filename_name_2 = 'dubian_shuju.xlsx';
filename_2 = [filename_lujin_2,filename_name_2];
sheet_2 = '';
t_1 = 21;
t_2 = 'W';
% 改变写入 列
range_2 = [t_2,'3:',t_2,'18'];
for i = 1:31dubian_bianliang_2 = dubian_bianliang_1(:,i);xlswrite(filename_2,dubian_bianliang_2,sheet_2,range_2)range_2 = [t_2,num2str(3+i*t_1),':',t_2,num2str(18+i*t_1)];
end
range_2_1 = [t_2,'1'];
xlswrite(filename_2,tex_2,sheet_2,range_2_1)

26.5. 线性规划

lingo 还是比 matlab 强，下面的方法好久不用了

% function 函数解多元方程 最小值
% function 基本语法
% [x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
% x的返回值是决策向量x的取值，fval的返回值是目标函数f(x)的取值
% fun是用 M 文件定义的函数f(x),代表了(非)线性目标函数
% x0是x的初始值
% A,b,Aeq,beq 定义了线性约束 ,如果没有线性约束，则A=[],b=[],Aeq=[],beq=[]
% lb和ub是变量x的下界和上界，如果下界和上界没有约束，则lb=[],ub=[],也可以写成lb的各分量都为 -inf,ub的各分量都为inf
% nonlcon是用 M 文件定义的非线性向量函数约束
% options 定义了优化参数，不填写表示使用 Matla b默认的参数设置
clc,clear;close all;
options = optimset('Algorithm','active-set'); % 可以不设置
[X,Y,EXITFLAG] = fmincon(@(x)(   x(1)^3 +  x(2)^3 + x(3)^3 +x(4)^3  ),...
[0 0 0 0],[],[],[],[],...
[-2 -2 -2 -2],[-1 -1 -1 -1]... % 设置范围
,[],options);
disp(['函数的满足最小值的 X 解为 : ',num2str(X)]);
disp(['函数的最小值的 Y 解为 : ',num2str(Y)]);