2元分割表HGMの関数説明書について
******************

この説明書では HGM(holonomic gradient method) を用いた2元分割表の関数について説明する. ChangeLog の項目は www.openxm.org の cvsweb でソースコードを読む時の助けになる情報が書かれている. このパッケージは下記のようにロードする.

     load("gtt_ekn3.rr");

gtt_ekn3.rr は gtt_ekn.rr を置き換える大きく改良されたパッケージである. 以下のモジュール名 gtt_ekn はすべて gtt_ekn3 と読み替えてほしい.

最新版の asir-contrib package を取得するには, 下記のように更新関数を呼び出す.

     import("names.rr");
     asir_contrib_update(|update=1);

本文中で引用している文献を列挙する.

  * [GM2016] Y.Goto, K.Matsumoto, Pfaffian equations and contiguity relations of the hypergeometric function of type (k+1,k+n+2) and their applications, arxiv:1602.01637 (version 1)

  * [T2016] Y.Tachibana, 差分ホロノミック勾配法のモジュラーメソッドによる計算の高速化, 2016, 神戸大学修士論文.

  * [GTT2016] Y.Goto, Y.Tachibana, N.Takayama, 2元分割表に対する差分ホロノミック勾配法の実装, 数理研講究録.

  * [TGKT] Y.Tachibana, Y.Goto, T.Koyama, N.Takayama, Holonomic Gradient Method for Two Way Contingency Tables, arxiv:1803.04170

  * [TKT2015] N.Takayama, S.Kuriki, A.Takemura, $A$-hypergeometric distributions and Newton polytopes. arxiv:1510.02269

このマニュアルで説明する関数を用いたプログラム例は gtt_ekn/test-t1.rr など.

2元分割表HGMの関数
***********

超幾何関数E(k,n)
===========

gtt_ekn.gmvector
----------------

gtt_ekn.gmvector(BETA,P)
     :: 周辺和 BETA, セルの確率 P の二元分割表に付随する超幾何関数 E(k,n) の値およびその微分の値を戻す.

gtt_ekn.ekn_cBasis_2(BETA,P)
     の別名である.

return
     ベクトル, 超幾何関数の値とその微分. 詳しくは下記.

beta
     行和, 列和のリスト. 成分はすべて正であること.

p
     二元分割表のセルの確率のリスト

  * gmvector は Gauss-Manin vector の略である [GM2016].

  * gmvector の戻り値は [GM2016] の 6章 p.23 のベクトル Sである. これは [GM2016] の４章で定義されているベクトル F の定数倍であり, その定数は第一成分が [GM2016] の６章で定義されている級数 S の値と等しくなるように決められている.

  * r1 x r2 分割表を考える. m+1=r1, n+1=r2 とおく. 正規化定数 Z は分割表 u を (m+1) × (n+1) 行列とするとき p^u/u! の和である. ここで和は行和列和が BETA であるような u 全体でとる [TKT2015], [GM2016]. S はこの多項式 Z の p を
       [[1,y11,...,y1n],
        [1,y21,...,y2n],...,
        [1,ym1, ...,ymn],
        [1,1, ..., 1]]
     　(1 が L 字型に並ぶ), と正規化した級数である.

  * 2x(n+1)分割表で, gmvector の戻り値を Lauricella  F_D で書くことが以下のようにできる (b[2][1]-b[1][1] >= 0 の場合). ここで b[1][1], b[1][2] は, それぞれ 1 行目の行和, 2 行目の行和, b[2][i] は i 列目の列和である.
     S=F_D(-b[1,1], [-b[2,2],...,-b[2,n+1]], b[2,1]-b[1,1]+1 ; y)/C,
     C=b[1,1]! b[2,2]! ... b[2,n+1]! (b[2,1]-b[1,1])! とおく. 1/C は L 字型の分割表
     [[b[1,1],       0,      ..., 0       ],
      [b[2,1]-b[1,1],b[2,2], ..., b[2,n+1]]]
     に対応. gmvector は
     [S,(y11/a2) d_11 S,(y12/a3) d_12 S, ..., (y1n/a_(n+1)) d_1n S]
     である. ここで d_ij は yij についての微分,
       [a0,     a1, ...                      ,a_(n+2)]
     = [-b[1,2],-b[1,1],b[2,2], ..., b[2,n+1],b[2,1]]
     である.

  * 周辺和 BETAの時の正規化定数のセル確率 P に対する値は 多項式に退化した E(k,n) の値で表現できる. 文献 [TKT2015], [GM2016] 参照.

  * option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう [T2016]. 分散計算用の各種パラメータの設定は gtt_ekn.setup で行なう.

