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

news2024/12/24 18:00:10

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

1 分解代码

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

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

%% 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是为了多次循环,防止最优结果受到随机干扰;大家如果不需要,可以将这句话删除。

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

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

运行后得到下图。

image

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

1.2 循环准备

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

%% Cycle Preparation

RFScheduleBar=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 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);

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

1.5 精度衡量

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

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

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轴中每一个变量的名称显示为一个正方形,如下图红色圈内。这里比较复杂,因此如果大家这一部分没有搞明白或者是一直报错,在本文下方直接留言就好~

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

至此,大功告成。

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

Redis常用数据结构与应用场景

常用数据结构 StringHashListSetZset String常用操作 String应用场景 Hash常用操作 hash应用场景 Hash结构优缺点 优点 同类数据归类整合存储,方便数据管理相比String操作消耗内存与spu更小相比string更节省空间 缺点 过期功能不能使用在field上,只用用在key上Redis集群…

Navicate 连接云服务器MySQL

Navicate 连接云服务器MySQL 1.打开Navicate,点击左上角的连接,选择MySQL 第一步:第一个页面是常规,按照图上的标注填写 第二步,点击 SSH ,进入下面的页面 第三步&#xff0c;点击测试连接

FreeRTOS任务相关的API函数

本篇文章记录我学习FreeRTOS的任务相关的API函数。主要涉及FreeRTOS的任务创建和删除函数、任务挂起和恢复函数。希望我的分享对你有所帮助。 读者如果需要实战FreeRTOS动态/静态任务创建和删除&#xff0c;可以参考以下文章&#xff1a; FreeRTOS动态 / 静态创建和删除任务-CS…

C++初阶 内存管理和模板

目录 一、new 1.1什么是new&#xff1f; 1.2为什么要有new&#xff1f; 1.3使用new 1.4 new的超级好处 二、delete 2.1什么是delete&#xff1f; 2.2为什么要有delete&#xff1f; 2.3使用delete 三、 malloc / free和new / delete的共同点和区别 四、浅谈模板 4.1什…

嵌入式学习第三篇——51单片机

目录 1&#xff0c;嵌入式系统 1&#xff0c;嵌入式系统的定义 2&#xff0c;单片机的定义 2&#xff0c;51单片机 1&#xff0c;开发环境 2&#xff0c;开发板使用的基本思路 1&#xff0c;查看原理图&#xff0c;查看芯片手册 2&#xff0c;获得调用硬件的管…

使用Arcgis对欧洲雷达高分辨率降水数据重投影

当前需要使用欧洲高分辨雷达降水数据&#xff0c;但是这个数据的投影问题非常头疼。实际的投影应该长这样&#xff08;https://gist.github.com/kmuehlbauer/645e42a53b30752230c08c20a9c964f9?permalink_comment_id2954366https://gist.github.com/kmuehlbauer/645e42a53b307…

4K Video Downloader forMac/win:畅享高清视频下载的终极利器!

在如今的数字时代&#xff0c;高清视频已经成为人们生活中不可或缺的一部分。无论是观看精彩的电影、音乐视频&#xff0c;还是学习教育类的在线课程&#xff0c;我们都希望能够以最清晰流畅的方式来欣赏。而为了满足这一需求&#xff0c;我们需要一款功能强大的高清视频下载软…

一文读懂C++的类和对象以及多态的原理

现实生活中&#xff0c;关于类和对象最好的例子是自然界的动物类&#xff0c;本文将以此为场景逐步引入C的概念&#xff0c;达到学习的目的。因为C这门语言本身有很多繁杂的内容&#xff0c;而网上的资源也是参差不齐&#xff0c;有的人见山谈山遇水聊水&#xff0c;有多人故弄…

Vue3 - 从 vue2 到 vue3 过渡,这一套就够了(案例 + 效果演示)(二)

目录 一、组合式 API 的使用 1.1、watch 函数 1.2、watchEffect 函数 1.3、toRef 和 toRefs 1.3.1、toRef 1.3.2、toRefs 1.4、vue3 的声明周期 一、组合式 API 的使用 1.1、watch 函数 与 vue2.x 中的 watch 配置功能一致&#xff0c;但是多了一些坑&#xff1a; 这…

大数据环境搭建(一)-Hive

1 hive介绍 由Facebook开源的,用于解决海量结构化日志的数据统计的项目 本质上是将HQL转化为MapReduce、Tez、Spark等程序 Hive表的数据是HDFS上的目录和文件 Hive元数据 metastore&#xff0c;包含Hive表的数据库、表名、列、分区、表类型、表所在目录等。 根据Hive部署模…

10 分钟在K8s 中部署轻量级日志系统 Loki

转载至我的博客 https://www.infrastack.cn &#xff0c;公众号&#xff1a;架构成长指南 Loki 是什么&#xff1f; Loki是由Grafana Labs开源的一个水平可扩展、高可用性&#xff0c;多租户的日志聚合系统的日志聚合系统。它的设计初衷是为了解决在大规模分布式系统中&#x…

WebGL+Three.js入门与实战——绘制水平移动的点、通过鼠标控制绘制(点击绘制、移动绘制、模拟画笔)

个人简介 &#x1f440;个人主页&#xff1a; 前端杂货铺 &#x1f64b;‍♂️学习方向&#xff1a; 主攻前端方向&#xff0c;正逐渐往全干发展 &#x1f4c3;个人状态&#xff1a; 研发工程师&#xff0c;现效力于中国工业软件事业 &#x1f680;人生格言&#xff1a; 积跬步…

MATLAB矩阵的操作(第二部分)

师从清风 矩阵的创建方法 在MATLAB中&#xff0c;矩阵的创建方法主要有三种&#xff0c;分别是&#xff1a;直接输入法、函数创建法和导入本地文件中的数据。 直接输入法 输入矩阵时要以中括号“[ ]”作为标识符号&#xff0c;矩阵的所有元素必须都在中括号内。 矩阵的同行元…

IEC104 S帧超时判定客户与服务端不匹配造成的异常链接问题分析

2、通过ss命令发现确有链接端口变化&#xff0c;与设备约一天一次的重连&#xff0c;通过抓包&#xff08;tcpdump -vvv -nn port 1001 -w 0926.cap&#xff09;分析得以下现象 2.1、异常情况时未对设备的I帧均匀的回S帧进行确认&#xff0c;正常情况时均匀的回S帧进行确认 2.…

数据在内存中的存储(上)

1. 整数在内存中的存储 整数的2进制表示方法有三种&#xff1a;即原码、反码和补码 三种表示方法均有符号位和数值位两部分&#xff0c;符号位都是用0表示“正”&#xff0c;用1表示“负”&#xff0c;而数值位最 高位的一位是被当做符号位&#xff0c;剩余的都是数值位。 正…

深度学习水论文怎么缝模块?搭积木永不过时!(附80个即插即用模型)

深度学习如何创新&#xff1f;如何水模型&#xff1f;总结来说就八个字&#xff1a;排列组合&#xff0c;会讲故事。说直白点&#xff0c;就是缝模块。 先看看别人怎么做&#xff0c;然后根据自己的实际情况将这些模块来一波随机组合&#xff0c;这样效率会高很多。我这边已经…

冰冻天气恰逢春运,“观冰精灵”化身电力供应守护者

据中国路网&#xff0c;截至2月1日14时&#xff0c;受降雪及路面结冰影响&#xff0c;河北、山西、内蒙古、黑龙江、江苏、安徽、河南、山东、西藏、陕西、宁夏、甘肃、新疆共封闭路段66个&#xff08;涉及44条高速公路、5条普通国道、5条普通省道&#xff09;&#xff0c;关闭…

Pandas.Series.clip() 修剪数值范围 详解 含代码 含测试数据集 随Pandas版本持续更新

关于Pandas版本&#xff1a; 本文基于 pandas2.2.0 编写。 关于本文内容更新&#xff1a; 随着pandas的stable版本更迭&#xff0c;本文持续更新&#xff0c;不断完善补充。 传送门&#xff1a; Pandas API参考目录 传送门&#xff1a; Pandas 版本更新及新特性 传送门&…

前端JavaScript篇之map和Object的区别、map和weakMap的区别

目录 map和Object的区别map和weakMap的区别 map和Object的区别 Object是JavaScript的内置对象&#xff0c;用于存储键值对。Object的键必须是字符串或符号&#xff0c;值可以是任意类型。Map是ES6引入的新数据结构&#xff0c;用于存储键值对。Map的键可以是任意类型&#xff…

高中数学两面角习题练习1

用到的定理 2 第1问证明&#xff1a; 第2问用到的知识和例子&#xff1a; 二面角锐角钝角判断的快速技巧 https://www.bilibili.com/video/BV13P41157K1/?spm_id_from333.788.recommend_more_video.0&vd_source91b03ee59c462b7b3cfbd57346cf1001 叉乘的几何意义及应用 …