多元线性回归

leadersheng

基于矩阵运算的多元线性回归分析
x=[x1;x2;x3];                                %向量合并为矩阵
X=[ones(length(y),1) x'];                    %在矩阵中加入常数向量并转置
Y=y';                                        %因变量向量转置
[m,n]=size(x);                               %计算自变量矩阵的行列数
B=inv(X'*X)*X'*Y;                            %计算回归系数
Yp=X*B;                                      %建立预测模型

%计算用于检验的统计量
R2=(abs(B'*X'*Y)-n*mean(y)^2)/(abs(Y'*Y)-n*mean(y)^2);
                                             %计算复相关系数
Radj2=R2-(1-R2)*(m+1)/(n-m-1);               %计算校正相关系数平方
s=sqrt((Y'*Y-B'*X'*Y)/(n-m-1));              %计算标准误差
v=s/mean(y);                                 %计算变异系数
F=(abs(B'*X'*Y)-n*mean(y)^2)/(m*s^2);        %计算F统计量
e=Y-Yp;                                      %计算残差
i=1:n-1;                                     %残差编号
DW=sumsqr(e(i+1)-e(i))/sumsqr(e);            %计算Durbin-Watson统计量

%计算偏自相关系数
i=1:n;j=1:m+1;                               %定义矩阵元素编号
Xy=[x' Y];                                   %将变量合并为一个新矩阵
M=mean(Xy(:,j));                             %计算各个变量的均值
S=std(Xy(:,j));                              %计算各个变量的标准差
Mv=M(ones(n,1),:);                           %均值向量平移为矩阵
Sv=S(ones(n,1),:);                           %标准差向量平移为矩阵
Xs=(Xy-Mv)./Sv;                              %数据标准化
Rs=cov(Xs);                                  %计算简单相关系数矩阵
C=inv(Rs);                                   %计算简单相关系数矩阵的逆矩阵
Cjy=C(:,m+1);                                %提取逆矩阵的末列
Cjj=diag(C);                                 %提取逆矩阵的对角线元素
Pr=-Cjy./((Cjj*C(m+1,m+1)).^0.5);            %计算偏相关系数
Rjy=Pr(1:m);                                 %提取自变量的偏相关系数

%计算t统计量
i=1:n;j=1:m;                                 %定义矩阵元素编号
xt=x';                                       %自变量矩阵转置
M=mean(xt(:,j));                             %计算自变量均值
N=M(ones(n,1),:);                            %均值向量平移为矩阵
Z=xt-N;                                      %变量中心化
P=Z'*Z;                                      %计算自变量的交叉乘积和
D=inv(P);                                    %交叉乘积和矩阵求逆
d=diag(D);                                   %提取逆矩阵的对角线元素
sb=d.^0.5*s;                                 %计算参数标准误差
b=B(2:m+1);                                  %提取回归系数
T=b./sb;                                     %计算t统计量

%给出部分计算结果
B,R2,s,F,DW,T,Rjy                            %给出参数和统计量的计算值

%基于回归计算程序包的多元线性回归分析
x1=[57.82 58.05 59.15 63.83 65.36 67.26 66.92 67.79 75.65 80.57 79.02 80.52 86.88 95.48 109.71 126.5 138.89 160.56];
x2=[27.05 28.89 33.02 35.23    24.94 32.95 30.35 38.70 47.99 54.18 58.73 59.85 64.57 70.97 81.54 94.01 103.23 119.33];
x3=[14.54 16.83 12.26 12.87    11.65 12.87 10.80 10.93 14.71 17.56 20.32 18.67 25.34 25.06 29.69 43.86 48.90 60.98];
y=[3.09 3.40 3.88 3.90 3.22 3.76 3.59 4.03 4.34 4.65 4.78 5.04 5.59 6.01 7.03 10.03 10.83 12.9];
x=[x1;x2;x3];                                %向量合并为矩阵
X=[ones(length(y),1) x'];                    %在矩阵中加入常数向量并转置
Y=y';                                        %因变量向量转置
[m,n]=size(x);                               %计算自变量矩阵的行列数
[B Bint E Eint Stats]=regress(Y,X);          %调用回归分析程序包
s=Stats(4)^0.5;                              %给出标准误差
DW=sumsqr(E(2:n)-E(1:n-1))/sumsqr(E);        %因变量向量转置
rcoplot(E,Eint);                             %绘制残差序列的置信区间图
%计算t统计量
SB=Stats(4)*inv(X'*X);                       %计算一个特殊的方阵
D=diag(SB);                                  %提取方阵对角线元素
sj=D.^0.5;                                   %计算参数标准误差
T=B./sj;                                     %计算t统计量

%计算偏自相关系数
i=1:n;j=1:m+1;                               %定义矩阵元素编号
Xy=[x' Y];                                   %将变量合并为一个新矩阵
M=mean(Xy(:,j));                             %计算各个变量的均值
S=std(Xy(:,j));                              %计算各个变量的标准差
Mv=M(ones(n,1),:);                           %均值向量平移为矩阵
Sv=S(ones(n,1),:);                           %标准差向量平移为矩阵
Xs=(Xy-Mv)./Sv;                              %数据标准化
Rs=cov(Xs);                                  %计算简单相关系数矩阵
C=inv(Rs);                                   %计算简单相关系数矩阵的逆矩阵
Cjy=C(:,m+1);                                %提取逆矩阵的末列
Cjj=diag(C);                                 %提取逆矩阵的对角线元素
Pr=-Cjy./((Cjj*C(m+1,m+1)).^0.5);            %计算偏相关系数
Rjy=Pr(1:m);                                 %提取自变量的偏相关系数

%计算t统计量（2）
MSE=E'*E/(n-m-1);
SB=MSE*inv(X'*X);
D=diag(SB);
T=B./sqrt(D)

%计算共线性容忍度和相应VIF值的程序
x1=[57.82 58.05 59.15 63.83 65.36 67.26 66.92 67.79 75.65 80.57 79.02 80.52 86.88 95.48 109.71 126.5 138.89 160.56];
x2=[27.05 28.89 33.02 35.23    24.94 32.95 30.35 38.70 47.99 54.18 58.73 59.85 64.57 70.97 81.54 94.01 103.23 119.33];
x3=[14.54 16.83 12.26 12.87    11.65 12.87 10.80 10.93 14.71 17.56 20.32 18.67 25.34 25.06 29.69 43.86 48.90 60.98];
x=[x1;x2;x3]';                                      %向量合并为矩阵
[m,n]=size(x');                                     %计算矩阵的行列数
for j=1:m                                           %变量编号
    y=x(:,j);                                       %定义因变量
    for k=1:m                                       %变量再次编号
        x1=x(:,k);                                  %提取一个变量
        for l=1:m                                   %变量三次编号
            x2=x(:,l);                              %再次提取一个变量
            if x2~=x1&x2~=y&x1~=y                   %变量选取条件
                X=[ones(n,1) x1 x2];                %定义自变量矩阵
                [B,Bint,E,Eint,Stats]=regress(y,X); %线性回归
                Tol=1-Stats(1);                     %计算容忍度
                VIF=1/Tol;                          %计算VIF值
            end
        end
    end
   Col=[j Tol VIF]                                  %按照变量编号给出计算结果
end
x=[x1;x2;x3];                                %自变量向量合并为矩阵
[m,n]=size(x);                               %计算矩阵的行列数
i=1:n;j=1:m;                                 %定义矩阵元素编号
X=x';                                        %自变量矩阵转置
M=mean(X(:,j));                              %计算各个变量的均值
S=std(X(:,j));                               %计算各个变量的标准差
Mv=M(ones(n,1),:);                           %均值向量平移为矩阵
Sv=S(ones(n,1),:);                           %标准差向量平移为矩阵
Xs=(X-Mv)./Sv;                               %数据标准化
Rs=cov(Xs);                                  %计算简单相关系数矩阵
C=inv(Rs);                                   %计算简单相关系数矩阵的逆矩阵
VIF=diag(C);                                 %提取对角线上的VIF值
Tol=ones(m,1)./VIF;                          %计算容忍度Tol值
Col=[[j]' Tol VIF];                          %提取自变量的偏相关系数

cassienju

lz厚道，谢谢！

shaha

厉害