例: 次は2 x 2 分割表で行和が [5,1],  列和が [3,3], 各セルの確率が [[1/2,1/3],[1/7,1/5]] の場合の gmvector の値である.

     [3000] load("gtt_ekn.rr");
     [3001] ekn_gtt.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])
     [775/27783]
     [200/9261]

例: N を2以上の自然数とする時, Gauss の超幾何関数(この場合は多項式となる) F(-36N,-11N,2N,(1-1/N)/56) の値は T3 に代入される ( [TGKT] ).

     N=2;
     T2=gtt_ekn.gmvector([[36*N,13*N-1],[38*N-1,11*N]],[[1,(1-1/N)/56],[1,1]])[0][0];
     D=fac(36*N)*fac(11*N)*fac(2*N-1);
     T3=T2*D;

ちなみに同じ値を Mathematica に計算させるには

     n=2; Hypergeometric2F1[-36*n,-11*n,2*n,(1-1/n)/56]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な計算ができる. 守備範囲の異なるプログラム同士の比較, debug 用参考.

     [3080] import("tk_fd.rr");
     [3081] A=tk_fd.marginal2abc([4,5],[2,4,3]);
     [-4,[-4,-3],-1]  // 2変数 FD のパラメータ. a,[b1,b2],c
     [3082] tk_fd.fd_hessian2(A[0],A[1],A[2],[1/2,1/3]);
     Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
     [4483/124416,[ 1961/15552 185/1728 ],
      [ 79/288 259/864 ]
      [ 259/864 47/288 ]]
     // 戻値は [F=F_D, gradient(F), Hessian(F)]

     // ekn_gt での例と同じパラメータ.
     [3543] A=tk_fd.marginal2abc([5,1],[3,3]);
     [-5,[-3],-1]
     [3544] tk_fd.fd_hessian2(A[0],A[1],A[2],[(1/3)*(1/7)/((1/2)*(1/5))]);
     Computing Dmat(ca) for parameters B=[-3],X=[ 10/21 ]
     [775/27783,[ 20/147 ],[ 17/42 ]]

参考: 一般の A 分布の正規化定数についての Hessian の計算は実験的 package ot_hessian_ahg.rr で実装のテストがされている. (これはまだ未完成のテスト版なので出力形式等も将来的には変更される.)

     import("ot_hgm_ahg.rr");
     import("ot_hessian_ahg.rr");
     def  htest4() {
       extern C11_A;
       extern C11_Beta;
       Hess=newmat(7,7);
       A =C11_A;
       Beta0= [b0,b1,b2,b3];
       BaseIdx=[4,5,6];
       X=[x0,x1,x2,x3,x4,x5,x6];
       for (I=0; I<7; I++) for (J=0; J<7; J++) {
         Idx = [I,J];
         H=hessian_simplify(A,Beta0,X,BaseIdx,Idx);
         Hess[I][J]=H;
         printf("[I,J]=%a, Hessian_ij=%a\n",Idx,H);
       }
       return(Hess);
     }
     [2917] C11_A;
     [[0,0,0,1,1,1,1],[1,0,0,1,0,1,0],[0,1,1,0,1,0,1],[1,1,0,1,1,0,0]]
     [2918] C11_Beta;
     [166,36,290,214]
     [2919] Ans=htest4$
     [2920] Ans[0][0];
     [[((b1-b0-1)*x4)/(x0^2),[4]],[((b1-b0-1)*x6)/(x0^2),[6]],
      [(b1^2+(-2*b0-1)*b1+b0^2+b0)/(x0^2),[]],[(x6)/(x0),[6,0]],[(x4)/(x0),[4,0]]]

参照
     gtt_ekn.setup gtt_ekn.pfaffian_basis

ChangeLog

  * この関数は [GM2016] のアルゴリズムおよび [T2016] による modular method を用いた高速化を実装したものである.

  * 変更を受けたファイルは OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_pfaffian_8.rr

gtt_ekn.nc
----------

gtt_ekn.nc(BETA,P)
     :: 周辺和 BETA, セルの確率 P の二元分割表の条件付き確率の正規化定数 Z およびその微分の値を戻す.

