tic;%计时开始
clear all
close all
clc
format compact
%% 加载数据
% num1 = xlsread('a.xlsx',1,'A4:C68');
% num2 = xlsread('a.xlsx',1,'D4:F21');
% num3 = [num1(:,1:2);num2(:,1:2)];
% num4 = [num1(:,3);num2(:,3)];
num = xlsread('数据2.xls');
num2 = [];
for ii = 14500:15000
num2 = [num2;num(ii:ii+191,3)'];
end
num3 = num2(:,1:96);
num4 = num2(:,97:192);
%% 归一化
[input,inputps]=mapminmax(num3',0,1);
[output,outputns]=mapminmax(num4',0,1);
input=input';
output=output';
%% 划分数据集
% n = randperm(81);
n =1:200;
P=input((1:end-1),:); %训练输入
T=output((1:end-1),:);
P_test=input(end,:);%测试输入
T_test=output(end,:);
clear data m n input output
%% 训练样本构造，分块，批量
numcases=5;%每块数据集的样本个数
numdims=size(P,2);%单个样本的大小
numbatches=9;%将1000组训练样本，分成250批，每一批5组
% 训练数据
for i=1:numbatches
train=P((i-1)*numcases+1:i*numcases,:);
batchdata(:,:,i)=train;
end%将分好的10组数据都放在batchdata中
%% 2.训练RBM
%% rbm参数
maxepoch=1000;%训练rbm的次数
numhid=130; numpen=200; numpen2=10;%dbn隐含层的节点数
disp('构建一个3层的深度置信网络DBN用于特征提取');
%% 无监督预训练
fprintf(1,'Pretraining Layer 1 with RBM: %d-%d ',numdims,numhid);
restart=1;
rbm1;%使用cd-k训练rbm，注意此rbm的可视层不是二值的，而隐含层是二值的
vishid1=vishid;hidrecbiases=hidbiases;
fprintf(1,'\nPretraining Layer 2 with RBM: %d-%d ',numhid,numpen);
batchdata=batchposhidprobs;%将第一个RBM的隐含层的输出作为第二个RBM 的输入
numhid=numpen;%将numpen的值赋给numhid，作为第二个rbm隐含层的节点数
restart=1;
rbm1;
hidpen=vishid; penrecbiases=hidbiases; hidgenbiases=visbiases;
fprintf(1,'\nPretraining Layer 3 with RBM: %d-%d\n ',numpen,numpen2);%200-100
batchdata=batchposhidprobs;%显然，将第二哥RBM的输出作为第三个RBM的输入
numhid=numpen2;%第三个隐含层的节点数
restart=1;
rbm1;
hidpen2=vishid; penrecbiases2=hidbiases; hidgenbiases2=visbiases;
%%%% 将预训练好的RBM用于初始化DBN权重%%%%%%%%%
w1=[vishid1; hidrecbiases]; %
w2=[hidpen; penrecbiases]; %
w3=[hidpen2; penrecbiases2];%
%% 有监督回归层训练
%===========================训练过程=====================================%
%==========DBN无监督用于提取特征，需要加上有监督的回归层==================%
%由于含有偏执，所以实际数据应该包含一列全为1的数，即w0x0+w1x1+..+wnxn 其中x0为1的向量 w0为偏置b
N1 = size(P,1);
digitdata = [P ones(N1,1)];
w1probs = 1./(1 + exp(-digitdata*w1));
w1probs = [w1probs ones(N1,1)];
w2probs = 1./(1 + exp(-w1probs*w2));
w2probs = [w2probs ones(N1,1)];
w3probs = 1./(1 + exp(-w2probs*w3));
H = w3probs';
nn=size(T,2);
T=T';
lamda=inf;%正则化系数
OutputWeight=pinv(H'+1/lamda) *T';%加入正则化系数lamda，lamda=inf就是没有正则化
Y=(H' * OutputWeight)';
%%%%%%%%%% 计算训练误差，不重要，看看图就行
% 反归一化
T=(mapminmax('reverse',T,outputns));
Y=(mapminmax('reverse',Y,outputns));
figure
plot(T','-*')
hold on
plot(Y','o-')
legend('期望输出','实际输出')
firstline = '训练阶段';
secondline = '实际输出与理想输出的结果对照';
title({firstline;secondline},'Fontsize',12);
xlabel('训练样本数')
%===========================测试过程=====================================%
%=======================================================================%
N2 = size(P_test,1);
w1=[vishid1; hidrecbiases];
w2=[hidpen; penrecbiases];
w3=[hidpen2; penrecbiases2];
test = [P_test ones(N2,1)];
%激活函数是常用的sigmoid时，也可以换成其他函数，甚至每层的都可以不一样
w1probs = 1./(1 + exp(-test*w1));
w1probs = [w1probs ones(N2,1)];
w2probs = 1./(1 + exp(-w1probs*w2));
w2probs = [w2probs ones(N2,1)];
w3probs = 1./(1 + exp(-w2probs*w3));
H1 = w3probs';
TY=(H1' * OutputWeight)'; % TY: the actual output of the testing data
% 反归一化
T_test=(mapminmax('reverse',T_test',outputns));
TY=(mapminmax('reverse',TY,outputns));
%%%%%%%%%% 计算测试结果
result=[TY;T_test];
error=TY-T_test;
fprintf('测试集输出结果分析\n');
fprintf('均方误差MSE\n');
MSE=mse(error)
fprintf('平均绝对误差MAE\n');
% MAE=mean(abs(TY-mean(T_test)))
% fprintf('平均相对误差MRE\n');
% MRE=mean(abs(TY-mean(T_test))./T_test); %不能计算平均相对误差，因为相对误差=绝对误差/真实值，而你真实值里面有0
% fprintf('不能计算平均相对误差，因为相对误差=绝对误差/真实值，而你真实值里面有0\n');
fprintf('相关系数\n');
a=corrcoef(TY,T_test);%皮尔逊相关系数
corrcoeff=a(1,2)
t=toc %计时结束
figure
plot(T_test,'r-*')
hold on
plot(TY,'bo-')
firstline = 'DBN深度信念网络测试阶段';
secondline = '实际输出与理想输出的结果对照';
title({firstline;secondline},'Fontsize',12);
xlabel('类别')
ylabel('样本')
legend('期望输出','验证输出')
