yydy331a.sas

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

/*  yydy331a.sas */


data dy323;
  input y1-y4 group @@;
  cards;
260  75 40 18 1 200  72 34 17 1 240  87 45 18 1 170  65 39 17 1
270 110 39 24 1 205 130 34 23 1 190  69 27 15 1 200  46 45 15 1
250 117 21 20 1 200 107 28 20 1 225 130 36 11 1 210 125 26 17 1
170  64 31 14 1 270  76 33 13 1 190  60 34 16 1 280  81 20 18 1
310 119 25 15 1 270  57 31  8 1 250  67 31 14 1 260 135 39 29 1
310 122 30 21 2 310  60 35 18 2 190  40 27 15 2 225  65 34 16 2
170  65 37 16 2 210  82 31 17 2 280  67 37 18 2 210  38 36 17 2
280  65 30 23 2 200  76 40 17 2 200  76 39 20 2 280  94 26 11 2
190  60 33 17 2 295  55 30 16 2 270 125 24 21 2 280 120 32 18 2
240  62 32 20 2 280  69 29 20 2 370  70 30 20 2 280  40 37 17 2
320  64 39 17 3 260  59 37 11 3 360  88 28 26 3 295 100 36 12 3
270  65 32 21 3 380 114 36 21 3 240  55 42 10 3 260  55 34 20 3
260 110 29 20 3 295  73 33 21 3 240 114 38 18 3 310 103 32 18 3
330 112 21 11 3 345 127 24 20 3 250  62 22 16 3 260  59 21 19 3
225 100 34 30 3 345 120 36 18 3 360 107 25 23 3 250 117 36 16 3
;
 proc iml;
   n1=20;n2=20;n3=20;
   n=n1+n2+n3;k=3; p=4;
   use dy323(obs=20);
   xa={y1 y2 y3 y4};
   read all  var xa into x1; print x1;
   use dy323(firstobs=21 obs=40);
   read all  var xa into x2; print x2;
   use dy323(firstobs=41 obs=60);
   read all  var xa into x3; print x3;
   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;
   da=det(a/(n-k));
   da1=det(a1/(n1-1));
   da2=det(a2/(n2-1));
   da3=det(a3/(n3-1));
   m=(n-k)*log(da)-19*(log(da1)+log(da2)+log(da3));
   dd=(2*p*p+3*p-1)*(k+1)/(6*(p+1)*(n-k));
   df=p*(p+1)*(k-1)/2;
   kc=(1-dd)*m;
   print da da1 da2 da3;
   print m dd df ;
   p0=1-probchi(kc,df);
   print kc p0;
 run;