return
     ベクトル [Z,[[d_11 Z, d_12 Z, ...], ..., [d_m1 Z, d_m2 Z, ...., d_mn Z]]]

beta
     行和, 列和のリスト. 成分はすべて正であること.

p
     二元分割表のセルの確率のリスト

  * r1 x r2 分割表を考える. m=r1, n=r2 とおく. 正規化定数 Z は分割表 u を m × n 行列とするとき p^u/u! の和である. ここで和は行和列和が BETA であるような u 全体でとる [TKT2015], [GM2016]. p^u は p_ij^u_ij の積, u! は u_ij! の積である. d_ij Z で Z の変数 p_ij についての偏微分を表す.

  * nc は gmvector の値を元に, [GM2016] の Prop 7.1 に基づいて Z の値を計算する.

  * option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう. 分散計算用の各種パラメータの設定は gtt_ekn.setup で行なう.

例: 2x3 分割表での Z とその微分の計算.

     [2237] gtt_ekn.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
     [4483/124416,[ 353/7776 1961/15552 185/1728 ]
     [ 553/20736 1261/15552 1001/13824 ]]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な計算ができる.

     [3076] import("tk_fd.rr");
     [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
     [-4,[-4,-3],-1]
     [3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
     RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
     [ 1 1 1 ]
     Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
     [4483/124416,[[353/7776,1961/15552,185/1728],
                   [553/20736,1261/15552,1001/13824]]]
     // 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],
     //             [d_21 Z, d_22 Z, d_23 Z]]] の値.
     //           ここで d_ij は i,j 成分についての微分を表す.

参照
     gtt_ekn.setup gtt_ekn.lognc

ChangeLog

  * 変更を受けたファイルは OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_eval.rr

gtt_ekn.lognc
-------------

gtt_ekn.lognc(BETA,P)
     :: 周辺和 BETA, セルの確率 P の二元分割表の条件付き確率の正規化定数 Z の log の近似値およびその微分の近似値を戻す.

return
     ベクトル [log(Z), [[d_11 log(Z), d_12 log(Z), ...], [d_21 log(Z),...], ... ]

beta
     行和, 列和のリスト. 成分はすべて正であること.

p
     二元分割表のセルの確率のリスト

  * 条件付き最尤推定に利用する [TKT2015].

  * option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう. 分散計算用の各種パラメータの設定は gtt_ekn.setup で行なう.

例: 2 × 3 分割表での例. 第一成分のみ近似値.

     [2238] gtt_ekn.lognc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
     [-3.32333832422461674630,[ 5648/4483 15688/4483 13320/4483 ]
     [ 3318/4483 10088/4483 9009/4483 ]]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な計算ができる.

     [3076] import("tk_fd.rr");
     [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
     [-4,[-4,-3],-1]
     [3078] tk_fd.log_ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
     RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
     [ 1 1 1 ]
     Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
     [-3.32333832422461674639485797719209322217260539267246045320,
      [[1.2598706, 3.499442, 2.971224],
       [0.7401293, 2.250278, 2.009591]]]
     // 戻値は [log(Z),
     //          [[d_11 log(Z), d_12 log(Z), d_13 log(Z)],
     //           [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]]
     // の近似値.

参照
     gtt_ekn.setup gtt_ekn.nc

ChangeLog

  * 変更を受けたファイルは OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.

gtt_ekn.expectation
-------------------

gtt_ekn.expectation(BETA,P)
     :: 周辺和 BETA, セルの確率 P の二元分割表の期待値を計算する.

return
     二元分割表の各セルの期待値のリスト.

beta
     行和, 列和のリスト. 成分はすべて正であること.

p
     二元分割表のセルの確率のリスト

  * [GM2016] の Algorithm 7.8 の実装.

  * option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう. 分散計算用の各種パラメータの設定は gtt_ekn.setup で行なう.

  * option index を与えると, 指定された成分の期待値のみ計算する. たとえば 2 x 2 分割表で index=[[0,0],[1,1]] と指定すると, 1 のある成分の期待値のみ計算する.

2×2, 3×3 の分割表の期待値計算例.

     [2235] gtt_ekn.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);
     [ 2/3 1/3 ]
     [ 4/3 8/3 ]
     [2236] gtt_ekn.expectation([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
     [ 5648/4483 7844/4483 4440/4483 ]
     [ 3318/4483 10088/4483 9009/4483 ]

     [2442] gtt_ekn.expectation([[4,14,9],[11,6,10]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]]);
     [ 207017568232262040/147000422096729819 163140751505489940/147000422096729819
                                             217843368649167296/147000422096729819 ]
     [ 1185482401011137878/147000422096729819 358095302885438604/147000422096729819
                                              514428205457640984/147000422096729819 ]
     [ 224504673820628091/147000422096729819 360766478189450370/147000422096729819
                                             737732646860489910/147000422096729819 ]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な計算ができる.

     [3076] import("tk_fd.rr");
     [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
     [-4,[-4,-3],-1]
     [3078] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
     RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
     [ 1 1 1 ]
     Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
     [[5648/4483,7844/4483,4440/4483],
      [3318/4483,10088/4483,9009/4483]]
     // 各セルの期待値.

参考: 一般の A 分布の計算は ot_hgm_ahg.rr. まだ実験的なため, module 化されていない. ot_hgm_ahg.rr についての参考文献: K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Systems II --- Holonomic Gradient Method, arxiv:1505.02947

     [3237] import("ot_hgm_ahg.rr");
     // 2 x 2 分割表.
     [3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],
             [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
     oohg_native=0, oohg_curl=1
     [1376777025/625400597,1750225960/625400597,
      2375626557/625400597,3252978816/625400597]
     // 2 x 2 分割表の期待値.

     // 2 x 3 分割表.
     [3238] hgm_ahg_expected_values_contiguity(
      [[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
      [5,2,4,3],[1,1/2,1/3,1,1,1],[x1,x2,x3,x4,x5,x6]|geometric=1);
     [5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
     // 2 x 3 分割表の期待値. 上と同じ問題.

3 x 3 分割表. 構造的0が一つ.

     /*
       dojo, p.221 のデータ.  成績3以下の生徒は集めてひとつに.
       2 1 1
       8 3 3
       0 2 6

       row sum: 4,14,8
       column sum: 10,6,10
       0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
     */

     A=[[0,0,0,1,1,1, 0,0],
        [0,0,0,0,0,0, 1,1],
        [1,0,0,1,0,0, 0,0],
        [0,1,0,0,1,0, 1,0],
        [0,0,1,0,0,1, 0,1]];
     B=[14,8,10,6,10];
     hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1],
     　　　　　　　　　　　　　　　　[x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);

     // 答.
     [14449864949304/9556267369631,
      10262588586540/9556267369631, 13512615942680/9556267369631,
      81112808747006/9556267369631,
      21816297744346/9556267369631, 30858636683482/9556267369631,

      25258717886900/9556267369631,51191421070148/9556267369631]

3 x 3 分割表.

     /*
      上のデータで 0 を 1 に変更.
       2 1 1
       8 3 3
       1 2 6

       row sum: 4,14,9
       column sum: 11,6,10
     */
     A=[[0,0,0,1,1,1,0,0,0],
        [0,0,0,0,0,0,1,1,1],
        [1,0,0,1,0,0,1,0,0],
        [0,1,0,0,1,0,0,1,0],
        [0,0,1,0,0,1,0,0,1]];
     B=[14,9,11,6,10];
     hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1,1],
                                   [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);

     // 期待値, 答.   x9 を指定していないので, 9番目の期待値は出力してない.
     [207017568232262040/147000422096729819,
      163140751505489940/147000422096729819,217843368649167296/147000422096729819,
      1185482401011137878/147000422096729819,
      358095302885438604/147000422096729819,514428205457640984/147000422096729819,
      224504673820628091/147000422096729819,360766478189450370/147000422096729819]

     // Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
     // まだ書いてない.

参照
     gtt_ekn.setup gtt_ekn.nc

ChangeLog

  * 変更を受けたファイルは OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.

gtt_ekn.setup
-------------

gtt_ekn.setup()
     :: 分散計算用の環境設定をおこなう. 現在の環境を報告する.

return

  * 使用するプロセスと素数の個数, 最小の素数を表示する. 準備されていない場合はその旨を表示.

  * このパッケージでの分散計算は複数のcpuを搭載した計算機で実行されることを想定している.

  * option nps (または number_of_processes)を与えると指定した数だけプロセスを用意する.

  * option nprm (または number_of_primes)を与えるとnprmが文字列の場合指定された素数リストのファイルを読み込む. nprmが自然数の場合さらにoption minp (minp =MINimum Prime)を与えるとminpより大きな素数をnprm個生成する. その際option fgp (または file_of_generated_primes)を与えると生成した素数リストをファイル名をfgpとして保存する.

  * 上記のoption を指定しなかった場合次のデフォルト値が用いられる. nps=1. nprm=10. fgp=0.

  * option report=1を与えると現在の環境の報告のみを行う. setup(|report=1)の別名としてreport関数を使用することもできる.

  * option subprogs=[file1,file2,...] により分散計算の子供プロセスにロードすべきファイル file1, file2, ... を指定する. default は subprogs=["gtt_ekn/childprocess.rr"] である.

  * gtt_ekn.set_debug_level(Mode) で Ekn_debug の値を設定する.

例: 素数のリストを生成してファイル p.txt へ書き出す.

     gtt_ekn.setup(|nps=2,nprm=20,minp=10^10,fgp="p.txt")$

例: chinese remainder theorem (crt) を使って gmvector を計算.

     [2867] gtt_ekn.setup(|nprm=20,minp=10^20);
     [2868] N=2; T2=gtt_ekn.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                                     [[1,(1-1/N)/56],[1,1]] | crt=1)$

参照
     gtt_ekn.nc gtt_ekn.gmvector

ChangeLog

  * 変更を受けたファイルは OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/g_mat_fac.rr

gtt_ekn.upAlpha, gtt_ekn.downAlpha
----------------------------------

gtt_ekn.upAlpha(I,K,N)

gtt_ekn.downAlpha(I,K,N)
     ::

i  a_i を a_i+1 (a_i を a_i-1) と変化させる contiguity relation.

k  E(k+1,n+k+2)型の超幾何関数の k. 分割表では (k+1)×(n+1).

n  E(k+1,n+k+2)型の超幾何関数の n. 分割表では (k+1)×(n+1).

return  contiguity relation の pfaffian_basis についての行列表現を戻す. [GM2016] の Cor 6.3.

  * upAlpha は　[GM2016] の Cor 6.3 の行列 U_i を戻す.

  * 関連する各関数の簡潔な説明と例も加える.

  * a_i を a_i-1 と変化させたい場合は関数 downAlpha を用いる.

  * a_i と分割表の周辺和を見るには, 関数 marginaltoAlpha([行和,列和]) を用いる.

  * pfaffian_basis は [GM2016] の４章のベクトル F に対応する偏微分を戻す.

  * optional 引数 arule, xrule で a_i や x_i_j を数にしたものをより効率的に求めることができる. 変化をうけるパラメータを数にしてしまっても特にエラー表示はしない. a_0 で和の条件を調整しているので注意(Todo, double check).

例: 以下の例は 2×2分割表(E(2,4)), 2×3分割表(E(2,5))の場合である. [2225] までは出力を略している.

     [2221] gtt_ekn.marginaltoAlpha([[1,4],[2,3]]);
     [[a_0,-4],[a_1,-1],[a_2,3],[a_3,2]]
     [2222] gtt_ekn.upAlpha(1,1,1);  // E(2,4) の a_1 方向の
                                     //     contiguity を表現する行列
     [2223] gtt_ekn.upAlpha(2,1,1);  // E(2,4) の a_2 方向
     [2224] gtt_ekn.upAlpha(3,1,1);  // E(2,4) の a_3 方向
     [2225] function f(x_1_1);
     [2232] gtt_ekn.pfaffian_basis(f(x_1_1),1,1);
     [ f(x_1_1) ]
     [ (f{1}(x_1_1)*x_1_1)/(a_2) ]
     [2233] function f(x_1_1,x_1_2);
     f() redefined.
     [2234] gtt_ekn.pfaffian_basis(f(x_1_1,x_1_2),1,2); // E(2,5), 2*3 分割表
     [ f(x_1_1,x_1_2) ]
     [ (f{1,0}(x_1_1,x_1_2)*x_1_1)/(a_2) ]
     [ (f{0,1}(x_1_1,x_1_2)*x_1_2)/(a_3) ]

     [2235]   RuleA=[[a_2,1/3],[a_3,1/2]]$ RuleX=[[x_1_1,1/5]]$
       base_replace(gtt_ekn.upAlpha(1,1,1),append(RuleA,RuleX))
      -gtt_ekn.upAlpha(1,1,1 | arule=RuleA, xrule=RuleX);

     [ 0 0 ]
     [ 0 0 ]

参照
     gtt_ekn.nc gtt_ekn.gmvector

ChangeLog

  * この関数は [GM2016] で与えられたアルゴリズムに従い contiguity relation を導出する.

  * 変更を受けたファイルは OpenXM/src/asir-contrib/packages/src/gtt_ekn/ekn_pfaffian_8.rr 1.1.

gtt_ekn.cmle
------------

gtt_ekn.cmle(U) u を観測データとするとき, P(U=u | row sum, column sum = these of U) を最大化する, 各セルの確率の近似値を求める.
     ::

u  観測データ(分割表)

return  セルの確率(分割表形式)

  * u を観測データとするとき, P(U=u | row sum, column sum = these of U) を最大化する, 各セルの確率の近似値を求める.

  * optional parameter で algorithm の振る舞い(たとえば有理数を近似して, 分母分子が小さい有理数にする, gradient descent の step幅)を調整すべきだが, これは作業中. 2017.03.03

例: 2 x 4 分割表.

     U=[[1,1,2,3],[1,3,1,1]];
     gtt_ekn.cmle(U);
      [[ 1 1 2 3 ]
       [ 1 3 1 1 ],[[7,6],[2,4,3,4]],   // Data, row sum, column sum
      [ 1 67147/183792 120403/64148 48801/17869 ]  // probability obtained.
      [ 1 1 1 1 ]]

例: 上の例は次の関数に.

     gtt_ekn.cmle_test3();

参照
     gtt_ekn.expectation

ChangeLog

  * gtt_ekn/mle.rr に本体がある.

  * gtt_ekn.rr の cmle 関数は wrapper.

gtt_ekn.set_debug_level, gtt_ekn.show_path, gtt_ekn.get_svalue, gtt_ekn.assert1, gtt_ekn.assert2, gtt_ekn.prob1
---------------------------------------------------------------------------------------------------------------

gtt_ekn.set_debug_level(M) debug メッセージのレベルを設定.

gtt_ekn.contiguity_mat_list_2  使用する contiguity を構成.

gtt_ekn.show_path()  どのように contiguity を適用したかの情報.

gtt_ekn.get_svalue()  static 変数の値を得る.

gtt_ekn.assert1(N)  2x2 分割表で動作チェック.

gtt_ekn.assert2(N)  3x3 分割表で動作チェック.

gtt_ekn.prob1(R1,R2,SIZE)  R1 x R2 分割表用のテストデータを作る.
     ::

m  レベル.

  * (M & 0x1) == 0x1 の時 g_mat_fac_test_plain と g_mat_fac_itor の両方を呼び出し値を比較する (gtt_ekn.setup した状態で).

  * (M & 0x2) == 0x2 の時 g_mat_fac_test への引数を tmp-input-数.ab として保存.

  * (M & 0x4) == 0x4 の時 matrix factorial の計算の呼び出し引数を表示.

  * N は問題の周辺和のサイズ.

  * get_svalue の戻り値は [Ekn_plist,Ekn_IDL,Ekn_debug,Ekn_mesg,XRule,ARule,Verbose,Ekn_Rq] の値.

例.

     [2846] gtt_ekn.set_debug_level(0x4);
     [2847] N=2; T2=gtt_ekn.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                                     [[1,(1-1/N)/56],[1,1]])$
     [2848] level&0x4: g_mat_fac_test([ 113/112 ]
     [ 1/112 ],[ (t+225/112)/(t^2+4*t+4) (111/112*t+111/112)/(t^2+4*t+4) ]
     [ (1/112)/(t^2+4*t+4) (111/112*t+111/112)/(t^2+4*t+4) ],0,20,1,t)
     Note: we do not use g_mat_fac_itor. Call gtt_ekn.setup(); to use the crt option.
     level&0x4: g_mat_fac_test([ 67/62944040755546030080000 ]
     [ 1/125888081511092060160000 ],[ (t+24)/(t^2+25*t+46) (2442)/(t^2+25*t+46) ]
     [ (1)/(t^2+25*t+46) (-111*t-111)/(t^2+25*t+46) ],0,73,1,t)
     level&0x4: g_mat_fac_test ------  snip

例.

     [2659] gtt_ekn.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]])$
     [2660] L=matrix_transpose(gtt_ekn.show_path())$
     [2661] L[2];
     [1 4 3 2]

[1 4 3 2] の index をもつパラメーター alpha の方向の contigity を求めそれを掛けて計算したことがわかる.  L[0] は用いた contiguity の行列. L[1] はcontiguity を適用する step 数.

例. 値を計算せずに path のみ求めたい場合.

     A=gtt_ekn.marginaltoAlpha_list([[400,410,1011],[910,411,500]])$
     [2666] gtt_ekn.contiguity_mat_list_2(A,2,2)$
     [2667] L=matrix_transpose(gtt_ekn.show_path())$
     [2668] L[2];
     [ 2 1 5 4 3 ]

例. 値を計算せずに path のみ求めたい場合. gtt_ekn3 による新しいアルゴリズムによる path の表示.

     A=gtt_ekn3.marginaltoAlpha_list([[10,20],[15,15]])$
     [2666] gtt_ekn3.contiguity_mat_list_3(A,1,1 | xrule=[[x_1_1,1/2]])$
     [t,[[ (-t-43/2)/(t-2) (-15/2)/(t-2) ]
     [ 1/2 -1/2 ],-9]]

例. 0 が戻れば g_mat_fac_plain と指定した計算方法の結果が一致したことがわかる. option を書かないと g_mat_fac_int との比較となる.

     [8859] gtt_ekn.assert2(1);
     Marginal=[[130,170,353],[90,119,444]]
     P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]
     Try g_mat_fac_test_int: Note: we do not use g_mat_fac_itor. Call gtt_ekn.setup(); to use the crt option.
     Timing (int) =0.413916 (CPU) + 0.590723 (GC) = 1.00464 (total), real time=0.990672

     Try g_mat_fac_test_plain: Note: we do not use g_mat_fac_itor. Call gtt_ekn.setup(); to use the crt option.
     Timing (rational) =4.51349 (CPU) + 6.32174 (GC) = 10.8352 (total)
     diff of both method =
     [ 0 0 0 ]
     [ 0 0 0 ]
     [ 0 0 0 ]
     [8860]

     [8863] gtt_ekn.setup(|nprm=100,minp=10^50);
     Number of processes = 1.
     Number of primes = 100.
     Min of plist = 100000000000000000000000000000000000000000000000151.
     0
     [8864] gtt_ekn.assert2(1 | crt=1);
     Marginal=[[130,170,353],[90,119,444]]
     P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]
     Try [[crt,1]]
     ----  snip

