yydy332.sas

高老师的经典多元统计分析教材源代码非常经典使用

/* yydy332.sas */

 proc iml;
   n1=20;n2=20;n3=20;
   n=n1+n2+n3;k=3; p=4;
   x1={
      260  75 40 18 , 200  72 34 17 , 240  87 45 18 , 170  65 39 17 ,
      270 110 39 24 , 205 130 34 23 , 190  69 27 15 , 200  46 45 15 ,
      250 117 21 20 , 200 107 28 20 , 225 130 36 11 , 210 125 26 17 ,
      170  64 31 14 , 270  76 33 13 , 190  60 34 16 , 280  81 20 18 ,
      310 119 25 15 , 270  57 31  8 , 250  67 31 14 , 260 135 39 29 } ;
    x2={
       310 122 30 21 ,310  60 35 18 ,190  40 27 15 ,225  65 34 16 ,
       170  65 37 16 ,210  82 31 17 ,280  67 37 18 ,210  38 36 17 ,
       280  65 30 23 ,200  76 40 17 ,200  76 39 20 ,280  94 26 11 ,
       190  60 33 17 ,295  55 30 16 ,270 125 24 21 ,280 120 32 18 ,
       240  62 32 20 ,280  69 29 20 ,370  70 30 20 ,280  40 37 17 };
    x3={
       320  64 39 17 ,260  59 37 11 ,360  88 28 26 ,295 100 36 12 ,
       270  65 32 21 ,380 114 36 21 ,240  55 42 10 ,260  55 34 20 ,
       260 110 29 20 ,295  73 33 21 ,240 114 38 18 ,310 103 32 18 ,
       330 112 21 11 ,345 127 24 20 ,250  62 22 16 ,260  59 21 19 ,
       225 100 34 30 ,345 120 36 18 ,360 107 25 23 ,250 117 36 16 };
     xx=x1//x2//x3;
   ln={[20] 1} ;
   x10=(ln*x1)/n1; print x10;
   mm1=i(n1)-j(n1,n1,1)/n1;
   mm=i(n)-j(n,n,1)/n;
   a1=x1`*mm1*x1; print a1;
   a2=x2`*mm1*x2; print a2;
   a3=x3`*mm1*x3; print a3;
   tt=xx`*mm*xx;  print tt;
   a=a1+a2+a3;
   print a;
   dt=det(tt/(n-k));
   da1=det(a1/(n1-1));
   da2=det(a2/(n2-1));
   da3=det(a3/(n3-1));
   m5=(n-k)*log(dt)-19*(log(da1)+log(da2)+log(da3));
   b=(2*p*p+3*p-1)*(k+1)/(6*(p+3)*(n-k))-(p-k+2)/((n-k)*(p+3));
   df=p*(p+3)*(k-1)/2;
   kc=(1-b)*m5;
   print dt da1 da2 da3;
   print m5 b df ;
   p0=1-probchi(kc,df);
   print kc p0;
 run;