matlab小白的代码库——多元回归问题

• 拟合用到的数据是老师给的鲍鱼数据<abalone.csv>，想做一个“整体重量”关于“长度”、“直径”、“高度”的多元回归分析。下面是多元回归的基础代码：

clc;clear 
A1=importdata('abalone.csv');

%   A1.data(101:end,:)=[];     %%  弄十几个数据看一眼
X1=A1.data(:,1:3);                 %%  三个自变量：length，diametar，heigth
y=A1.data(:,4);                       %%  一个因变量：whole weigth

X=[ones(length(y),1),X1(:,1).^2,X1(:,2).^2.5,X1(:,3)];       %%点乘构造高次项
%%研究过图像，比较符合"Y=a*X1^2+b*X2^2.5+X3"的形状

[b, bint, r, rint, stats]=regress(y,X);
stats;              %%stats，最关键的度量

rcoplot(r,rint)      %%刻画目前data的拟合残差

• 下面贴原来的残差图

• 下面打算删除残差不合格的点，思考了半天构造了一个删除函数（untitled.m），功能是将因变量y，自变量矩阵X放进去，以残差为“canshu”为阈值，删除不合格点，返回删除后的y，X矩阵。

function [laterX,latery]=untitled(y,X,canshu)

[~, ~, r, ~, ~]=regress(y,X);
[chang,~]=size(X);
   
    for i=1:chang              %%数据的指针删除code
       if (abs(r(i))>canshu)     %%绝对值的设定也是无奈之举，看上面的残差图，不合格点在0点上下都有
             X(i,:)=[];
             y(i,:)=[];                %%这两句删除
             laterX=X(:,:);
             latery=y(:,:);            %%这两句提取结果
             chang=chang-1;
       end
       if(i>=chang)
          break                       %%指针操作防止 i 溢出
       end
    end 
end

• 构造了空矩阵box存放观察y函数参数b的所有变化

number=1;
box=zeros(700,7);

• 这个是主函数里面调用function删除残差的

for i=0.01:0.001:0.7        %%0.01-0.7，由上面残差图设置的，超过0.9就不合理了，步长0.001还是有点误差的
     
    [boxX1,boxy]=untitled(y,X,i);        %%boxX1和boxy存放去噪后的数据
    
    [b1,~, ~, ~,~]=regress(boxy,boxX1);
  
    lost=((b1(1)-b(1))^2+(b1(2)-b(2))^2*(boxX1(2).^2)...
        +(b1(3)-b(3))^2*(boxX1(3).^2.5)+(b1(4)-b(4))^2*boxX1(4)).^2;
        
    %%构造损失函数lost，经验风险最小化，未考虑结构最小化，也就是没考虑过拟合的问题。
    %%因为第一轮box观察了去噪后的stats值，从数理分析的角度都是合适的，mse与R2都在允许范围内
    %%所以第二轮只观察参数b的变化，来使折线图更加逼近真实，具体思辨过程请看我前一篇文章
    
    ZZH=log(10000*lost);              %%log正则化，因为lost有10^8阶的变化
    
    box(number,1)=i;                    %%box放参数b开始规划
    box(number,2:5)=b1;
    box(number,6)=ZZH;
    
    number=number+1;
    
end

• 下面就是为了求lost的min

    figure(50)                          %%看ZZH图像
    plot(1:length(box(:,5)),box(:,6),'r');   

• 贴ZZH图，可以看出有一个最小值

    
    [h,hh]=min(box(:,6))   %%用min函数求ZZH的min，因为ZZH的min也就是lost的min，及想求的lost的min
    %%h，hh存放这个min值时的阈值，求最优阈值下的y，X矩阵，再次拟合运算及为所求的最优参数b11
    
    [boxX11,boxy1]=untitled(y,X,box(hh,1));
    [b11,~, r11, rint11,~]=regress(boxy1,boxX11);
    
    figure(60)
    Z0=b(1)+b(2)*boxX11(:,1).^2+b(3)*boxX11(:,2).^2.5+b(4)*boxX11(:,3); 
     %%未去噪时的“拟合—原数据对比图”
    Z=b11(1)+b11(2)*boxX11(:,1).^2+b11(3)*boxX11(:,2).^2.5+b11(4)*boxX11(:,3);
    %%去噪后，求最优参数后的    “拟合—原数据对比图”
    plot(boxX11,Z,'b',boxX11,Z0,'y') 
    
    figure(70)    
    rcoplot(r11,rint11)
    %%去噪后，求最优参数后的残差图

• 下面贴图对比
  （两个函数我设置了两种颜色的，看不出来是因为拟合的太好而盖住了）

感觉残差不合格点少了一点哈，想完全删除是不可能的，效果图的大致覆盖说明了这一次拟合是非常好的，做到了拟合现在，预测未来。

（希望我的代码能帮到大家，由于我还是个小白，code有错，设计思路，数理分析有问题也很正常，欢迎大家提出来。有错误才有进步）

最后凑不要脸的求一个赞hhhh

- （2020/4/30）更新：
上午老师看了我的博客给我提了几点意见
1.模型的级数是怎么看出来的，2.数据是不能随便删除的

我想了一下确实。首先讲第二点：

data不能随便删除，以为每个数据都有存在的意义。我的拟合优度已经达到0.94了已经非常好了，所以data的额外删除是没必要的，可以容忍。你把部分data删除而留下自己需要的数据，强求模型的优度而这样做，不是鸵鸟行为么？

而如果拟合优度R2过小时，可以考虑删除误差比较大的data来达到比较好的拟合结果。我写的delete函数就有用武之地了。

然后讲第一点：

X1适合平方项，X2我给了2.5次方，X3当成了直线来拟合。

而用直线拟合Y=a*X+b来拟合时优度只有0.81，还是比较合适的。