なお二番目の例の timing (total) [例では省略] は mod 計算を subprocess がやっているので正しい値ではない. real time が計算時間の目安になる.

例.

     [9054] L=gtt_ekn.prob1(3,5,10 | factor=1, factor_row=3);
     [[[10,20,420],[30,60,90,120,150]],[[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1,1,1,1]]]
     [9055] number_eval(gtt_ekn.expectation(L[0],L[1]));
     [ 0.434161208918863  ... snip ]

参照
     gtt_ekn.nc

ChangeLog

  * gtt_ekn/ekn_eval.rr で matrix factorial の計算の呼び出し引数を表示する.

  * grep 'iand(Ekn_debug,0x1)' *.rr でソースコードの該当の位置をさがす.

modular計算
*********

中国剰余定理とitor
===========

gtt_ekn.chinese_itor
--------------------

gtt_ekn.chinese_itor(DATA,IDLIST)
     :: mod p で計算した結果(ベクトル)から chinese remainder theorem, itor(integer to rational) で有理数ベクトルを得る.

return  [val, n]  ここで val は答え. また, n = n1*n2*...

data   [[val1,n1],[val2,n2], ...], ここで val mod n1 = val1, val mod n2 = val2,...

idlist  chinese, itor を実行するサーバIDのリスト.

  * 中国剰余定理を用いて val0 mod n1 = val1, val0 mod n2 = val2, ... となる val0 を求める. val に algorithm itor を適用する.

  * sqrt(n) より val0 が大きい時は itor が適用されて val0 が有理数 val=a/b に変換される. つまり b*x =1 mod n となる逆数 x を考えて, x*a % n = val0 となる数 val を戻す. 見つからないときは failure を戻す.

