现在的位置: 首页 > 综合 > 正文

iml回归

2018年10月22日 ⁄ 综合 ⁄ 共 3281字 ⁄ 字号 评论关闭

FROM 《SAS/IML USER'S GUIDE》 

/* 已知X,Y,利用最小二乘法估计回归方程Y=XB+E */
*ods trace output;
proc iml;
start reg;
   n = nrow(x);  /* 观测值数目 */
   k = ncol(x);    /* 变量数目 */
   xpx = x`*x;         /* 内积 */
   xpy = x`*y;
   xpxi = inv(xpx); /* 内积矩阵的逆矩阵 */
   b = xpxi*xpy;    /* 参数b的估计 */
   yhat = x*b;      /* y的预测值 */
   resid = y-yhat;  /* 真实值预测值残差 */
   sse = resid`*resid; /* 残差平方和*/
   dfe = n-k; /* 残差的自由度 */
   mse = sse/dfe; /* 均方误差 */
   rmse = sqrt(mse); /*均方标准差 */
   covb = xpxi#mse; /* 估计值b的协方差 */
   stdb = sqrt(vecdiag(covb)); /* 标准误差 */
   t = b/stdb; /* 估计值的t检验,零假设b=0 */
   probt = 1-probf(t#t,1,dfe); /* t检验有显著性的p值 */
   print name b stdb t probt;
   s = diag(1/stdb);
   corrb = s*covb*s; /* 估计值的相关系数矩阵 */
   print ,"Covariance of Estimates", covb[r=name c=name] ,
          "Correlation of Estimates",corrb[r=name c=name] ;
   if nrow(tval)=0 then return; /* 设定t值 */
      projx = x*xpxi*x`; /* hat matrix */
      vresid = (i(n)-projx)*mse; /* covariance of residuals */
      vpred = projx#mse; /* covariance of predicted values */
      h = vecdiag(projx); /* hat leverage values */
      lowerm = yhat-tval#sqrt(h*mse); /* low. conf lim for mean */
      upperm = yhat+tval#sqrt(h*mse); /* upper lim. for mean */
      lower = yhat-tval#sqrt(h*mse+mse); /* lower lim. for indiv*/

      upper = yhat+tval#sqrt(h*mse+mse);/* upper lim. for indiv */
      print ,,"Predicted Values, Residuals, and Limits" ,,
              y yhat resid h lowerm upperm lower upper;
finish reg;

/* Routine to test a linear combination of the estimates */
/* given L, this routine tests hypothesis that LB = 0. */
start test;
   dfn=nrow(L);
   Lb=L*b;
   vLb=L*xpxi*L`;
   q=Lb`*inv(vLb)*Lb /dfn;
   f=q/mse;
   prob=1-probf(f,dfn,dfe);
print ,f dfn dfe prob;
finish test;
/* Run it on population of U.S. for decades beginning 1790 */
x= { 1 1 1,
     1 2 4,
     1 3 9,
     1 4 16,
     1 5 25,
     1 6 36,
     1 7 49,
     1 8 64 };
y= {3.929,5.308,7.239,9.638,12.866,17.069,23.191,31.443};
name={"Intercept", "Decade", "Decade**2" };
tval=2.57; /* for 5 df at 0.025 level to get 95% conf. int. */
reset fw=7;
run reg;
do;
print ,"TEST Coef for Linear";
L={0 1 0 };
run test;
print ,"TEST Coef for Linear,Quad";
L={0 1 0,0 0 1};
run test;
print ,"TEST Linear+Quad = 0";
L={0 1 1 };
run test;
end;

线性代数知识梳理
1.Y=XB+E,参数B结果:

PROC IML;
   x={1 2,2 3,4 3};
   y={-1,-2,-2};
   b=INV(x‘*x)*x‘*y;
   PRINT b;
QUIT;

2.方差的最小二乘估计:

PROC IML;
    x={1 2,2 3,4 3};
    y={-1,-2,-2};
   n=NROW(x); /*样本容量n*/
   p=ROUND(TRACE(GINV(x)*x)); /*矩阵稚次*/
   in=I(n);
   s2=(1/(n-p))*y‘*(in-x*GINV(x‘*x)*x‘)*y;
   s=SQRT(s2);
   PRINT s2 s;
QUIT;

3.矩阵的分解:a.奇异阵分解:,CALL SVD( )

Proc IML;
   a = {1 2 3 4,
           5 6 7 8,
           9 0 1 2};
   RESET FUZZ;
   CALL SVD(u,q,v,a);
   upu = u‘*u;
   vpv = v‘*v;
   vvp = v*v‘;
   uup = u*u‘;
   a2 = u*DIAG(q)*v‘;
   ginva=GINV(a);
   m = NROW(a); n = NCOL(a);
   DO i=1 TO n;
      IF q[i] <= 1E-12 * q[1] then q[i] = 0;
      ELSE q[i] = 1 / q[i];
   END;
   ga = v*DIAG(q)*u‘;
   PRINT a a2„ u q v„ upu uup„ vpv vvp„ ginva ga;
QUIT;

b.)施密特正交化,CALL GSORTH()

PROC IML;
   a={1 2 3,2 3 5,3 7 11};
   CALL GSORTH(p,t,lindep,a);
   p_inv=INV(p);
   b=p*t;
   PRINT a b,p p_inv,t lindep;
QUIT;

4.线性方程求解

• HOMOGEN Function: 解其次线性方程
• SOLVE Function: 解普通线性方程
• TRISOLV Function: 三角形矩阵解方程

a).Solving AX=0 for X:

PROC IML;
   RESET FUZZ;
   a={10 5 15, 12 6 18, 14 7 21, 16 8 24};
   x=HOMOGEN(a);
   zero1=a*x[,1];
   zero2=a*x[,2];
   PRINT a x zero1 zero2;
QUIT;

b).Solving AX=B for X: 区别:x=INV(a)*b;

PROC IML;
   a={1 2 3,6 5 4,0 7 8};
   b={1, 2, 3};
   x=SOLVE(a,b);
   ax=a*x;
   PRINT a x b ax ;
QUIT;

5.特征值与特征向量

• EIGEN CALL计算特征值和特征向量为方阵

• EIGVAL Function计算特征值为一个方阵

• EIGVEC Function计算特征向量为一个方阵
• GENEIG Call 计算特征值和特征向量为一般矩阵

Proc IML;
   a={ 1 0 1,
          0 4 -1,
          1 -1 2};
   RESET FUZZ;
   eigval=EIGVAL(a);
   rank=ROUND(TRACE(GINV(a)*a)); /** 矩阵A的稚 **/
   PRINT a , ,eigval, ,rank;

   RESET FUZZ;
   CALL EIGEN(eigval,eigvec,a);
   ax1=a*eigvec[,1];
   ax2=a*eigvec[,2];
   ax3=a*eigvec[,3];
   lx1=eigval[1]*eigvec[,1];
   lx2=eigval[2]*eigvec[,2];
   lx3=eigval[3]*eigvec[,3];
   PRINT a„eigval„eigvec„ax1 lx1„ax2 lx2„ax3 lx3;
QUIT;

抱歉!评论已关闭.