clc;clear;close all;
data_select = xlsread('多元回归数据集.xlsx');
%% 数据划分
x_feature_label=data_select(:,1:end-1); %x特征
y_feature_label=data_select(:,end); %y标签
index_label1=1:(size(x_feature_label,1));
index_label=randperm(103); % 数据索引 ; % 数据索引
if isempty(index_label)
index_label=index_label1;
end
spilt_ri=[8 1 1]; %划分比例 训练集:验证集:测试集
train_num=round(spilt_ri(1)/(sum(spilt_ri))*size(x_feature_label,1)); %训练集个数
vaild_num=round((spilt_ri(1)+spilt_ri(2))/(sum(spilt_ri))*size(x_feature_label,1)); %验证集个数
%训练集，验证集，测试集
train_x_feature_label=x_feature_label(index_label(1:train_num),:);
train_y_feature_label=y_feature_label(index_label(1:train_num),:);
vaild_x_feature_label=x_feature_label(index_label(train_num+1:vaild_num),:);
vaild_y_feature_label=y_feature_label(index_label(train_num+1:vaild_num),:);
test_x_feature_label=x_feature_label(index_label(vaild_num+1:end),:);
test_y_feature_label=y_feature_label(index_label(vaild_num+1:end),:);
%Zscore 标准化
%训练集
x_mu = mean(train_x_feature_label); x_sig = std(train_x_feature_label);
train_x_feature_label_norm = (train_x_feature_label - x_mu) ./ x_sig; % 训练数据标准化
y_mu = mean(train_y_feature_label); y_sig = std(train_y_feature_label);
train_y_feature_label_norm = (train_y_feature_label - y_mu) ./ y_sig; % 训练数据标准化
%验证集
vaild_x_feature_label_norm = (vaild_x_feature_label - x_mu) ./ x_sig; %验证数据标准化
vaild_y_feature_label_norm=(vaild_y_feature_label - y_mu) ./ y_sig; %验证数据标准化
%测试集
test_x_feature_label_norm = (test_x_feature_label - x_mu) ./ x_sig; % 测试数据标准化
test_y_feature_label_norm = (test_y_feature_label - y_mu) ./ y_sig; % 训练数据标准化
%% 算法处理块
disp('MLP回归')
t1=clock;
hidden_size=[15 9]; %神经网络隐藏层
num_iter=10; %迭代次数
Mdl=fitrnet(train_x_feature_label_norm,train_y_feature_label_norm,'OptimizeHyperparameters','auto','HyperparameterOptimizationOptions',struct('AcquisitionFunctionName','expected-improvement-plus','MaxObjectiveEvaluations',num_iter)) ;
y_train_predict_norm=predict(Mdl,train_x_feature_label_norm); %训练集预测结果
y_vaild_predict_norm=predict(Mdl,vaild_x_feature_label_norm); %验证集预测结果
y_test_predict_norm=predict(Mdl,test_x_feature_label_norm); %测试集预测结果
t2=clock;
Time=t2(3)*3600*24+t2(4)*3600+t2(5)*60+t2(6)-(t1(3)*3600*24+t1(4)*3600+t1(5)*60+t1(6));
y_train_predict=y_train_predict_norm*y_sig+y_mu; %反标准化操作
y_vaild_predict=y_vaild_predict_norm*y_sig+y_mu;
y_test_predict=y_test_predict_norm*y_sig+y_mu;
train_y=train_y_feature_label;
train_MAE=sum(abs(y_train_predict-train_y))/length(train_y) ; disp(['训练集平均绝对误差MAE：',num2str(train_MAE)])
train_MAPE=sum(abs((y_train_predict-train_y)./train_y))/length(train_y); disp(['训练集平均相对误差MAPE：',num2str(train_MAPE)])
train_MSE=(sum(((y_train_predict-train_y)).^2)/length(train_y)); disp(['训练集均方根误差MSE：',num2str(train_MSE)])
train_RMSE=sqrt(sum(((y_train_predict-train_y)).^2)/length(train_y)); disp(['训练集均方根误差RMSE：',num2str(train_RMSE)])
train_R2= 1 - (norm(train_y - y_train_predict)^2 / norm(train_y - mean(train_y))^2); disp(['训练集均方根误差R2：',num2str(train_R2)])
vaild_y=vaild_y_feature_label;
vaild_MAE=sum(abs(y_vaild_predict-vaild_y))/length(vaild_y) ; disp(['验证集平均绝对误差MAE：',num2str(vaild_MAE)])
vaild_MAPE=sum(abs((y_vaild_predict-vaild_y)./vaild_y))/length(vaild_y); disp(['验证集平均相对误差MAPE：',num2str(vaild_MAPE)])
vaild_MSE=(sum(((y_vaild_predict-vaild_y)).^2)/length(vaild_y)); disp(['验证集均方根误差MSE：',num2str(vaild_MSE)])
vaild_RMSE=sqrt(sum(((y_vaild_predict-vaild_y)).^2)/length(vaild_y)); disp(['验证集均方根误差RMSE：',num2str(vaild_RMSE)])
vaild_R2= 1 - (norm(vaild_y - y_vaild_predict)^2 / norm(vaild_y - mean(vaild_y))^2); disp(['验证集均方根误差R2: ',num2str(vaild_R2)])
test_y=test_y_feature_label;test_y(find(test_y==0))=0.0001; %防止分母为0
test_MAE=sum(abs(y_test_predict-test_y))/length(test_y) ; disp(['测试集平均绝对误差MAE：',num2str(test_MAE)])
test_MAPE=sum(abs((y_test_predict-test_y)./test_y))/length(test_y); disp(['测试集平均相对误差MAPE：',num2str(test_MAPE)])
test_MSE=(sum(((y_test_predict-test_y)).^2)/length(test_y)); disp(['测试集均方根误差MSE：',num2str(test_MSE)])
test_RMSE=sqrt(sum(((y_test_predict-test_y)).^2)/length(test_y)); disp(['测试集均方根误差RMSE：',num2str(test_RMSE)])
test_R2= 1 - (norm(test_y - y_test_predict)^2 / norm(test_y - mean(test_y))^2); disp(['测试集均方根误差R2：',num2str(test_R2)])
disp(['算法运行时间Time: ',num2str(Time)])
%% 绘图块
figure
XX=1:length(test_y_feature_label);
index_show=1:50;
plot(gca,XX,test_y_feature_label, 'r-*', 'LineWidth', 1)
hold (gca,'on')
plot(gca, XX,y_test_predict, 'b-o', 'LineWidth', 1)
hold (gca,'on')
set(gca,'LineWidth',1)
xlabel('样本编号')
ylabel('回归结果')
title(gca,'贝叶斯优化MLP')
legend('真实值','预测值')
grid on