例: [3!, 5^3*3!]=[6,750] が戻り値. 6 mod 109 =6, 750 mod 109=96 が最初の引数の [[6,96],109]. 以下同様.

     gtt_ekn.setup(|nps=2,nprm=3,minp=101,fgp="p_small.txt");
     SS=gtt_ekn.get_svalue();
     SS[0];
       [103,107,109]   // list of primes
     SS[1];
       [0,2]           // list of server ID's
     gtt_ekn.chinese_itor([[[ 6,96 ],109],[[ 6,29 ],103],[[ 6,1 ],107]],SS[1]);
       [[ 6 750 ],1201289]

     // 引数はスカラーでもよい.
     gtt_ekn.chinese_itor([[96,109],[29,103]],SS[1]);
       [[ 750 ],11227]

例: gtt_ekn/childprocess.rr (server で実行される) の関数 chinese (chinese remainder theorem) と euclid.

     load("gtt_ekn/childprocess.rr");
     chinese([newvect(2,[6,29]),103],[newvect(2,[6,750]),107*109]);
       // mod 103 で [6,29], mod (107*109) で [6,750] となる数を mod 103*(107*109)
       // で求めると,
       [[ 6 750 ],1201289]
     euclid(3,103);  // mod 103 での 3 の逆数. つまり 1/3
       -34
     3*(-34) % 103; // 確かに逆数.
        1

