MATLAB的code练习（多元回归）

（2020.4.27    更新中：有空慢慢写）

（已经做完了，想看结果的请移步去另一个帖子，这个贴将写一写整个学习过程的思考细节故会比较琐碎。而结果贴只会写最优化的结果）

matlab小白的代码库——多元回归问题_PhD#Chen的博客-CSDN博客

Day.1    今天想学习一下，写一个多元回归问题交作业，先在CSDN找了一段源代码。

clc,clear 

x1=[7 1 11 11 7 11 3 1 2 21 1 11 10];
x2=[26 29 56 31 52 55 71 31 54 47 40 66 68];
x3=[6 15 8 8 6 9 17 22 18 4 23 9 8];
x4=[60 52 20 47 33 22 6 44 22 26 34 12 12];
y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4];
 
X= [ones(length(y),1),x1',x2',x3',x4'];Y=y';
%把行向里转秩为列向里Y=y' ;%把行向量转秩为列向量！！！！最后应该是一个长条形

[b, bint, r, rint, stats]=regress(Y,X) ;

figure(1)
rcoplot(r,rint)     %残差

figure(2)
z=b(1)+b(2)*x1+b(3)*x2+b(4)*x3+b(5)*x4;
plot(X,Y,'k',X,z,'r')   %折线图；观察z曲线的拟合程度

Day.2-3   依葫芦画瓢  思考模型原理，优化了一下

%%自行仿写.02
clc;clear 
A1=importdata('abalone.csv');

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

figure(1)
subplot(221);plot(X1(:,1),y,'k+')
subplot(222);plot(X1(:,2),y,'k+')
subplot(223);plot(X1(:,3),y,'k+')   
%为了观察三个自变量适合用什么次数的函数

figure(2)
plot(X1(:,1),y,'k+',X1(:,2),y,'b+',X1(:,3),y,'g+') 
%k=black  b=blue g=green
%为了观察三个自变量的阶数

X=[ones(length(y),1),X1(:,1).^2,X1(:,2).^2.5,X1(:,3)];
%借鉴了点乘的思路来构造高次项

[b, bint, r, rint, stats]=regress(y,X);

b;      % B：回归系数，是个向量。
bint;   % BINT：回归系数的95%区间估计。
r;    % R：残差。
rint;  % RINT：置信区间。
stats  % STATS：用于检验回归模型的统计量。
% 有4个数值：判定系数R^2，F统计量观测值，检验的p的值!!!!，误差方差的估计。

figure(3);
rcoplot(r,rint)

figure(4)
z=b(1)+b(2)*X1(:,1).^2+b(3)*X1(:,2).^2.5+b(4)*X1(:,3);
plot(X,z,'r',X,y,'k')  

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%下面第一次删除残差不合格的点

    for circulate=2

    [chang,kuan]=size(X1)

    for i=1:chang              %%数据的指针删除code
       if (abs(r(i))>canshu)     %%域值的大小合适，注意stats，b的改变防止过拟合
             r(i)
             X1(i,:)=[];
             y(i,:)=[];
             chang=chang-1;
             [chang,kuan]=size(X1)    %%观察验证比大小
             i                        %%观察验证比大小     
       end

       if(i>=chang)
          break                       %%指针操作防止i溢出
       end

    end    

    X2=[ones(length(y),1),X1(:,1).^2,X1(:,2).^2.5,X1(:,3).^4];
    [b1,bint1, r1, rint1, stats1]=regress(y,X2);
    r1;

    figure(5);
    rcoplot(r1,rint1)   %%不合格的点少了一点，依然很多？？？？？？

    figure(6)
    z1=b1(1)+b1(2)*X1(:,1).^2+b1(3)*X1(:,2).^2.5+b1(4)*X1(:,3);
    plot(X2,z1,'r',X2,y,'k') 
 

    end




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%看不懂的数学分析:
    for yingcang=1
    [n,k]=size(X2)
    A=X2'*X2;             %求算信息阵A，
    C=inv(A);           %求算信息阵的逆阵，
    b1=X2\y              % 求算回归统计数向量，其中第一行为回归截距a,
    SSE=y'*y-b1'*X2'*y   %求算离回归平方和，
    MSe=SSE/(n-k-1)     %求算离回归方差，
    MSr=b1.*b1./diag(C)  %求算偏回归平方和，其中第一行是a与0差异的偏平方和，
    F=MSr/MSe            %F测验，其中第一行为a与0差异的F值，
    sb=sqrt(MSe*diag(C));     %求算回归统计数标准误，
    t=b1./sb            % 回归统计数的 t 测验，其中第一行为a与0差异的t测验值。
    yz=[t, t.^2, F]        %验证t^2=F
    SST=var(y)*(n-1)
    R2=(SST-SSE)/SST     %SST=SSR+SSE   
    SSR=SST-SSE


    %   注：r^2(决定系数)越大（接近于1），它们之间的关系越密切，拟合的效果越好。
    %   当然，对于复杂的多变数非线性关系的分析，统计上应以离回归方差（MSe）最小为佳。
    %   MSe=RSS/(n-k-1), RSS为离回归平方和，n为观察值组数，k为模型的效应项数（不包括常数）
    

