在MATLAB中实现Kriging代理模型

%% 清除环境变量
clear
clc
clf

%% 2.读取数据
data=xlsread('.xlsx','A1:D201'); % Matlab2021版本以上无法使用xlsread函数，可用Load函数替代 

% 设置神经网络的输入和输出
input=data(:,1:end-1);    %第1列至倒数第2列为输入
output=data(:,end);       %最后1列为输出

num_size = 0.8;                              % 训练集占数据集比例
outdim = 1;                                  % 最后一列为输出
num_samples = size(data, 1);                  % 样本个数
data = data(randperm(num_samples), :);         % 打乱数据集（不希望打乱时，注释该行）
num_train_s = round(num_size * num_samples); % 训练集样本个数
f_ = size(data, 2) - outdim;                  % 输入特征维度
%% 3.设置训练集和测试集
X1 = data(1: num_train_s, 1: f_)';
X1 = X1';
Y1 = data(1: num_train_s, f_ + 1: end)';
Y1 = Y1';
M = size(X1, 2);

X2 = data(num_train_s + 1: end, 1: f_)';
X2 = X2';
Y2 = data(num_train_s + 1: end, f_ + 1: end)';
Y2 = Y2';
N = size(X2, 2);

theta = [10.*ones(1,3)];%初始参数值，参数个数为自变量数
lob = [1E-5.*ones(1,3)];%参数下界
upb = [20.*ones(1,3)];%参数上界


%% 克里金模型训练（regpoly1是1阶多项式，regpoly0是0阶，regpoly2是2阶）
[dmodel,perf] = dacefit(X1,Y1, @regpoly2, @corrgauss, theta, lob, upb);

%% 预测和评价指标
[P1] = predictor(X1,dmodel);%训练结果
[P2] = predictor(X2,dmodel);%测试结果

R1 = round(max(0,1-(sum((P1-Y1).^2)/sum((Y1-mean(Y1)).^2))),4);
RMSE1 = round(sqrt(mean((P1-Y1).^2)),4);

R2 = round(max(0,1-(sum((P2-Y2).^2)/sum((Y2-mean(Y2)).^2))),4);
RMSE2 = round(sqrt(mean((P2-Y2).^2)),4);

%%  训练集绘图
set(figure,'Position',[100,100,1200,500]);
subplot(1,2,1)
hold on
plot(Y1,'b-*')
plot(P1,'r-o')
legend('真实值','预测值')
xlim([1,length(Y1)])
title(["模型训练集结果"    "R^2="+num2str(R1)+"，RMSE="+num2str(RMSE1)])
subplot(1,2,2)
err=abs(Y1-P1);
plot(err)
legend('绝对误差')
%% 测试集绘图
set(figure,'Position',[100,100,1200,500]);
subplot(1,2,1)
hold on
plot(Y2,'b-*')
plot(P2,'r-o')
legend('真实值','预测值')
xlim([1,length(Y2)])
title(["模型测试集结果"    "R^2="+num2str(R2)+"，RMSE="+num2str(RMSE2)])
subplot(1,2,2)
err=abs(Y2-P2);
plot(err)
legend('绝对误差')