遷移行列の冪乗根 の変更点

状態遷移確率行列を2乗することで,2時点先の状態遷移確率行列を得ます.ですから,1年間の状態遷移確率行列の12乗根を求めれば,1ヶ月間の状態遷移確率行列になります.
ですが,SASの行列演算言語IMLでは,行列の冪乗演算子**は,-1以上の整数しか受付けません.

そこで,以下に,行列を固有値と固有ベクトル行列に分解して,冪乗根を求める方法を示します.
ここでは,平方根を求めていますが,##(1/2)を,##(1/12)にすれば12乗根が得られます.
(注:##は,IMLでの,行列の「各要素」の冪乗演算子です.)

警告!:すみません.固有値が複素数になるケースに対応していませんでした.
そのうち,直しますが,どなたか直してくれるとうれしいです(翔)
→たねやさんが回答してくれました.(感謝感謝)

#pre{{
*************************************************
rootmatrix.sas
翔
*************************************************;

proc iml;

A={
 0.5 0.3  0.2
,0.1 0.85 0.05
, 0    0    1
};
print A;


/*
固有値は実数であるとみなしてますが,
複素数の場合もあり,その計算ができていません!
申し訳ない.
*/
D=diag(eigval(A)[,1]);
U=eigvec(A);

RA=U*(D##(1/2))*inv(U);
print RA;

/*確認のため,2乗してみる*/
RA2=RA**(-1);
print RA2;

run;

/*
              A

      0.5       0.3       0.2
      0.1      0.85      0.05
        0         0         1


              RA

0.6989223 0.1858028 0.1152748
0.0619343 0.9156923 0.0223734
        0         0         1


             RA2

      0.5       0.3       0.2
      0.1      0.85      0.05
        0         0         1
 */

}}

必ずしも解が見つかるとは限りませんが、その他にもIMLの非線形最適化サブルーチンを応用して利用することができます。最適化の結果で、目的関数が0に収束していることを確認してください。
#pre{{
*************************************************
rootmatrixnlp.sas
たねや
*************************************************;

proc iml;

A={
 0.5 0.3  0.2
,0.1 0.85 0.05
, 0    0    1
};



n=ncol(a);

start obj(x) global(a,n);
    f=(btran(x,1,n)**2-a)[##];
    return(f);
finish obj; 

x0 = btran(a,1,n);
optn = {0 1};
call nlptr(rc,xres,"obj",x0,optn);
RA=btran(xres,1,n);
print RA;

/*
ABSGCONV convergence criterion satisfied.


                                   0.6989223 0.1858028 0.1152748
                                   0.0619343 0.9156923 0.0223734
                                   -4.881E-9 -6.064E-9         1
 */

}}