Matlab利用隨機(jī)森林(RF)算法實(shí)現(xiàn)回歸預(yù)測(cè)詳解

更新時(shí)間：2023年02月17日 14:10:34 作者：瘋狂學(xué)習(xí)GIS

這篇文章主要為大家詳細(xì)介紹了Matlab如何利用隨機(jī)森林（RF）算法實(shí)現(xiàn)回歸預(yù)測(cè)，以及自變量重要性排序的操作，感興趣的小伙伴可以了解一下

本文分為兩部分，首先是對(duì)代碼進(jìn)行分段、詳細(xì)講解，方便大家理解；隨后是完整代碼，方便大家自行嘗試。另外，關(guān)于基于MATLAB的神經(jīng)網(wǎng)絡(luò)（ANN）代碼與詳細(xì)解釋，我們將在后期博客中介紹。

1 分解代碼

1.1 最優(yōu)葉子節(jié)點(diǎn)數(shù)與樹數(shù)確定

首先，我們需要對(duì)RF對(duì)應(yīng)的葉子節(jié)點(diǎn)數(shù)與樹的數(shù)量加以擇優(yōu)選取。

%% Number of Leaves and Trees Optimization

for RFOptimizationNum=1:5
    
RFLeaf=[5,10,20,50,100,200,500];
col='rgbcmyk';
figure('Name','RF Leaves and Trees');
for i=1:length(RFLeaf)
    RFModel=TreeBagger(2000,Input,Output,'Method','R','OOBPrediction','On','MinLeafSize',RFLeaf(i));
    plot(oobError(RFModel),col(i));
    hold on
end
xlabel('Number of Grown Trees');
ylabel('Mean Squared Error') ;
LeafTreelgd=legend({'5' '10' '20' '50' '100' '200' '500'},'Location','NorthEast');
title(LeafTreelgd,'Number of Leaves');
hold off;

disp(RFOptimizationNum);
end

其中，RFOptimizationNum是為了多次循環(huán)，防止最優(yōu)結(jié)果受到隨機(jī)干擾；大家如果不需要，可以將這句話刪除。

RFLeaf定義初始的葉子節(jié)點(diǎn)個(gè)數(shù)，我這里設(shè)置了從5到500，也就是從5到500這個(gè)范圍內(nèi)找到最優(yōu)葉子節(jié)點(diǎn)個(gè)數(shù)。

Input與Output分別是我的輸入（自變量）與輸出（因變量），大家自己設(shè)置即可。

運(yùn)行后得到下圖。

首先，我們看到MSE最低的線是紅色的，也就是5左右的葉子節(jié)點(diǎn)數(shù)比較合適；再看各個(gè)線段大概到100左右就不再下降，那么樹的個(gè)數(shù)就是100比較合適。

1.2 循環(huán)準(zhǔn)備

由于機(jī)器學(xué)習(xí)往往需要多次執(zhí)行，我們就在此先定義循環(huán)。

%% Cycle Preparation

RFScheduleBar=waitbar(0,'Random Forest is Solving...');
RFRMSEMatrix=[];
RFrAllMatrix=[];
RFRunNumSet=10;
for RFCycleRun=1:RFRunNumSet

其中，RFRMSEMatrix與RFrAllMatrix分別用來存放每一次運(yùn)行的RMSE、r結(jié)果，RFRunNumSet是循環(huán)次數(shù)，也就是RF運(yùn)行的次數(shù)。

1.3 數(shù)據(jù)劃分

