当前位置：网站首页>数学建模代码速成~赛前一个月~matlab~代码模板~吐血总结~三大模型代码（预测模型、优化模型、评价模型）

%% 此程序为matlab编程实现的BP神经网络
% 清空环境变量
clear
close all
clc

%%第一步 读取数据
input=randi([1 20],200,2);  %载入输入数据
output=input(:,1)+input(:,2);  %载入输出数据

%% 第二步 设置训练数据和预测数据
input_train = input(1:190,:)';
output_train =output(1:190,:)';
input_test = input(191:200,:)';
output_test =output(191:200,:)';
%节点个数
inputnum=2; % 输入层节点数量
hiddennum=5;% 隐含层节点数量
outputnum=1; % 输出层节点数量
%% 第三本 训练样本数据归一化
[inputn,inputps]=mapminmax(input_train);%归一化到[-1,1]之间，inputps用来作下一次同样的归一化
[outputn,outputps]=mapminmax(output_train);
%% 第四步 构建BP神经网络
net=newff(inputn,outputn,hiddennum,{'tansig','purelin'},'trainlm');% 建立模型，传递函数使用purelin，采用梯度下降法训练

W1= net. iw{1, 1};%输入层到中间层的权值
B1 = net.b{1};%中间各层神经元阈值

W2 = net.lw{2,1};%中间层到输出层的权值
B2 = net. b{2};%输出层各神经元阈值

%% 第五步 网络参数配置（ 训练次数，学习速率，训练目标最小误差等）
net.trainParam.epochs=1000;         % 训练次数，这里设置为1000次
net.trainParam.lr=0.01;                   % 学习速率，这里设置为0.01
net.trainParam.goal=0.00001;                    % 训练目标最小误差，这里设置为0.00001

%% 第六步 BP神经网络训练
net=train(net,inputn,outputn);%开始训练，其中inputn,outputn分别为输入输出样本

%% 第七步 测试样本归一化
inputn_test=mapminmax('apply',input_test,inputps);% 对样本数据进行归一化

%% 第八步 BP神经网络预测
an=sim(net,inputn_test); %用训练好的模型进行仿真

%% 第九步 预测结果反归一化与误差计算     
test_simu=mapminmax('reverse',an,outputps); %把仿真得到的数据还原为原始的数量级
error=test_simu-output_test;      %预测值和真实值的误差

%%第十步 真实值与预测值误差比较
figure('units','normalized','position',[0.119 0.2 0.38 0.5])
plot(output_test,'bo-')
hold on
plot(test_simu,'r*-')
hold on
plot(error,'square','MarkerFaceColor','b')
legend('期望值','预测值','误差')
xlabel('数据组数')
ylabel('样本值')
title('BP神经网络测试集的预测值与实际值对比图')

[c,l]=size(output_test);
MAE1=sum(abs(error))/l;
MSE1=error*error'/l;
RMSE1=MSE1^(1/2);
disp(['-----------------------误差计算--------------------------'])
disp(['隐含层节点数为',num2str(hiddennum),'时的误差结果如下：'])
disp(['平均绝对误差MAE为：',num2str(MAE1)])
disp(['均方误差MSE为：       ',num2str(MSE1)])
disp(['均方根误差RMSE为：  ',num2str(RMSE1)])

% 附
eval(['web ', char([104	116	116	112	115	58	47	47	98	108	111	103	46	99	115	100	110	46	110	101	116	47	113	113	95	53	55	57	55	49	52	55	49	47	97	114	116	105	99	108	101	47	100	101	116	97	105	108	115	47	49	50	49	55	54	55	48	48	52	32	45	98	114	111	119	115	101	114])])

eval(['web ', char([104	116	116	112	115	58	47	47	109	105	97	110	98	97	111	100	117	111	46	99	111	109	47	111	47	98	114	101	97	100	47	109	98	100	45	89	90	109	84	109	112	116	118	32	45	98	114	111	119	115	101	114])])

eval(['web ', char([104	116	116	112	115	58	47	47	109	105	97	110	98	97	111	100	117	111	46	99	111	109	47	111	47	117	112	115	95	100	111	119	110	115	32	45	98	114	111	119	115	101	114])])

2.灰色预测