Day.4  重新构造，去噪与调参

由于存在残差不合格的点，构造了一个删除函数，以“参数”为阈值来删除残差不合格的点，将剩下的数据继续拟合得到新的b，stats

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

[b, bint, r, rint, stats]=regress(y,X);
stats;

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

for i=0.01:0.001:0.7
     
    [boxX1,boxy]=untitled(y,X,i);
    
    [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;
    
    ZZH=log(10000*lost);
    
    box(number,1)=i;
    box(number,2:5)=b1;
    box(number,6)=ZZH;
    
    number=number+1;
    
end

    figure(50)
    plot(1:length(box(:,5)),box(:,6),'r');   
    
    
    [h,hh]=min(box(:,6))
    
    [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)

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

[~, ~, r, ~, ~]=regress(y,X);
[chang,~]=size(X);
   
    for i=1:chang              %%数据的指针删除code
       if (abs(r(i))>canshu)     %%域值的大小合适，注意stats，b的改变防止过拟合
             X(i,:)=[];
             y(i,:)=[];
             laterX=X(:,:);
             latery=y(:,:);
             chang=chang-1;
       end
       if(i>=chang)
          break                       %%指针操作防止i溢出
       end
    end 
end

仿写的各种问题：

1. csv选题与导入：csv导入，A是结构体数据，如何把data取出来？

2. 数据清洗，A结构体有text与data区分   只选取了四组data：自变量：length，diameter，height  ；因变量：shucked weight
3. 三个自变量可以用一个大X1表示，一拖三
4. y表示自变量，现在回想其实day2的数据清洗意思不大，就是给大X1方便了而已
5. 大X   error了半天，矩阵大小不合格，又翻书，观察源代码，发现大X是一个长方形矩阵，所以矩阵格式如何设定，源代码y=y'的意义要搞清楚。然而任然没彻底看懂数学公式。。。。。。。。。
6. 突发奇想用了几个size和罗列来辅助观察大小
7. b，bint，states 之前有nan过，原因在于XY格式大小的问题，nan肯定是错的，但是第一次没有报错
8. Z的形式改了好几下，两个矩阵的乘法怎么一一对应

Day.3的优化

1. 要有逻辑，先看自变量适合什么形式的函数，在构造多次项回归（一次和二次）
2. 比阶来显示次数
3. 点乘构造高次项！！！！！！！！！
4. 折线图看不懂啥意思
5. 重要收获！！！！！！！

       Data的条件删除，指针操作和size，i两项的观察，以及最后的break

1. 两次去噪效果差不多，删除了第三轮去噪
2. 残差不合格去噪后第二轮发现还有漏网之鱼，也有新的不合格项？？？？？？？？？？直接导致了优化失败
3. 去噪阈值的选择，手工选的，经常过拟合，需要观察stats，b的改变来手工调参
4. 取了一次成功的降噪之后发现函数图像改变不大？？？？？？
5. 建模的需求已经可以满足了，剩下的研究生数理统计就没看

这几天把思考的细节补上，最后凑不要脸的求各位看官给个免费的赞鼓励一下哈哈哈哈