我偶爾會去逛一個叫做「統計生活館」的網站。這個大概是我覺得國內最完整的一個在介紹統計學的網站。裡面有一個討論區,經常有人會上去問些問題。如果我知道答案,也會在那邊幫點小忙。
前幾日有位仁兄問了一個問題。他用亂數生成的方法模擬了一個資料庫包含1068筆資料,然後用 SAS 的 PROC CLUSTER 做 Cluster analysis,最後再把樹狀圖畫出來。這個流程其實是很簡單的,但是他要重複上述的動作 10000 次。他問程式該怎樣寫。
我本來想說用個 Macro 來搞定,但是卡在怎樣在亂數模擬變數那一階段請 SAS 生出 10000 個資料庫,因為 Do loop 不能掛在 Data procedure 的外面。後來想到一個簡單的方法,不用 Macro 就可以搞定。
原始程式如下:
data simulation;
do no=0001 to 1068;
array Q[*] Q01-Q25;
do j=1 to 25;
k=rand('uniform');
Q(j)=ranbin(1,1,k );
end;
output ;
end;
proc distance data=simulation method=djaccard absent=0 out=distjacc;
var anominal(Q01-Q25);
run;
proc cluster data=distjacc method=ward rsq rmsstd pseudo outtree=tree;
var _numeric_;
run;
proc tree data=tree out=out;
run;
修改後的程式如下:
data simulation;
do index=1 to 10000;
do no=0001 to 1068;
array Q[*] Q01-Q25;
do j=1 to 25;
k=rand('uniform');
Q(j)=ranbin(1,1,k );
end;
output ;
end;
end;
run;
proc distance data=simulation method=djaccard absent=0 out=distjacc;
by index;
var anominal(Q01-Q25);
run;
proc cluster data=distjacc method=ward rsq rmsstd pseudo outtree=tree;
by index;
var Dist1-Dist1068;
run;
proc tree data=tree;
by index;
run;
差別主要在生成資料那一步驟,只需要在外面多掛一個迴圈,讓裡面生亂數的迴圈再跑 10000 次,然後資料裡面會多了一個叫做 index 的變數。這個變數從 1 到 10000,就等於是我把 10000 個資料庫全部放在一起。之後的任何 Procedure 就可以用 by statement 來搞定。稍微需要注意的是,在 PROC CLUSTER 的 Var statement 中不能再用通用變數 _numeric_,而要改成Dist1-Dist1068。原因是若用 _numeric_ 通用變數的話,會把 index 也放進去分析,這樣會導致距離矩陣不是 n x n,而是 (n+1) x n。目前唯一還不算解決的是最後的 Proc tree 程序,不能再使用 out option,我想這是可以解決的,只是懶得再想而已。:P
-----
留言列表