接下來，我們需要將數(shù)據(jù)劃分為訓(xùn)練集與測(cè)試集。這里要注意：RF其實(shí)一般并不需要?jiǎng)澐钟?xùn)練集與測(cè)試集，因?yàn)槠淇梢圆捎么庹`差（Out of Bag Error，OOB Error）來衡量自身的性能。但是因?yàn)槲沂亲隽硕喾N機(jī)器學(xué)習(xí)方法的對(duì)比，需要固定訓(xùn)練集與測(cè)試集，因此就還進(jìn)行了數(shù)據(jù)劃分的步驟。

%% Training Set and Test Set Division

RandomNumber=(randperm(length(Output),floor(length(Output)*0.2)))';
TrainYield=Output;
TestYield=zeros(length(RandomNumber),1);
TrainVARI=Input;
TestVARI=zeros(length(RandomNumber),size(TrainVARI,2));
for i=1:length(RandomNumber)
    m=RandomNumber(i,1);
    TestYield(i,1)=TrainYield(m,1);
    TestVARI(i,:)=TrainVARI(m,:);
    TrainYield(m,1)=0;
    TrainVARI(m,:)=0;
end
TrainYield(all(TrainYield==0,2),:)=[];
TrainVARI(all(TrainVARI==0,2),:)=[];

其中，TrainYield是訓(xùn)練集的因變量，TrainVARI是訓(xùn)練集的自變量；TestYield是測(cè)試集的因變量，TestVARI是測(cè)試集的自變量。

因?yàn)槲疫@里是做估產(chǎn)回歸的，因此變量名稱就帶上了Yield，大家理解即可。

1.4 隨機(jī)森林實(shí)現(xiàn)

這部分代碼其實(shí)比較簡(jiǎn)單。

%% RF

nTree=100;
nLeaf=5;
RFModel=TreeBagger(nTree,TrainVARI,TrainYield,...
    'Method','regression','OOBPredictorImportance','on', 'MinLeafSize',nLeaf);
[RFPredictYield,RFPredictConfidenceInterval]=predict(RFModel,TestVARI);

其中，nTree、nLeaf就是本文1.1部分中我們確定的最優(yōu)樹個(gè)數(shù)與最優(yōu)葉子節(jié)點(diǎn)個(gè)數(shù)，RFModel就是我們所訓(xùn)練的模型，RFPredictYield是預(yù)測(cè)結(jié)果，RFPredictConfidenceInterval是預(yù)測(cè)結(jié)果的置信區(qū)間。

1.5 精度衡量

在這里，我們用RMSE與r衡量模型精度。

%% Accuracy of RF

RFRMSE=sqrt(sum(sum((RFPredictYield-TestYield).^2))/size(TestYield,1));
RFrMatrix=corrcoef(RFPredictYield,TestYield);
RFr=RFrMatrix(1,2);
RFRMSEMatrix=[RFRMSEMatrix,RFRMSE];
RFrAllMatrix=[RFrAllMatrix,RFr];
if RFRMSE<400
    disp(RFRMSE);
    break;
end
disp(RFCycleRun);
str=['Random Forest is Solving...',num2str(100*RFCycleRun/RFRunNumSet),'%'];
waitbar(RFCycleRun/RFRunNumSet,RFScheduleBar,str);
end
close(RFScheduleBar);

在這里，我定義了當(dāng)RMSE滿足<400這個(gè)條件時(shí)，模型將自動(dòng)停止；否則將一直執(zhí)行到本文1.2部分中我們指定的次數(shù)。其中，模型每一次運(yùn)行都會(huì)將RMSE與r結(jié)果記錄到對(duì)應(yīng)的矩陣中。

1.6 變量重要程度排序

接下來，我們結(jié)合RF算法的一個(gè)功能，對(duì)所有的輸入變量進(jìn)行分析，去獲取每一個(gè)自變量對(duì)因變量的解釋程度。

%% Variable Importance Contrast

VariableImportanceX={};
XNum=1;
% for TifFileNum=1:length(TifFileNames)
%     if ~(strcmp(TifFileNames(TifFileNum).name(4:end-4),'MaizeArea') | ...
%             strcmp(TifFileNames(TifFileNum).name(4:end-4),'MaizeYield'))
%         eval(['VariableImportanceX{1,XNum}=''',TifFileNames(TifFileNum).name(4:end-4),''';']);
%         XNum=XNum+1;
%     end
% end

for i=1:size(Input,2)
    eval(['VariableImportanceX{1,XNum}=''',i,''';']);
    XNum=XNum+1;
end

figure('Name','Variable Importance Contrast');
VariableImportanceX=categorical(VariableImportanceX);
bar(VariableImportanceX,RFModel.OOBPermutedPredictorDeltaError)
xtickangle(45);
set(gca, 'XDir','normal')
xlabel('Factor');
ylabel('Importance');

這里代碼就不再具體解釋了，大家會(huì)得到一幅圖，是每一個(gè)自變量對(duì)因變量的重要程度，數(shù)值越大，重要性越大。

其中，我注釋掉的這段是依據(jù)我當(dāng)時(shí)的數(shù)據(jù)情況來的，大家就不用了。

更新：

這里請(qǐng)大家注意，上述代碼中我注釋掉的內(nèi)容，是依據(jù)每一幅圖像的名稱對(duì)重要性排序的X軸（也就是VariableImportanceX）加以注釋（我當(dāng)時(shí)做的是依據(jù)遙感圖像估產(chǎn)，因此每一個(gè)輸入變量的名稱其實(shí)就是對(duì)應(yīng)的圖像的名稱），所以使得得到的變量重要性柱狀圖的X軸會(huì)顯示每一個(gè)變量的名稱。大家用自己的數(shù)據(jù)來跑的時(shí)候，可以自己設(shè)置一個(gè)變量名稱的字段元胞然后放到VariableImportanceX，然后開始figure繪圖；如果在輸入數(shù)據(jù)的特征個(gè)數(shù)（也就是列數(shù)）比較少的時(shí)候，也可以用我上述代碼中間的這個(gè)for i=1:size(Input,2)循環(huán)——這是一個(gè)偷懶的辦法，也就是將重要性排序圖的X軸中每一個(gè)變量的名稱顯示為一個(gè)正方形，如下圖紅色圈內(nèi)。這里比較復(fù)雜，因此如果大家這一部分沒有搞明白或者是一直報(bào)錯(cuò)，在本文下方直接留言就好~

1.7 保存模型

接下來，就可以將合適的模型保存。

%% RF Model Storage

RFModelSavePath='G:\CropYield\02_CodeAndMap\00_SavedModel\';
save(sprintf('%sRF0410.mat',RFModelSavePath),'nLeaf','nTree',...
    'RandomNumber','RFModel','RFPredictConfidenceInterval','RFPredictYield','RFr','RFRMSE',...
    'TestVARI','TestYield','TrainVARI','TrainYield');

其中，RFModelSavePath是保存路徑，save后的內(nèi)容是需要保存的變量名稱。

2 完整代碼

完整代碼如下：

%% Number of Leaves and Trees Optimization
for RFOptimizationNum=1:5
    
RFLeaf=[5,10,20,50,100,200,500];
col='rgbcmyk';
figure('Name','RF Leaves and Trees');
for i=1:length(RFLeaf)
    RFModel=TreeBagger(2000,Input,Output,'Method','R','OOBPrediction','On','MinLeafSize',RFLeaf(i));
    plot(oobError(RFModel),col(i));
    hold on
end
xlabel('Number of Grown Trees');
ylabel('Mean Squared Error') ;
LeafTreelgd=legend({'5' '10' '20' '50' '100' '200' '500'},'Location','NorthEast');
title(LeafTreelgd,'Number of Leaves');
hold off;

disp(RFOptimizationNum);
end

%% Notification
% Set breakpoints here.

%% Cycle Preparation
RFScheduleBar=waitbar(0,'Random Forest is Solving...');
RFRMSEMatrix=[];
RFrAllMatrix=[];
RFRunNumSet=50000;
for RFCycleRun=1:RFRunNumSet

%% Training Set and Test Set Division
RandomNumber=(randperm(length(Output),floor(length(Output)*0.2)))';
TrainYield=Output;
TestYield=zeros(length(RandomNumber),1);
TrainVARI=Input;
TestVARI=zeros(length(RandomNumber),size(TrainVARI,2));
for i=1:length(RandomNumber)
    m=RandomNumber(i,1);
    TestYield(i,1)=TrainYield(m,1);
    TestVARI(i,:)=TrainVARI(m,:);
    TrainYield(m,1)=0;
    TrainVARI(m,:)=0;
end
TrainYield(all(TrainYield==0,2),:)=[];
TrainVARI(all(TrainVARI==0,2),:)=[];

%% RF
nTree=100;
nLeaf=5;
RFModel=TreeBagger(nTree,TrainVARI,TrainYield,...
    'Method','regression','OOBPredictorImportance','on', 'MinLeafSize',nLeaf);
[RFPredictYield,RFPredictConfidenceInterval]=predict(RFModel,TestVARI);
% PredictBC107=cellfun(@str2num,PredictBC107(1:end));

%% Accuracy of RF
RFRMSE=sqrt(sum(sum((RFPredictYield-TestYield).^2))/size(TestYield,1));
RFrMatrix=corrcoef(RFPredictYield,TestYield);
RFr=RFrMatrix(1,2);
RFRMSEMatrix=[RFRMSEMatrix,RFRMSE];
RFrAllMatrix=[RFrAllMatrix,RFr];
if RFRMSE<1000
    disp(RFRMSE);
    break;
end
disp(RFCycleRun);
str=['Random Forest is Solving...',num2str(100*RFCycleRun/RFRunNumSet),'%'];
waitbar(RFCycleRun/RFRunNumSet,RFScheduleBar,str);
end
close(RFScheduleBar);

%% Variable Importance Contrast
VariableImportanceX={};
XNum=1;
% for TifFileNum=1:length(TifFileNames)
%     if ~(strcmp(TifFileNames(TifFileNum).name(4:end-4),'MaizeArea') | ...
%             strcmp(TifFileNames(TifFileNum).name(4:end-4),'MaizeYield'))
%         eval(['VariableImportanceX{1,XNum}=''',TifFileNames(TifFileNum).name(4:end-4),''';']);
%         XNum=XNum+1;
%     end
% end

for i=1:size(Input,2)
    eval(['VariableImportanceX{1,XNum}=''',i,''';']);
    XNum=XNum+1;
end

figure('Name','Variable Importance Contrast');
VariableImportanceX=categorical(VariableImportanceX);
bar(VariableImportanceX,RFModel.OOBPermutedPredictorDeltaError)
xtickangle(45);
set(gca, 'XDir','normal')
xlabel('Factor');
ylabel('Importance');

%% RF Model Storage
RFModelSavePath='G:\CropYield\02_CodeAndMap\00_SavedModel\';
save(sprintf('%sRF0410.mat',RFModelSavePath),'nLeaf','nTree',...
    'RandomNumber','RFModel','RFPredictConfidenceInterval','RFPredictYield','RFr','RFRMSE',...
    'TestVARI','TestYield','TrainVARI','TrainYield');

至此，大功告成。

到此這篇關(guān)于Matlab利用隨機(jī)森林(RF)算法實(shí)現(xiàn)回歸預(yù)測(cè)詳解的文章就介紹到這了,更多相關(guān)Matlab回歸預(yù)測(cè)內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

您可能感興趣的文章: