Matlab实现随机森林、神经网络、Lasso回归

实现之前所说的上传Matlab实现随机森林、神经网络、Lasso回归的承诺。

Lasso

lasso具有降维的功能，但区别于PCA，lasso直接减少特征数，做的是特征选择，PCA是通过空间转换将特征空间从高维空间转换到低维空间，是降维（PCA的缺点是丢失语意）。当特征有很强的语意的时候，用LASSO更好，后续的分析会更高的保持可解释性，可以给出y=wx表达式（y、w、x均为向量）；反之，数据的语意性不强，PCA更好。

%清空命令窗口
clc;
%清空工作区的变量和值
clear;
%关闭所有程序运行过程中打开的绘图窗口，此处为绘图窗口。
close all;
%读取excel数据，
x1=xlsread("C:\xxx.xlsx",'C4:J328');
y1=xlsread("C:\xxx.xlsx",'K4:K328');
%调用lasso函数，其中参数为（自变量x，因变量y，使用交叉验证，10折交叉验证，α，α=0为岭回归，α=1为lasso回归）；返回值为（b：权重系数，fitinfo:模型信息）
[b,fitinfo] = lasso(x1,y1,'CV',10,'Alpha',1); 
%画图
axTrace = lassoPlot(b,fitinfo,'PlotType','Lambda','XScale','log');
%为图像添加图例，位置在图的右侧外面
legend('show','Location','EastOutside');
axCV = lassoPlot(b,fitinfo,'PlotType','CV');
%寻找最小误差对应的迭代次数
lam1SE = fitinfo.Index1SE;
%最小误差的值
mse_1=fitinfo.MSE(lam1SE);
%取最小误差对应的系数；b矩阵lam1SE列所有行
mat=b(:,lam1SE);
%寻找系数中的非零项（~=0为不等于0）
[row1SE, ] = find(b(:,lam1SE)~=0);
%计算原来的最小均方误差
rhat = x1\y1;
res = x1*rhat - y1;   
MSEmin_real= res'*res/325; 
%最小均方误差对应的迭代次数，上面误差是1se这里是mse
lamMinMSE = fitinfo.IndexMinMSE;
%主成分的系数
matMinMSE = b(:,lamMinMSE);  
%寻找非零自变量的下标（即主成分的下标）
[rowMinMSE, ] = find(b(:,lamMinMSE)~=0);
%两种计算误差方式使得所降成的维度不同，根据自己需求比较两个误差计算方式的差异选择留几个变量。

随机森林RF

随机森林可以计算对于最后的结果贡献率最高的特征，对于回归可以得出线性模型。

