| 
 | 
沙发
 
 
 楼主 |
发表于 2013-5-15 03:56:43
|
只看该作者
 
 
 
Re: 如何用SAS生成一阶马尔科夫转移矩阵?
抛块砖头,机器上没有IML 
[code:21u2mewb]%let chain= 1,1,0,0,0 ,0 ,1 ,1,0,1,1,1,-1,-1,0,-1,1,1,-1,-1,-1,0,1,0,0,0 ,-1,1 ,-1,-1,-1,0 
; 
%let N = %sysfunc(countw(%superq(chain), %str(,))); 
proc fcmp; 
        array X[&N] (&chain); 
        array mat[3, 3]; 
        array J[3,1]; 
        array prob[3, 3]; 
         
        call zeromatrix(mat); 
        do i=1 to &N-1; 
                mat[X[i]+2, X[i+1]+2] = mat[X[i]+2, X[i+1]+2] + 1; 
        end; 
        mat[X[&N]+2, x[1]+2] = mat[X[&N]+2, x[1]+2] + 1; 
        put mat=; 
         
        call fillmatrix(J,1); 
        call mult(mat, J, J); 
 
        do i=1 to 3; 
                do k=1 to 3; 
                        prob[i,k] = mat[i,k]/J[i]; 
                end; 
        end;         
        put prob=; 
run; 
[/code:21u2mewb] 
 
[code:21u2mewb]%let circular = 1,1,0,0,0 ,0 ,1 ,1,0,1,1,1,-1,-1,0,-1,1,1,-1,-1,-1,0,1,0,0,0 ,-1,1 ,-1,-1,-1,0; 
%let N = %sysfunc(countw(%superq(circular), %str(,))); 
data xxx/view=xxx; 
        array x[0:%eval(&N-1)] _temporary_ (&circular); 
        do _n_=lbound(x) to hbound(x); 
                x1 = x[_n_]; 
                x2 = x[mod(_n_+1, &N)]; 
                output; 
        end; 
proc format; 
        picture mypct 
        low-high = '09.9999'(mult=100); 
 
proc tabulate data=xxx out=yyy; 
        class x1 x2; 
        table x1=' ', x2=' '*rowpctn=' '*f=mypct.; 
run; 
options validvarname=any; 
proc transpose data=yyy out=zzz(drop=_name_); 
        by x1; 
        id x2; 
        var pct:; 
run; 
data zzz; 
        set zzz; 
        array x[*] _numeric_; 
        do _n_=1 to dim(x); 
                if lowcase(vname(x[_n_])) ne 'x1' then 
                        x[_n_] = coalesce(x[_n_],0)/100; 
        end; 
run;[/code:21u2mewb] |   
 
 
 
 |