例: gtt_ekn/childprocess.rr (server で実行される) の関数 itor (integer to rational) の例. itor(Y,Q,Q2,Idx) では Y < Q2 なら Y がそのまま戻る.  Idx は 内部用の index で好きな数でよい. 戻り値の第２成分となる.

     load("gtt_ekn/childprocess.rr");
     for (I=1;I<11; I++) print([I,itor(I,11,3,0)]);
     [1,[1,0]]
     [2,[2,0]]
     [3,[-2/3,0]] //euclid(3,11); ->4,  4*(-2)%11 -> 3 なので確かに -2/3 は元の数の候補
     [4,[failure,0]]
     [5,[-1/2,0]]
     [6,[1/2,0]]
     [7,[-1/3,0]]
     [8,[failure,0]]
     [9,[-2,0]]
     [10,[-1,0]]

参照
     gtt_ekn.setup

ChangeLog

  * 関連ファイルは gtt_ekn/g_mat_fac.rr gtt_ekn/childprocess.rr

binary splitting
****************

matrix factorial
================

gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit
-----------------------------------------------------------------------

gtt_ekn3.init_bsplit(|minsize=16,levelmax=1);
     :: binary split の実行のためのパラメータを設定する.

gtt_ekn3.init_dm_bsplit(|bsplit_x=0, bsplit_reduce=0)
     :: binary split の分散実行のためのパラメータを設定する.

