SAS中文论坛

标题: 如何用SAS生成一阶马尔科夫转移矩阵？ [打印本页]

---

作者: shiyiming    时间: 2013-5-14 16:09
标题: 如何用SAS生成一阶马尔科夫转移矩阵？
原始数据如下：
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
因为在32个数据中，-1有10个，0有11个，1有11个且以0结尾。又因-1 →－1有5次；－1→0有3次；－1→1有2次；0→－1有2次；0 →0有5次；0 →1有3次，1 →－1有3次，1 →0有3次，1 →1有5次，所以Xn构成一个齐次马尔科夫链，且它的一步转移概率矩阵为：
1/2 3/10 1/5
2/11 5/11 4/11
3/11 3/11 5/11

如何用SAS生成上面的矩阵？

---

作者: shiyiming    时间: 2013-5-15 03:56
标题: 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]

---

欢迎光临 SAS中文论坛 (http://www.mysas.net/forum/)

Powered by Discuz! X3.2