2元分割表HGM関数

Risa/Asir 2元分割表HGM関数説明書

1.1 版

2017 年 3 月 3 日

by Y.Goto, Y.Tachibana, N.Takayama

Copyright © Risa/Asir committers 2004–2010. All rights reserved.


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

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

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

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

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


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2 2元分割表HGMの関数


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

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


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2.1.1 gtt_ekn.gmvector

gtt_ekn.gmvector(beta,p)

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

gtt_ekn.ekn_cBasis_2(beta,p)

の別名である.

return

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

beta

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

p

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

例: 次は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]

参考: 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 @ref{gtt_ekn.pfaffian_basis}

ChangeLog


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2.1.2 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

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

例: 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


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2.1.3 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

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

例: 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


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2.1.4 gtt_ekn.expectation

gtt_ekn.expectation(beta,p)

:: 周辺和 beta, セルの確率 p の二元分割表の期待値を計算する.

return

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

beta

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

p

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

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


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2.1.5 gtt_ekn.setup

gtt_ekn.setup()

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

return

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

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

gtt_ekn.nc gtt_ekn.gmvector

ChangeLog


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2.1.6 gtt_ekn.upAlpha

gtt_ekn.upAlpha(i,k,n)

::

i 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.

例: 以下の例は 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) ]
参照

gtt_ekn.nc gtt_ekn.gmvector

ChangeLog


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

2.1.7 gtt_ekn.cmle

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

::

u 観測データ(分割表)
return セルの確率(分割表形式)

例: 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


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

3 modular計算


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

3.1 中国剰余定理とitor


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

3.1.1 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のリスト.

例: [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 で好きな数でよい. 戻り値の第2成分となる.

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


[ << ] [ < ] [?] [ > ] [ >> ]         [??] [??] [???] [ ? ]

Index

??:   G  
?????  ?

G
gtt_ekn.chinese_itor 中国剰余定理とitor 3.1.1 gtt_ekn.chinese_itor
gtt_ekn.cmle 2.1.7 gtt_ekn.cmle
gtt_ekn.expectation 2.1.4 gtt_ekn.expectation
gtt_ekn.gmvector 2.1.1 gtt_ekn.gmvector
gtt_ekn.lognc 2.1.3 gtt_ekn.lognc
gtt_ekn.nc 2.1.2 gtt_ekn.nc
gtt_ekn.setup 2.1.5 gtt_ekn.setup
gtt_ekn.upAlpha 2.1.6 gtt_ekn.upAlpha

??:   G  

[??] [??] [???] [ ? ]

??


[??] [??] [???] [ ? ]

???????


[??] [??] [???] [ ? ]

????????

?????10? 24, 2017?texi2html 5.0????????????

????????????????????????????

??? ?? ??? 1.2.3???????
[ << ] FastBack Beginning of this chapter or previous chapter 1
[ < ] Back Previous section in reading order 1.2.2
[?] Up Up section 1.2
[ > ] Forward Next section in reading order 1.2.4
[ >> ] FastForward Next chapter 2
[??] ?? Cover (top) of document  
[??] ?? Table of contents  
[???] ??? ???  
[ ? ] About About (help)  

?????????????????1.2.3???????????????


?????10? 24, 2017?texi2html 5.0????????????