SAS/IML过程，读取矩阵每行数据计算马氏距离问题

shiyiming

下面计算马氏距离的程序
data score;
   input num x1-x5@;
   label x1='政治' x2='语文' x3='外语' x4='数学' x5='物理';
cards;
1  99  94  93 100 100
2  99  88  96  99  97
3 100  98  81  96 100
4  93  88  88  99  96
5 100  91  72  96  78
6  90  78  82  75  97
7  75  73  88  97  89
8  93  84  83  68  88
9  87  73  60  76  84
10 95  82  90  62  39
11 76  72  43  67  78
12 85  75  50  34  37
;
proc iml;
    n=12;p=5;
        xx={x4 x5};
        use score;
        read all var xx into x;
        e={[12] 1};
        x0=(e*x)/n;
        mm=i(12)-j(12,12,1)/n;
        a=x`*mm*x;
        s=a/(n-1);
        si=inv(s);print x0 s si; /*si为s的逆矩阵*/
                 [color=#FF0000:35e970rp] use score(obs=1);
        read all var xx into xx1;
               d1=(xx1-x0)*si*(xx1-x0)`; /*d 为马氏距离*/
                use score(firstobs=2 obs=2);
        read all var xx into xx2;
        d2=(xx2-x0)*si*(xx2-x0)`;
                use score(firstobs=3 obs=3);
        read all var xx into xx3;
        d3=(xx3-x0)*si*(xx3-x0)`;
                use score(firstobs=4 obs=4);
         read all var xx into xx4;
        d4=(xx4-x0)*si*(xx4-x0)`;
               use score(firstobs=5 obs=5);
        read all var xx into xx5;
        d5=(xx3-x0)*si*(xx5-x0)`;
        use score(firstobs=6 obs=6);
        read all var xx into xx6;
        d6=(xx6-x0)*si*(xx6-x0)`;
        use score(firstobs=7 obs=7);
        read all var xx into xx7;
        d7=(xx7-x0)*si*(xx7-x0)`;
        use score(firstobs=8 obs=8);
        read all var xx into xx8;
        d8=(xx8-x0)*si*(xx8-x0)`;
        use score(firstobs=9 obs=9);
        read all var xx into xx9;
        d9=(xx9-x0)*si*(xx9-x0)`;
        use score(firstobs=10 obs=10);
        read all var xx into xx10;
        d10=(xx10-x0)*si*(xx10-x0)`;
        use score(firstobs=11 obs=11);
        read all var xx into xx11;
        d11=(xx11-x0)*si*(xx11-x0)`;
        use score(firstobs=12 obs=12);
        read all var xx into xx12;
        d12=(xx12-x0)*si*(xx12-x0)`;[/color:35e970rp]print d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12;
run;quit;
我不会直接从矩阵x中读取每行数据，只能用 use score(obs=1);read all var xx into xx1;d1=(xx1-x0)*si*(xx1-x0)`; 来计算在第一个观测（即x矩阵的第一行）下的马氏距离，共12个同样的步骤，很繁琐！请哪位高手给个循环模式，可直接一起计算这十二个马氏距离！谢谢！！

shiyiming

我不敢不懂装懂。不过好像在这里的IML并不是很复杂。把你红色的代码用以下的代替应该可以：

[code:w4zpmdad]_t =(x-x0)*si*(x-x0)`;
f =j(1, 12);
do i =1 to 12;
        f[1,i] =_t[i, i];
end;
print f;[/code:w4zpmdad]
另外
[quote:w4zpmdad]d5=(xx3-x0)*si*(xx5-x0)`;[/quote:w4zpmdad]
这里的xx3应该是错误的。
希望有所帮助。
京剧

shiyiming

我调试了，可结果出现这样   
                                   F
1      1              1          1             1            1              1             1            1           1            1          1
是不是那里出错了，请再修改下，非常感谢！

另外就是，在求出了d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 这些值后，怎样把它们存放在一个新的数据集中，以备接下其它过程用（如用freq过程分析这12个数据频率等）？

shiyiming

我不修改你修改。
问题很显然是你没有使用loop更新F矩阵。
最后的结果当然和你原来的是一模一样的。
京剧

shiyiming

[code:1io7xa07]proc iml;
n=12;p=5;
xx={x4 x5};
use score;
read all var xx into x;
e={[12] 1};
x0=(e*x)/n;
mm=i(12)-j(12,12,1)/n;
a=x`*mm*x;
s=a/(n-1);
si=inv(s);print x0 s si; /*si为s的逆矩阵*/

d=j(1,nrow(x),.);
do i=1 to nrow(x);
d[i]=(x[i,]-x0)*si*(x[i,]-x0)`;
end;
print d;
run;quit;
[/code:1io7xa07]