%PCA_BP神經網絡 2008.12.8 
%采用先主成分分析取得主元，有12個變為8個。后用BP神經網絡預測。結果
%數據擬合
clf reset; 
close all; 
clear; 
clc; 
echo on;

%-------------------------------------------------------------------
load paperdata;
beishao3=paperdata(:,1:10);
%數據獲取，預處理。beishao3是一個355*12的矩陣。355個樣本。

%樣本數據歸一化
%beishao3_G = beishao3./(ones(size(beishao3,1),1)*sqrt(diag(beishao3'*beishao3))');
%法二
[m,n]=size(beishao3);
for i=1:n
    beishao3_G(:,i)=0.1+(beishao3(:,i)-min(beishao3(:,i)))/(max(beishao3(:,i))-min(beishao3(:,i)))*(0.9-0.1);
end
%beishao3_G=zscore(beishao3);
%-------------------------------------------------------------------
%對數據進行princomp主成分分析

% [coefs,scores,variances,t2] = princomp(beishao3_G(:,1:8));
% 
% zcf = beishao3_G(:,1:8)*coefs;
% 
% 
% %-------------------------------------------------------------------
% % 選擇的主成分占整個的百分比
% percent_explained = 100*variances/sum(variances);
% pareto(percent_explained)
% title('篩選出的主成分');
% xlabel('Principal Component');
% ylabel('Variance Explained (%)');
% legend('主成分個數','每一個百分比',1);

%%-------------------------------------------------------------------
%%選擇訓練+測試樣本
sample_data = beishao3_G(:,1:5);

%-------------------------------------------------------------------
% 利用BP網絡實現函數逼近
%NEWFF--生成一個新的前向神經網絡
%TRAIN---對BP網絡進行訓練
%SIM----對BP網絡進行仿真

%-------------------------------------------------------------------
% 選取訓練樣本加測試樣本 
%弱精品位

myTraindata = sample_data(1:250,:);%選取的測試樣本數據
P = myTraindata'; 
T = beishao3_G(1:250,9)';%T為訓練樣本弱精品位的目標輸出

myTestdata = sample_data(251:353,:);%選取的訓練樣本數據
p1 = myTestdata';
t1 = beishao3_G(251:353,9)';%t1為訓練樣本的目標輸出

%———————————————————————————————————————
%畫出訓練樣本數據及目標輸出
i=1:length(myTraindata);

%創建BP神經網絡
%net=newff(minmax(P),[6 14 1],{'purelin','tansig','purelin'},'trainoss');
net = newff(minmax(P),[10 1],{'tansig','purelin'},'trainscg');

%———————————————————————————————————————
% net.trainParam.epochs=2500;
% %Maximum number of epochs to train
% net.trainParam.goal=1e-4;
% %Performance goal
% net.trainParam.mu=0.005;
% %Marquardt adjustment parameter
% net.trainParam.mu_dec=0.1;
% %Decrease factor for mu
% net.trainParam.mu_inc=5;
% %Increase factor for mu
% net.trainParam.mu_max=1e-10;
% %Maximum value for mu
% net.trainParam.max_fail=10;
% %Maximum validation failures
% net.trainParam.mem_reduc=1;
% %Factor to use for memory/speed tradeoff
% net.trainParam.min_grad=1e-10;
% %Minimum performance gradient
% net.trainParam.show=100;
% %Epochs between displays (NaN for no displays)
% net.trainParam.time=inf;
% %Maximum time to train in seconds

%———————————————————————————————————————
%神經網絡參數設置
% %net.performFcn='sse';
% %net.trainFcn = 'trainlm';%用這個的結果比較好
%net.trainFcn = 'trainscg';
%net.trainFcn = 'trainbr'; 
net.trainParam.goal=1e-3;
net.trainParam.show=100;
net.trainParam.lr=0.9;
net.trainParam.epochs=3000;
net.trainParam.Mu_max=1e-12;
net.trainParam.min_grad=1e-10;
%%%%%%%%%%%%%%%%%%%%%%%%下面一句不懂的話可以參考helpnet.trainParam.mc里面有說明
net.trainParam.mc=0.9;          
net.trainParam.mem_reduc=1;

%———————————————————————————————————————
% 重新初始化網絡
net = init(net);            

%———————————————————————————————————————
% 對BP神經網絡進行訓練
[net,tr] = train(net,P,T);

%———————————————————————————————————————
% 對BP神經網絡進行仿真分析
Yn=sim(net,P);

%———————————————————————————————————————
% 恢復被歸一化的數據
%[Y]=MAPMINMAX(Yn); 

%———————————————————————————————————————
%計算均方誤差
e=sqrt((Yn-T).^2)';
ee = mean(e);
MSE1 =  (1/(300-1))*sum((Yn-T).^2)';

%———————————————————————————————————————
%神經網絡訓練后的輸出結果
%畫圖描繪仿真結果
figure;
subplot(2,1,1)
plot(i,T,'b',i,Yn,'r');
title('BP神經網絡的函數逼近結果');
xlabel('訓練樣本數據');
ylabel('訓練結果');
legend('訓練樣本數據','訓練結果',1);
subplot(2,1,2)
plot(i,e(i),'r');
title('訓練誤差曲線');
xlabel('訓練樣本個數');
ylabel('誤差');
legend('訓練樣本個數','誤差',1);


%———————————————————————————————————————
%用未經訓練的數據對訓練成功的BP網絡進行驗證測試
%對測試數據進行歸一化處理
%tramnmx：利用預先計算的最大和最小值對數據進行變換
%p1n=tramnmx(p1,minp,maxp);
%對變換后的測試數據進行仿真
y = sim(net,p1);

% 將仿真結果還原成原始數據
%[y]=postmnmx(yn); 

%———————————————————————————————————————
%畫圖
figure;
j=1:length(myTestdata);
subplot(2,1,1)
plot(j,t1,'b',j,y,'r');
title('BP神經網絡的函數逼近結果檢驗');
xlabel('測試樣本數據');
ylabel('測試結果');
legend('測試樣本數據','測試結果',1);
E(j)=sqrt((y(j)-t1(j)).^2)';
EE = mean(E);
MSE = (1/52)*sum((y-t1).^2)';
subplot(2,1,2)
plot(j,E(j),'r')
title('測試誤差曲線');
xlabel('測試樣本個數');
ylabel('誤差');
legend('測試樣本個數','測試誤差',1);