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

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

1 分解代码

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

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

%% Number of Leaves and Trees Optimizationfor RFOptimizationNum=1:5RFLeaf=[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是为了多次循环,防止最优结果受到随机干扰;大家如果不需要,可以将这句话删除。

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

InputOutput分别是我的输入(自变量)与输出(因变量),大家自己设置即可。

运行后得到下图。

image

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

1.2 循环准备

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

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

其中,RFRMSEMatrixRFrAllMatrix分别用来存放每一次运行的RMSEr结果RFRunNumSet是循环次数,也就是RF运行的次数。

1.3 数据划分

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

%% Training Set and Test Set DivisionRandomNumber=(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 随机森林实现

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

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

其中,nTreenLeaf就是本文1.1部分中我们确定的最优树个数与最优叶子节点个数,RFModel就是我们所训练的模型,RFPredictYield是预测结果,RFPredictConfidenceInterval是预测结果的置信区间。

1.5 精度衡量

在这里,我们用RMSEr衡量模型精度。

%% Accuracy of RFRFRMSE=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<400disp(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部分中我们指定的次数。其中,模型每一次运行都会将RMSEr结果记录到对应的矩阵中。

1.6 变量重要程度排序

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

%% Variable Importance ContrastVariableImportanceX={};
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
% endfor i=1:size(Input,2)eval(['VariableImportanceX{1,XNum}=''',i,''';']);XNum=XNum+1;
endfigure('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轴中每一个变量的名称显示为一个正方形,如下图红色圈内。这里比较复杂,因此如果大家这一部分没有搞明白或者是一直报错,在本文下方直接留言就好~

1.7 保存模型

接下来,就可以将合适的模型保存。

%% RF Model StorageRFModelSavePath='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:5RFLeaf=[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<1000disp(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
% endfor i=1:size(Input,2)eval(['VariableImportanceX{1,XNum}=''',i,''';']);XNum=XNum+1;
endfigure('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');

至此,大功告成。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/448883.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

算法练习-二叉树的节点个数【完全/普通二叉树】(思路+流程图+代码)

难度参考 难度&#xff1a;中等 分类&#xff1a;二叉树 难度与分类由我所参与的培训课程提供&#xff0c;但需要注意的是&#xff0c;难度与分类仅供参考。且所在课程未提供测试平台&#xff0c;故实现代码主要为自行测试的那种&#xff0c;以下内容均为个人笔记&#xff0c;旨…

2024年美赛B题潜水器定位和搜救建模代码和完整论文文档

目前已完成2024年美赛B题潜水器定位和搜救的建模代码和论文编写&#xff0c;部分文章内容和代码如下&#xff1a; 摘要 在海洋探险和搜救领域&#xff0c;潜水器的定位和搜救任务具有重要意义。本文旨在开发一系列模型来预测潜水器位置、分析不确定性、确定信息传递策略、建议…

[UI5 常用控件] 05.FlexBox, VBox,HBox,HorizontalLayout,VerticalLayout

文章目录 前言1. FlexBox布局控件1.1 alignItems 对齐模式1.2 justifyContent 对齐模式1.3 Direction1.4 Sort1.5 Render Type1.6 嵌套使用1.7 组件等高显示 2. HBox,VBox3. HorizontalLayout&#xff0c;VerticalLayout 前言 本章节记录常用控件FlexBox,VBox,HBox,Horizontal…

说说RDB和AOF

简介&#xff1a; 众所周知&#xff0c;redis是一个内存数据库&#xff0c;当机器重启后&#xff0c;内存中数据都会丢失。所以redis提供了两种持久化方式&#xff0c;即&#xff1a;RDB(保存一个时间点前的数据)和AOF(保存redis服务器端执行的每一条命令)。 RDB: RDB有两种…

Web实战丨基于django+hitcount的网页计数器

文章目录 写在前面Django简介主要程序运行结果系列文章写在后面 写在前面 本期内容 基于djangohitcount的网页计数器 所需环境 pythonpycharm或vscodedjango 下载地址 https://download.csdn.net/download/m0_68111267/88795611 Django简介 Django 是一个开源的、基于 …

DES加密原理

DES加密算法综合运用了置换、代替、代数等多种密码技术&#xff0c;具有设计精 巧、实现容易、使用方便等特点。DES加密算法的明文、密文和密钥的分组长度 都是64位&#xff0c;详细的DES加密算法结构如图6-10所示。 图6-10 DES加密算法结构图 DES加密过程如下所示&#xff…

Linux系统漏洞一键检测与修复工具

支持检测及修复漏洞的列表 OpenSSL CVE-2021-3712 OpenSSH CVE-2021-41617 sudo CVE-2021-3156 glibc CVE-2018-11236 polkit CVE-2021-4034 wget CVE-2017-13090 kernel CVE-2016-5195 bash CVE-2016-7543 samba CVE-2021-…

代码随想录第二十四天| 回溯算法● 理论基础 ● 77. 组合

文章目录 理论基础![在这里插入图片描述](https://img-blog.csdnimg.cn/direct/09da30301c104f02baf792ccbf39da15.png)效率回溯法解决的问题如何理解回溯法回溯法模板 77.组合思路&#xff1a;回溯法三部曲 代码&#xff1a;思路-剪枝代码&#xff1a; 理论基础 效率 虽然回…

后端——go系统学习笔记(不断更新中......)

数组 固定大小 初始化 arr1 : [3]int{1, 2, 3} arr2 : [...]int{1, 2, 3} var arr3 []int var arr4 [4]int切片 长度是动态的 初始化 arr[0:3] slice : []int{1,2,3} slice : make([]int, 10)len和cap len是获取切片、数组、字符串的长度——元素的个数cap是获取切片的容量—…

Elasticsearch-内存结构

ElasticSearch的内存从大的结构可以分堆内存&#xff08;On Heap&#xff09;和堆外内存&#xff08;Off Heap&#xff09;。Off Heap部分由Lucene进行管理。On Heap部分存在可GC部分和不可GC部分&#xff0c;可GC部分通过GC回收垃圾对象&#xff0c;从而释放内存。不可GC部分不…

手机云控制发电机组 有网络随时随地操控监控运行

GenCloudTM 发电机组云控系统简介 Ver2.0 目录 公司简介…… …………………………… ………………………………………………1概 述…… …………………………… ………………………………………………1主要功能及特点………… …………… ………… ………………………………

Halcon C++ 环境与配置

Halcon C 环境与配置 1、环境设置相关 头文件路径添加 D:\MVTec\HALCON-22.11-Steady\include D:\MVTec\HALCON-22.11-Steady\include\halconcpp D:\MVTec\HALCON-22.11-Steady\include\halconclab文件添加 D:\MVTec\HALCON-22.11-Steady\lib\x64-win64link添加路径 D:\MV…