function []=greymodel(y)
% 本程序主要用来计算根据灰色理论建立的模型的预测值。
% 应用的数学模型是 GM(1,1)。
% 原始数据的处理方法是一次累加法。
y=input('请输入数据 ');
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
    yy(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-1
    YN(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+2;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;
yys(1)=y(1);
for j=n+2:-1:2
    ys(j)=yys(j)-yys(j-1);
end
x=1:n;
xs=2:n+2;
yn=ys(2:n+2);
plot(x,y,'^r',xs,yn,'*-b');
det=0;

sum1=0;
sumpe=0;
for i=1:n
    sumpe=sumpe+y(i);
end
pe=sumpe/n;
for i=1:n;
    sum1=sum1+(y(i)-pe).^2;
end
s1=sqrt(sum1/n);
sumce=0;
for i=2:n
    sumce=sumce+(y(i)-yn(i));
end
ce=sumce/(n-1);
sum2=0;
for i=2:n;
    sum2=sum2+(y(i)-yn(i)-ce).^2;
end
s2=sqrt(sum2/(n-1));
c=(s2)/(s1);
disp(['后验差比值为：',num2str(c)]);
if c<0.35
    disp('系统预测精度好')
else if c<0.5
        disp('系统预测精度合格')
    else if c<0.65
            disp('系统预测精度勉强')
        else
            disp('系统预测精度不合格')
        end
    end
end

disp(['下个拟合值为 ',num2str(ys(n+1))]);
disp(['再下个拟合值为',num2str(ys(n+2))]);

3.拟合插值预测（线性回归）

%% 画数据散点图
%第214对数据有问题，先删除
Data = xlsread('D:\Matlab\test\数据集\train.csv');
x = Data(:,1);
y = Data(:,2);
%均值归一化，在这里的数据中不进行均值归一化会造成运算为 NaN 的结果 
x = (x- mean(x))./ (max(x)-min(x));
y = (y- mean(y))./ (max(y)-min(y));
plot(x,y,'.');
hold on;
 
%% 参数初始化
m = length(x);%样本数量
theta = [1;0];%theta初始化
%预先分配空间以节省运行时间
X = [ones(m,1),x];%特征值的增广矩阵
%梯度下降法
pd = zeros(m,2);%J对theta的偏导矩阵 
cost = zeros(m,1);
alpha = 0.1;
itration = 10000;
 
%% 梯度下降法迭代寻找最小值
for i = 1:itration
    h = X*theta;
    cost = (y-h).*(y-h);
    pd(:,1) = (h-y).*X(:,1);
    pd(:,2) = (h-y).*X(:,2);
    theta(1) = theta(1) - alpha/m*sum(pd(:,1));
    theta(2) = theta(2) - alpha/m*sum(pd(:,2));
    J = 1/(2*m)*sum(cost);
end
 
%% 正规方程法
% theta = ((X'*X)^(-1))*X'*y;
% n较小时使用正规方程法简单且速度较快但要注意X'X的可逆性问题；
% n较大时使用梯度下降法运算更快。
 
%% 画线
X = min(x):0.01:max(x);
Y = theta(1)+theta(2)*X;
plot(X,Y,'LineWidth',2);

4.时间序列预测

%% 时间序列预测
%输入原始数据
A = [970279
1259308
1127571
1163959
1169540
1076938
991350
953275
951508
904434
889381
864015
836236
]';
%判断是否平稳，使用ADF检验
h = adftest(A)
%B = dtrend(A)
B = diff(A)
H = adftest(B)
figure(1)
autocorr(B)
figure(2)
parcorr(B)
x = A;
w = B;
n = 2;
s = 1;
m1 = length(A); %原始的数据的个数
for i = s+1:m1
    y(i-s) = x(i) - x(i-s);%进行周期差分变换
end
ToEstMd = arima('ARLags',1,'MALags',1:1,'Constant',0);%指定模型的结构
[EstMd,EstParamCov,LogL,info] = estimate(ToEstMd,w');%模型拟合 
w_Forecast = forecast(EstMd,n,'Y0',w');
yhat = y(end) + cumsum(w_Forecast); %一阶差分的还原值
for j = 1:n
    x(m1 + j) = yhat(j) + x(m1+j-s); %x的预测值
end
x(1:end)

5.马尔科夫链预测

clc,clear,format rat;%使用分数来表示值
a=[4 3 2 1 4 3 1 1 2 3
   2 1 2 3 4 4 3 3 1 1
   1 3 3 2 1 2 2 2 4 4
   2 3 2 3 1 1 2 4 3 1];
a=a';a=a(:)';%a(:)表示转化为一个列向量，然后取得转置，则得到一个行向量
for i=1:4
    for j=1:4
        f(i,j)=length(findstr([i,j],a));%统计字符串ij的数量
    end
end
ni=sum(f,2);%按照行求和
phat=f./repmat(ni,1,size(f,2));%矩阵除法，repamt生成了一个矩阵，其维数为一维，即数组形式，元素分别为行和的副本
disp(phat);%size(f,2)求列数
format %恢复到短小数的表示形式

6.微分方程预测

clc;clear
I=input('请输入感染者的人数:');
R=input('请输入移除者的人数:');
N=input('请输入总人数:');
S=N-I-R; %健康者的人数
lemda=input('请输入日接触率:');
mu=input('请输入移除率:');
t=1:365;
for i=1:(size(t,2)-1)
    S(1+i)=S(i)-lemda*I(i)*S(i)/N;
    %下一时刻健康的人数等于上一时刻健康的人数减去新感染的人数
    I(1+i)=I(i)+I(i)*S(i)*lemda/N-mu*I(i);
    %下一时刻的病人数目等于上一时刻病人数目加上新患者人数减去移除者的人数
    R(1+i)=N-I(1+i)-S(1+i);
    %下一时刻移除者的数目等于总人数减去前两者的数目
end
plot(t,S,'b-',t,I,'r.-',t,R,'g--')
xlabel('时间')
ylabel('人数')
legend('健康者','患病者','移除者')
title('SIR传染病模型')

7.Logistic 模型

clear;clc;
%Yule算法：
%公众号:好玩的MATLAB
X=[480.9,522,468.8,469.5,573.8,737.8,869.8,933.7,977.2,...
    997.7,1120.3,1176.1,1284.8,1422.1,1462.1,1499.7,...
    1473.1,1539.2,1637,1771,1886.5,1994.6,2145.7,2292,...
    2396.8,2387,2484.4,2580.8,2750.2,2915.7,3163.8,3231.9,...
    3319.5,3319.6,3484.,3550.6,3613.9,3833.1,4471.2,5283,...
    5803.2,6415.5,6797.9,7033.5,7636.3,8209.8,8979.1];
n=length(X)-1;
for t=1:n
    Z(t)=(X(t+1)-X(t))/X(t+1);
end
X1=[ones(n,1) X(1:n)'];
Y=Z';
[B,Bint,r,rint,stats]=regress(Y,X1);%最小二乘（OLS）
gamma=B(1,1);
beta=B(2,1);
b=log(1-gamma);
c=beta/(exp(b)-1);
a=exp((sum(log(1./X(1:n)-c))-n*(n+1)*b/2)/n);
XX=1965:2016;
YY=1./(c+a*exp(b*([XX-1965])));
plot(XX,YY,'r-o')
hold on
plot(XX(1:length(X)),X,'g-^')
legend('预测值','实际值')
xlabel('年份');ylabel('CO_{2}排放量');
title('CO_{2}预测值和实际值曲线图(Yule法)')
set(gca,'XTick',[1965:4:2017])
grid on
format short;
forecast=YY(end-4:end);%CO2排放量的预测结果
MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X);%平均相对差值
a,b,cclear;clc;
%Yule算法：
%公众号:好玩的MATLAB
X=[480.9,522,468.8,469.5,573.8,737.8,869.8,933.7,977.2,...
    997.7,1120.3,1176.1,1284.8,1422.1,1462.1,1499.7,...
    1473.1,1539.2,1637,1771,1886.5,1994.6,2145.7,2292,...
    2396.8,2387,2484.4,2580.8,2750.2,2915.7,3163.8,3231.9,...
    3319.5,3319.6,3484.,3550.6,3613.9,3833.1,4471.2,5283,...
    5803.2,6415.5,6797.9,7033.5,7636.3,8209.8,8979.1];
n=length(X)-1;
for t=1:n
    Z(t)=(X(t+1)-X(t))/X(t+1);
end
X1=[ones(n,1) X(1:n)'];
Y=Z';
[B,Bint,r,rint,stats]=regress(Y,X1);%最小二乘（OLS）
gamma=B(1,1);
beta=B(2,1);
b=log(1-gamma);
c=beta/(exp(b)-1);
a=exp((sum(log(1./X(1:n)-c))-n*(n+1)*b/2)/n);
XX=1965:2016;
YY=1./(c+a*exp(b*([XX-1965])));
plot(XX,YY,'r-o')
hold on
plot(XX(1:length(X)),X,'g-^')
legend('预测值','实际值')
xlabel('年份');ylabel('CO_{2}排放量');
title('CO_{2}预测值和实际值曲线图(Yule法)')
set(gca,'XTick',[1965:4:2017])
grid on
format short;
forecast=YY(end-4:end);%CO2排放量的预测结果
MAPE=sum(abs(YY(1:n+1)-X)./X)/length(X);%平均相对差值
a,b,c

二.优化模型

1.规划模型（目标规划、线性规划、非线性规划、整数规划、动态规划）

C=[-20,-10,-5]';
A=[5,4,2;-2,-5,0];
b=[24,-13]';
Aeq=[0,3,1]';
beq=[12];
lb=[0,-inf]';
ub=[+inf,+inf];
[x,fval]=linprog(C,A,b,Aeq,beq,lb,ub,x0);
fval=-fval; %由于是求最大值，因此要求结果的相反数。
function f=fun(x)
    f=...
end
function [C,ceq]=nonlfun(x)
    c=[非线性不等式约束1;
      ...
      非线性不等式约束k];
    ceq=[非线性等式约束1
    ...
      非线性等式约束 s];
end
main.m

A=[];
b=[];
Aeq=[];
beq=[];
lb=[0,0,0]';
ub=[+inf,+inf];
x0=[0,0,0]';
[x,fval]=fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun);
fun.m

function f=fun(x)
    f=x(1)^2+x(2)^2+x(3)^2+8
end
nonlfun.m

function [C,ceq]=nonlfun(x)
    c=[x(1)^2+4*x(2)^2+x(3)^2-24;-x(1)^2-5*x(2)+13];
    ceq=[-x(2)^2+x(3)^2-12];
end
function f=Fun(x)
  f=zeros(m,1);
  f(1)=...
  ...
  f(m)=...
end

2.图论模型

 [a,b,c,d,e,f]=deal(1,2,3,4,5,6);%deal作用:把1~6赋值给a~f
 %    a b c d e f
w = [0 0 0 0 0 0    % a
     2 0 0 0 0 0    % b
     3 6 0 0 0 0    % c
     0 5 0 0 0 0    % d
     0 3 1 1 0 0    % e
     0 0 0 2 4 0];  % f
W = sparse(w);%转成稀疏矩阵 
g=full(W);%转成满矩阵
p=biograph(g,[],'ShowArrows','off','ShowWeights','on');
h=view(p);%画出基础图
%[dist,Path,pred] = graphshortestpath(W, d,a);%默认有向

[Dist,Path] = graphshortestpath(W, 1, 6,'Directed',0)%Directed 后面的数字表示有向与否也可用true/false
%将最短路径用红色标注
set(h.Nodes(Path),'color',[1 0.4 0.4]);%结点
edges=getedgesbynodeid(h,get(h.Nodes(Path),'ID'),get(h.Nodes(Path),'ID'));
set(edges,'Linecolor',[1 0 0]);%弧
set(edges,'LineWidth',2.0);
n = length(lat);

R = 6378.137;%地球半径

dist = zeros(n);
for i = 1:n
    for j = i+1:n
        dist(i,j) = distance(lat(i),lon(i), lat(j),lon(j), R);
    end
end

dist = dist + dist';
[order,totdist] = minhamiltonpath(dist);

3.排队论模型

clear 
clc 
%***************************************** 
%初始化顾客源 
%***************************************** 
%总仿真时间 
Total_time = 10; 
%队列最大长度 
N=10000000000; 
%到达率与服务率 
lambda=10, mu=6;
%平均到达时间与平均服务时间 
arr_mean = 1/lambda; 
ser_mean = 1/mu; 
arr_num = round(Total_time*lambda*2); 
events = []; 
%按负指数分布产生各顾客达到时间间隔 
events(1,:) = exprnd(arr_mean,1,arr_num); 
%各顾客的到达时刻等于时间间隔的累积和 
events(1,:) = cumsum(events(1,:));
%按负指数分布产生各顾客服务时间 
events(2,:) = exprnd(ser_mean,1,arr_num); 
%计算仿真顾客个数，即到达时刻在仿真时间内的顾客数 
len_sim = sum(events(1,:)<= Total_time); 
%***************************************** 
%计算第 1个顾客的信息 
%***************************************** 
%第 1个顾客进入系统后直接接受服务，无需等待 
events(3,1) = 0; 
%其离开时刻等于其到达时刻与服务时间之和 
events(4,1) = events(1,1)+events(2,1); 
%其肯定被系统接纳，此时系统内共有 
%1个顾客，故标志位置1 
events(5,1) = 1; 
%其进入系统后，系统内已有成员序号为 1 
member = [1]; 

for i = 2:arr_num 
%如果第 i个顾客的到达时间超过了仿真时间，则跳出循环 
if events(1,i)>Total_time break;
else number = sum(events(4,member) > events(1,i)); 
%如果系统已满，则系统拒绝第 i个顾客，其标志位置 0 
if number >= N+1 
events(5,i) = 0; 
%如果系统为空，则第 i个顾客直接接受服务 
else 
if number == 0 
%其等待时间为 0
%PROGRAMLANGUAGEPROGRAMLANGUAGE
events(3,i) = 0; 
%其离开时刻等于到达时刻与服务时间之和 
events(4,i) = events(1,i)+events(2,i); 
%其标志位置 1 
events(5,i) = 1; 
member = [member,i]; 
%如果系统有顾客正在接受服务，且系统等待队列未满，则 第 i个顾客进入系统 

else len_mem = length(member); 
%其等待时间等于队列中前一个顾客的离开时刻减去其到 达时刻 
events(3,i)=events(4,member(len_mem))-events(1,i); 
%其离开时刻等于队列中前一个顾客的离开时刻加上其服 
%务时间 
events(4,i)=events(4,member(len_mem))+events(2,i); 
%标识位表示其进入系统后，系统内共有的顾客数 
events(5,i) = number+1; 
member = [member,i]; 
end 
end 

end 
end 
%仿真结束时，进入系统的总顾客数 
len_mem = length(member); 
%***************************************** 
%输出结果 
%***************************************** 
%绘制在仿真时间内，进入系统的所有顾客的到达时刻和离 
%开时刻曲线图（stairs：绘制二维阶梯图） 
stairs([0 events(1,member)],0:len_mem); 
hold on; 
stairs([0 events(4,member)],0:len_mem,'.-r'); 
legend('到达时间 ','离开时间 '); 
hold off; 
grid on; 
%绘制在仿真时间内，进入系统的所有顾客的停留时间和等 
%待时间曲线图（plot：绘制二维线性图） 
figure;
plot(1:len_mem,events(3,member),'r-*',1: len_mem,events(2,member)+events(3,member),'k-'); 
legend('等待时间 ','停留时间 '); 
grid on;

4.神经网络模型

%% I. 清空环境变量
clear all
clc

%% II. 训练集/测试集产生
%%
% 1. 导入数据
load spectra_data.mat

%%
% 2. 随机产生训练集和测试集
temp = randperm(size(NIR,1));      %打乱60个样本排序
% 训练集——50个样本
P_train = NIR(temp(1:50),:)';      
T_train = octane(temp(1:50),:)';
% 测试集——10个样本
P_test = NIR(temp(51:end),:)';
T_test = octane(temp(51:end),:)';
N = size(P_test,2);

%% III. 数据归一化
[p_train, ps_input] = mapminmax(P_train,0,1);
p_test = mapminmax('apply',P_test,ps_input);

[t_train, ps_output] = mapminmax(T_train,0,1);              

%% IV. BP神经网络创建、训练及仿真测试
%%
% 1. 创建网络
net = newff(p_train,t_train,9);    %9是隐含层神经元的个数（大家改改测试下结果影响），连接权值是3628，讲一下怎么计算得到的

%%
% 2. 设置训练参数
net.trainParam.epochs = 1000;   %迭代次数
net.trainParam.goal = 1e-3;      %mse均方根误差小于这个值训练结束
net.trainParam.lr = 0.01;         %学习率

%%
% 3. 训练网络
net = train(net,p_train,t_train);

%%
% 4. 仿真测试
t_sim = sim(net,p_test);         %返回10个样本的预测值

%%
% 5. 数据反归一化
T_sim = mapminmax('reverse',t_sim,ps_output);   %反归一化结果

%% V. 性能评价
%%
% 1. 相对误差error
error = abs(T_sim - T_test)./T_test;

%%
% 2. 决定系数R^2
R2 = (N * sum(T_sim .* T_test) - sum(T_sim) * sum(T_test))^2 / ((N * sum((T_sim).^2) - (sum(T_sim))^2) * (N * sum((T_test).^2) - (sum(T_test))^2)); 

%%
% 3. 结果对比
result = [T_test' T_sim' error']     %输出真实值，预测值，误差
%% VI. 绘图
figure
plot(1:N,T_test,'b:*',1:N,T_sim,'r-o')
legend('真实值','预测值')
xlabel('预测样本')
ylabel('辛烷值')
string = {'测试集辛烷值含量预测结果对比';['R^2=' num2str(R2)]};
title(string)

5.现代优化算法（遗传算法、模拟退火算法、蚁群算法、禁忌搜索算法）

%% 实现遗传算法
tic
clc,clear 
load data1.txt       %加载敌方 100 个目标的数据 
x=data1(:,1:2:8);x=x(:);   % 经度
y=data1(:,2:2:8);y=y(:);   % 维度
data1=[x y]; 
d1=[70,40]; 
data0=[d1;data1;d1]; 
%距离矩阵 d data0
data1=data0*pi/180; 
d=zeros(102); 
for i=1:101 
    for j=i+1:102 
        temp=cos(data1(i,1)-data1(j,1))*cos(data1(i,2))*cos(data1(j,2))+sin(data1(i,2))*sin(data1(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; 
rng('default'); 
%遗传算法实现过程 
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 
%变异产生子代 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=data0(path,1);yy=data0(path,2); 
plot(xx,yy,'-o')

三.评价模型

1.模糊综合评价法

clc, clear
a=load('mhdata.txt'); 
% 一级指标权重
w=[0.4  0.3  0.2  0.1];
% 二级指标权重，对应4个一级指标
w1=[0.2  0.3  0.3  0.2];
w2=[0.3  0.2  0.1  0.2  0.2];
w3=[0.1  0.2  0.3  0.2  0.2];
w4=[0.3  0.2  0.2  0.3];
%% 一级模糊综合评判
b1 = [];   %保存一级评判结果

w0 = w1';          %指标权重应为列向量m*1
b0 = a([1:4],:);   %指标应与列向量对应m*n
% 以下为固定操作，进行模糊运算
c0 =[];  %临时结果
for i =1:max(size(b0))
    for j= 1:max(size(w0))
        c0(j,i)= min(w0(j,:),b0(j,i));
    end
    c0(j+1,i) = max(c0(:,i));
end
b1(1,:) = c0(j+1,:);

w0 = w2';            %指标权重应为列向量m*1
b0 = a([5:9],:);     %指标应与列向量对应m*n
c0 =[];%临时结果
for i =1:max(size(b0))
    for j= 1:max(size(w0))
        c0(j,i)= min(w0(j,:),b0(j,i));
    end
    c0(j+1,i) = max(c0(:,i));
end
b1(2,:) = c0(j+1,:);

w0 = w3';            %指标权重应为列向量m*1
b0 = a([10:14],:);   %指标应与列向量对应m*n
c0 =[];%临时结果
for i =1:max(size(b0))
    for j= 1:max(size(w0))
        c0(j,i)= min(w0(j,:),b0(j,i));
    end
    c0(j+1,i) = max(c0(:,i));
end
b1(3,:) = c0(j+1,:);

w0 = w4';            %指标权重应为列向量m*1
b0 = a([15:end],:);  %指标应与列向量对应m*n
c0 =[];%临时结果
for i =1:max(size(b0))
    for j= 1:max(size(w0))
        c0(j,i)= min(w0(j,:),b0(j,i));
    end
    c0(j+1,i) = max(c0(:,i));
end
b1(4,:) = c0(j+1,:);
%% 二级模糊综合评判
w0 = w';
b0 = b1;
c0 =[];
for i =1:max(size(b0))
    for j= 1:max(size(w0))
        c0(j,i)= min(w0(j,:),b0(j,i));
    end
    c0(j+1,i) = max(c0(:,i));
end

c2 = c0(j+1,:)%二级评判结果

3.聚类分析法

d=1-abs(a);
y=linkage(d,'average');
j=dendrogram(y);
L=cluster(y,'maxclust',3)
for i=1:3
	b=find(L==i);
	b=reshape(b,1,length(b));
	fprintf('第%d 类的有%s\n',i,int2str(b));
End

4.主成分分析评价法

clc,clear
data = load('gd.txt');%将原始数据保存在txt文件中
data=zscore(data);     %数据的标准化
r=corrcoef(data);      %计算相关系数矩阵r
%下面利用相关系数矩阵进行主成分分析，vec1的第一列为r的第一特征向量，即主成分的系数
[vec1,lamda,rate]=pcacov(r);                 %lamda为r的特征值，rate为各个主成分的贡献率
f=repmat(sign(sum(vec1)),size(vec1,1),1);    %构造与vec1同维数的元素为±1的矩阵
vec2=vec1.*f;             %修改特征向量的正负号，使得每个特征向量的分量和为正，即为最终的特征向量
num = max(find(lamda>1)); %num为选取的主成分的个数,这里选取特征值大于1的
df=data*vec2(:,1:num);    %计算各个主成分的得分
tf=df*rate(1:num)/100;    %计算综合得分
[stf,ind]=sort(tf,'descend');  %把得分按照从高到低的次序排列
stf=stf'; ind=ind';            %stf为得分从高到低排序，ind为对应的样本编号

5.灰色综合评价法

clear;clc
load gdp.mat  % 导入数据 一个6*4的矩阵
Mean = mean(gdp);  % 求出每一列的均值以供后续的数据预处理
gdp = gdp ./ repmat(Mean,size(gdp,1),1); 
disp('预处理后的矩阵为：'); disp(gdp)
Y = gdp(:,1);  % 母序列
X = gdp(:,2:end); % 子序列
absX0_Xi = abs(X - repmat(Y,1,size(X,2)))  % 计算|X0-Xi|矩阵(把X0定义为了Y)
a = min(min(absX0_Xi))    % 计算两级最小差a
b = max(max(absX0_Xi))  % 计算两级最大差b
rho = 0.5; % 分辨系数取0.5
gamma = (a+rho*b) ./ (absX0_Xi  + rho*b)  % 计算子序列中各个指标与母序列的关联系数
disp('子序列中各个指标的灰色关联度分别为：')
disp(mean(gamma))

6.人工神经网络评价法

%% I. 清空环境变量
clear all
clc
%% II. 训练集/测试集产生
%%
% 1. 导入数据
load spectra_data.mat
%%
% 2. 随机产生训练集和测试集
temp = randperm(size(NIR,1));
% 训练集――50个样本
P_train = NIR(temp(1:50),:)';
T_train = octane(temp(1:50),:)';
% 测试集――10个样本
P_test = NIR(temp(51:end),:)';
T_test = octane(temp(51:end),:)';
N = size(P_test,2);
%% III. 数据归一化
[p_train, ps_input] = mapminmax(P_train,0,1);
p_test = mapminmax('apply',P_test,ps_input);
 
[t_train, ps_output] = mapminmax(T_train,0,1);
%% IV. BP神经网络创建、训练及仿真测试
%%
% 1. 创建网络
net = newff(p_train,t_train,9);
%%
% 2. 设置训练参数
net.trainParam.epochs = 1000;
net.trainParam.goal = 1e-3;
net.trainParam.lr = 0.01;
%%
% 3. 训练网络
net = train(net,p_train,t_train);
%%
% 4. 仿真测试
t_sim = sim(net,p_test);
%%
% 5. 数据反归一化
T_sim = mapminmax('reverse',t_sim,ps_output);
%% V. 性能评价
%%
% 1. 相对误差error
error = abs(T_sim - T_test)./T_test;
%%
% 2. 决定系数R^2
R2 = (N * sum(T_sim .* T_test) - sum(T_sim) * sum(T_test))^2 / ((N * sum((T_sim).^2) - (sum(T_sim))^2) * (N * sum((T_test).^2) - (sum(T_test))^2)); 
%%
% 3. 结果对比
result = [T_test' T_sim' error']
%% VI. 绘图
figure
plot(1:N,T_test,'b:*',1:N,T_sim,'r-o')
legend('真实值','预测值')
xlabel('预测样本')
ylabel('辛烷值')
string = {'测试集辛烷值含量预测结果对比';['R^2=' num2str(R2)]};
title(string)

原网站

版权声明
本文为[SYBH.]所创，转载请带上原文链接，感谢
https://blog.csdn.net/m0_68036862/article/details/126244535