clc;
clear;
in=xlsread("C:\xxx.xlsx",'C4:J328');
Out=xlsread("C:\xxx.xlsx",'K4:K328');
%此处未进行归一化，切记，使用归一化的数据进行训练，后面要记得反归一化，且数据多处进行归一化，不可遗漏
nTree = 80;%树的个数
%训练模型
Factor = TreeBagger(nTree, In, Out,'Method','regression','NumPredictorsToSample',122,'OOBpredictorImportance','on','MinLeafSize',5);%
%性能评估，k-fold交叉验证法
[Predict_label,Scores] = predict(Factor, In);
%相关系数
cct=corrcoef(Out,Predict_label');
cct=cct(2,1);

plot(Out,Out,'LineWidth',3);

hold on
scatter(Out,Predict_label,'filled');
hold off
grid on
view(Factor.Trees{1},'Mode','graph')
%变量重要性直方图
weights=Factor.OOBPermutedVarDeltaError;
%重要性降序
[B,iranked] = sort(weights,'descend');
figure
xlabel('变量重要性','FontSize',30,'Interpreter','latex');
ylabel('排序','FontSize',30,'Interpreter','latex');
title('变量相对重要性展示','FontSize',17,'Interpreter','latex');
hold on
%45个数据为黄色
barh(weights(iranked(1:45)),'y');
%15个最重要属性为红色
barh(weights(iranked(1:15)),'r');

grid on 
xt = get(gca,'XTick');    %坐标轴刻度
xt_spacing=unique(diff(xt));
xt_spacing=xt_spacing(1);    
yt = get(gca,'YTick');    
ylim([0 40]); %y轴显示数目上限
xl=xlim;
xlim([0 2.5]);
%为直方图的每个线条添加文字
 for ii=1:40
    text( max([0 weights(iranked(ii))+0.02*max(weights)]),ii,['变量' num2str(iranked(ii))],'Interpreter','latex','FontSize',10);
 end

set(gca,'FontSize',16)
set(gca,'XTick',0:2*xt_spacing:1.1*max(xl));
set(gca,'YTick',yt);
set(gca,'TickDir','out');%坐标轴刻度方向
set(gca, 'ydir', 'reverse' )%最大值在上方
set(gca,'LineWidth',2);   %设置线条粗细
drawnow

%误差变化图
figure
plot(Factor.oobError,'LineWidth',2);
xlabel('生长树数量','FontSize',30)
ylabel('袋外错误率','FontSize',30)
title('袋外错误率','FontSize',30)
set(gca,'FontSize',16)
set(gca,'LineWidth',2);   
grid on
drawnow

神经网络 NN

神经网络将数据进行糅合，最后所得为模型，具有不可解释性，只可实现输入数据并进行预测。
搬运工上线，很多解说该程序的博客，就不多bb啦~

clc;
clear;
close all;
sample=xlsread("C:\Users\CHANG\Desktop\c\sampleR_40.xlsx",'B2:AP326');
p=sample(:,1:end-1);
t=sample(:,end);
p=p';%神经网络输入的行为一个特征的所有样本值
t=t';%目标应为行向量
[p1,ps]=mapminmax(p); 
[t1,ts]=mapminmax(t);
[trainsample.p1,valsample.p1,testsample.p1] =dividerand(p1,0.7,0.15,0.15); 
[trainsample.t1,valsample.t1,testsample.t1] =dividerand(t1,0.7,0.15,0.15);
TF1='tansig';TF2='purelin';
net=newff(p1,t1,layer,{TF1 TF2},'traingdm');%网络创建
%网络参数的设置
net.trainParam.epochs=10000;%训练次数设置   
net.trainParam.goal=1e-7;%训练目标设置    
net.trainParam.lr=0.01;%学习率设置,应设置为较少值，太大虽然会在开始加快收敛速度，但临近最佳点时，会产生动荡，而致使无法收敛    
net.trainParam.mc=0.9;%动量因子的设置，默认为0.9     
net.trainParam.show=25;%显示的间隔次数       
net.trainFcn='trainlm';
[net,tr]=train(net,trainsample.p1,trainsample.t1);
%计算仿真，其一般用sim函数   
[normtrainoutput,trainPerf]=sim(net,trainsample.p1,[],[],trainsample.t1);%训练的数据，根据BP得到的结果    
[normvalidateoutput,validatePerf]=sim(net,valsample.p1,[],[],valsample.t1);%验证的数据，经BP得到的结果    
[normtestoutput,testPerf]=sim(net,testsample.p1,[],[],testsample.t1);     
trainoutput=mapminmax('reverse',normtrainoutput,ts);    
validateoutput=mapminmax('reverse',normvalidateoutput,ts);   
testoutput=mapminmax('reverse',normtestoutput,ts);
%正常输入的数据的反归一化的处理，得到其正式值    
trainvalue=mapminmax('reverse',trainsample.t1,ts);%正常的验证数据
validatevalue=mapminmax('reverse',valsample.t1,ts);%正常的验证的数据     
testvalue=mapminmax('reverse',testsample.t1,ts);%正常的测试数据    
pnew=[159.0     90.2    50.21   0.5680  1019.1993   206.8057    394.2145    102.1239    120.8178    51.6357     69.5808     60.4115     -0.1513     2.0921  -0.3742     49.9826     0.0000  282.6794    921079.3000     131.3584    25.8939     3.8522  24.4931     283.2769    -4.9646     66.4584     139.4476    70.4751     170.5868    0.2601  -5.1580     474.5039    12265800.0000   9.1833  18.7165     56.4848     3.1464  2.2574  31.8697     0.4928 ];   
pnew=pnew';
%将输入数据归一化  
pnewn=mapminmax(pnew);    
anewn=sim(net,pnewn);    
anew=mapminmax('reverse',anewn,ts);    
predictA(layer-layerBegin+1,trainNum)=anew;
%绝对误差的计算   
errorSquaresSum=zeros(layerEnd,trainEnd);
for i=1:(layerBegin-layerEnd+1)    
	for j=trainBegin:trainEnd       
		 errorSquaresSum(i,j)=(predictA(i,j)-88.80)^2;  
	end
end
errors=trainvalue-trainoutput;
%plotregression拟合图
figure,plotregression(trainvalue,trainoutput)
%误差图
figure,plot(1:length(errors),errors,'-b')
title('误差变化图')
%误差值的正态性的检验
figure,hist(errors);%频数直方图 
figure,normplot(errors);
figure,normplot(errors);%Q-Q图
[muhat,sigmahat,muci,sigmaci]=normfit(errors);
%参数估计 均值,方差,均值的0.95置信区间,方差的0.95置信区间
[h1,sig,ci]= ttest(errors,muhat);
%假设检验
figure, ploterrcorr(errors);
%绘制误差的自相关图 
figure, parcorr(errors);
%绘制偏相关图% %输入数据
save p.mat; 
save t.mat;
       

补充

当系数矩阵和y的维数不呼应时，求x*a=y最小二乘解，使用命令为：左除
代码：

x=xlsread("C:\xxx.xlsx",'B2:T326');
y=xlsread("C:\xxx.xlsx",'U2:U326');
//组成增广矩阵
plus=[x y];
//求x的秩
rX=rank(x);
//求增广矩阵的秩
rP=rank(plus);
k_Left=x\y;
y_new_Left=x*k_Left;
//计算误差。误差平方和
error_Left=sum((y-y_new_Left).^2);