MATLAB实现随机森林（RF）回归与自变量影响程度分析

妖狐艹你老母 2024-03-16 12:28 54阅读 0赞

本文介绍基于**MATLAB**，利用**随机森林**（**RF**）算法实现回归预测，以及自变量重要性排序的操作。

*  [1 分解代码][1]
    
     *  [1.1 最优叶子节点数与树数确定][1.1]
     *  [1.2 循环准备][1.2]
     *  [1.3 数据划分][1.3]
     *  [1.4 随机森林实现][1.4]
     *  [1.5 精度衡量][1.5]
     *  [1.6 变量重要程度排序][1.6]
     *  [1.7 保存模型][1.7]
 *  [2 完整代码][2]

本文分为两部分，首先是对代码进行分段、详细讲解，方便大家理解；随后是完整代码，方便大家自行尝试。另外，关于基于**MATLAB**的神经网络（**ANN**）代码与详细解释，我们将在后期博客中介绍。

## 1 分解代码 ##

### 1.1 最优叶子节点数与树数确定 ###

首先，我们需要对**RF**对应的叶子节点数与树的数量加以择优选取。

其中，`RFOptimizationNum`是为了多次循环，防止最优结果受到随机干扰；大家如果不需要，可以将这句话删除。

`RFLeaf`定义初始的叶子节点个数，我这里设置了从`5`到`500`，也就是从`5`到`500`这个范围内找到最优叶子节点个数。

`Input`与`Output`分别是我的输入（自变量）与输出（因变量），大家自己设置即可。

运行后得到下图。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3poZWJ1c2hpYmlhb3NoaWZ1_size_16_color_FFFFFF_t_70][]

首先，我们看到`MSE`最低的线是红色的，也就是`5`左右的叶子节点数比较合适；再看各个线段大概到`100`左右就不再下降，那么树的个数就是`100`比较合适。

### 1.2 循环准备 ###

由于机器学习往往需要多次执行，我们就在此先定义循环。

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

其中，`RFRMSEMatrix`与`RFrAllMatrix`分别用来存放每一次运行的**RMSE**、**r结果**，`RFRunNumSet`是循环次数，也就是**RF**运行的次数。

### 1.3 数据划分 ###

接下来，我们需要将数据划分为训练集与测试集。这里要注意：**RF**其实一般并不需要划分训练集与测试集，因为其可以采用袋外误差（Out of Bag Error，**OOB Error**）来衡量自身的性能。但是因为我是做了多种机器学习方法的对比，需要固定训练集与测试集，因此就还进行了数据划分的步骤。

%% 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`是训练集的因变量，`TrainVARI`是训练集的自变量；`TestYield`是测试集的因变量，`TestVARI`是测试集的自变量。

因为我这里是做估产回归的，因此变量名称就带上了`Yield`，大家理解即可。

### 1.4 随机森林实现 ###

这部分代码其实比较简单。

%% 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**部分中我们确定的最优树个数与最优叶子节点个数，`RFModel`就是我们所训练的模型，`RFPredictYield`是预测结果，`RFPredictConfidenceInterval`是预测结果的置信区间。

### 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);

在这里，我定义了当**RMSE**满足`<400`这个条件时，模型将自动停止；否则将一直执行到本文**1.2**部分中我们指定的次数。其中，模型每一次运行都会将**RMSE**与**r结果**记录到对应的矩阵中。

### 1.6 变量重要程度排序 ###

接下来，我们结合**RF**算法的一个功能，对所有的输入变量进行分析，去获取每一个自变量对因变量的解释程度。

%% 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');

这里代码就不再具体解释了，大家会得到一幅图，是每一个自变量对因变量的重要程度，数值越大，重要性越大。

其中，我注释掉的这段是依据我当时的数据情况来的，大家就不用了。

**更新**：这里请大家注意，上述代码中我注释掉的内容，是依据每一幅图像的名称对重要性排序的`X`轴（也就是`VariableImportanceX`）加以注释（我当时做的是依据遥感图像估产，因此每一个输入变量的名称其实就是对应的图像的名称），所以使得得到的变量重要性柱状图的`X`轴会显示每一个变量的名称。大家用自己的数据来跑的时候，可以自己设置一个变量名称的字段元胞然后放到`VariableImportanceX`，然后开始`figure`绘图；如果在输入数据的特征个数（也就是列数）比较少的时候，也可以用我上述代码中间的这个`for i=1:size(Input,2)`循环——这是一个偷懒的办法，也就是将重要性排序图的X轴中每一个变量的名称显示为一个正方形，如下图红色圈内。这里比较复杂，因此如果大家这一部分没有搞明白或者是一直报错，在本文下方直接留言就好~

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3poZWJ1c2hpYmlhb3NoaWZ1_size_16_color_FFFFFF_t_70 1][]

### 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`后的内容是需要保存的变量名称。

## 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');

至此，大功告成。

[1]: https://www.cnblogs.com/fkxxgis/p/17095368.html#1-%E5%88%86%E8%A7%A3%E4%BB%A3%E7%A0%81
[1.1]: https://www.cnblogs.com/fkxxgis/p/17095368.html#11-%E6%9C%80%E4%BC%98%E5%8F%B6%E5%AD%90%E8%8A%82%E7%82%B9%E6%95%B0%E4%B8%8E%E6%A0%91%E6%95%B0%E7%A1%AE%E5%AE%9A
[1.2]: https://www.cnblogs.com/fkxxgis/p/17095368.html#12-%E5%BE%AA%E7%8E%AF%E5%87%86%E5%A4%87
[1.3]: https://www.cnblogs.com/fkxxgis/p/17095368.html#13-%E6%95%B0%E6%8D%AE%E5%88%92%E5%88%86
[1.4]: https://www.cnblogs.com/fkxxgis/p/17095368.html#14-%E9%9A%8F%E6%9C%BA%E6%A3%AE%E6%9E%97%E5%AE%9E%E7%8E%B0
[1.5]: https://www.cnblogs.com/fkxxgis/p/17095368.html#15-%E7%B2%BE%E5%BA%A6%E8%A1%A1%E9%87%8F
[1.6]: https://www.cnblogs.com/fkxxgis/p/17095368.html#16-%E5%8F%98%E9%87%8F%E9%87%8D%E8%A6%81%E7%A8%8B%E5%BA%A6%E6%8E%92%E5%BA%8F
[1.7]: https://www.cnblogs.com/fkxxgis/p/17095368.html#17-%E4%BF%9D%E5%AD%98%E6%A8%A1%E5%9E%8B
[2]: https://www.cnblogs.com/fkxxgis/p/17095368.html#2-%E5%AE%8C%E6%95%B4%E4%BB%A3%E7%A0%81
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3poZWJ1c2hpYmlhb3NoaWZ1_size_16_color_FFFFFF_t_70]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/16/b5938c29d7ce4f15b8b6aef4b83c401b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3poZWJ1c2hpYmlhb3NoaWZ1_size_16_color_FFFFFF_t_70 1]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/16/a3534af27adc44768c277e722ebb8eaa.png