gtt_ekn3.setup_dm_bsplit(C)
     :: binary split の分散実行のために C 個のプロセスを立ち上げる.

C はlevelmax-1 に設定する. 特に levalmax=1 のときは分散計算を行わない.

bsplit_x=1 のとき, debug 用に各プロセスを xterm で表示.

  * expectation などの関数に bs=1 オプションを与えると matrix factorial を binary splitting method で計算する.

例: bs=1 と無い場合の比較.

     [4618] cputime(1)$
     [4619] gtt_ekn3.expectation(Marginal=[[1950,2550,5295],[1350,1785,6660]],
                               P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]|bs=1)$
     4.912sec(4.914sec)
     [4621] V2=gtt_ekn3.expectation(Marginal=[[1950,2550,5295],[1350,1785,6660]],
                               P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]])$
     6.752sec(6.756sec)

例: 分散計算する場合. 分散計算はかえって遅くなる場合が多いので注意. 下記の例での bsplit_x=1 option は debug windows を開くのでさらに遅くなる. gtt_ekn3.test_bs_dist(); でもテストできる.

     [3669] C=4$ gtt_ekn3.init_bsplit(|minsize=16,levelmax=C+1)$ gtt_ekn3.init_dm_bsplit(|bsplit_x=1)$
     [3670] [3671] [3672] gtt_ekn3.setup_dm_bsplit(C);
     [0,0]
     [3673] gtt_ekn3.assert2(10|bs=1)$

参照
     gtt_ekn3.gmvector gtt_ekn3.expectation gtt_ekn3.assert1 gtt_ekn3.assert2

ChangeLog

  * 関連ファイルは gtt_ekn3/mfac_bs.rr gtt_ekn3/dm_bsplit.rr