[BACK]Return to gtt_ekn-ja.texi CVS log [TXT][DIR] Up to [local] / OpenXM / src / asir-contrib / packages / doc / gtt_ekn

Annotation of OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi, Revision 1.4

1.4     ! takayama    1: %% $OpenXM: OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi,v 1.3 2016/03/28 23:14:54 takayama Exp $
1.3       takayama    2: %% ptex -kanji euc gtt_ekn.texi   (.texi までつける. platex でなく ptex)
                      3: %% 以下コメントは @comment で始める.  \input texinfo 以降は普通の tex 命令は使えない.
1.1       takayama    4: \input texinfo
                      5: @iftex
                      6: @catcode`@#=6
                      7: @def@fref#1{@xrefX[#1,,@code{#1},,,]}
                      8: @def@b#1{{@bf@gt #1}}
                      9: @catcode`@#=@other
                     10: @end iftex
                     11: @overfullrule=0pt
                     12: @c -*-texinfo-*-
                     13: @comment %**start of header
1.3       takayama   14: @comment --- おまじない終り ---
1.1       takayama   15:
1.3       takayama   16: @comment --- GNU info ファイルの名前 ---
1.1       takayama   17: @setfilename xyzman
                     18:
1.3       takayama   19: @comment --- タイトル ---
                     20: @settitle 2元分割表HGM
1.1       takayama   21:
                     22: @comment %**end of header
                     23: @comment %@setchapternewpage odd
                     24:
1.3       takayama   25: @comment --- おまじない ---
1.1       takayama   26: @ifinfo
                     27: @macro fref{name}
                     28: @ref{\name\,,@code{\name\}}
                     29: @end macro
                     30: @end ifinfo
                     31:
                     32: @iftex
                     33: @comment @finalout
                     34: @end iftex
                     35:
                     36: @titlepage
1.3       takayama   37: @comment --- おまじない終り ---
1.1       takayama   38:
1.3       takayama   39: @comment --- タイトル, バージョン, 著者名, 著作権表示 ---
                     40: @title 2元分割表HGM関数
                     41: @subtitle Risa/Asir 2元分割表HGM関数説明書
1.4     ! takayama   42: @subtitle 1.1 版
        !            43: @subtitle 2017 年 3 月 3 日
1.1       takayama   44:
                     45: @author  by Y.Goto, Y.Tachibana, N.Takayama
                     46: @page
                     47: @vskip 0pt plus 1filll
                     48: Copyright @copyright{} Risa/Asir committers
                     49: 2004--2010. All rights reserved.
                     50: @end titlepage
                     51:
1.3       takayama   52: @comment --- おまじない ---
1.1       takayama   53: @synindex vr fn
1.3       takayama   54: @comment --- おまじない終り ---
1.1       takayama   55:
1.3       takayama   56: @comment --- @node は GNU info, HTML 用 ---
                     57: @comment --- @node  の引数は node-name,  next,  previous,  up ---
1.1       takayama   58: @node Top,, (dir), (dir)
                     59:
1.3       takayama   60: @comment --- @menu は GNU info, HTML 用 ---
                     61: @comment --- chapter 名を正確に並べる ---
                     62: @comment --- この文書では chapter XYZ, Chapter Index がある.
                     63: @comment ---  Chapter XYZ には section XYZについて, section XYZに関する関数がある.
1.1       takayama   64: @menu
1.3       takayama   65: * 2元分割表HGMの関数説明書について::
                     66: * 2元分割表HGMの関数::
1.4     ! takayama   67: * modular計算
1.1       takayama   68: * Index::
                     69: @end menu
                     70:
1.3       takayama   71: @comment --- chapter の開始 ---
                     72: @comment --- 親 chapter 名を正確に. 親がない場合は Top ---
                     73: @node 2元分割表HGMの関数説明書について,,, Top
                     74: @chapter 2元分割表HGMの関数説明書について
                     75:
                     76: この説明書では
                     77: HGM(holonomic gradient method) を用いた2元分割表の関数について説明する.
                     78: ChangeLog の項目は www.openxm.org の cvsweb で
                     79: ソースコードを読む時の助けになる情報が書かれている.
1.1       takayama   80:
1.3       takayama   81: 本文中で引用している文献を列挙する.
1.1       takayama   82: @itemize @bullet
                     83: @item [GM2016]
                     84: 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)
                     85: @item [T2016]
1.3       takayama   86: Y.Tachibana, 差分ホロノミック勾配法のモジュラーメソッドによる計算の高速化,
                     87: 2016, 神戸大学修士論文.
1.1       takayama   88: @item [GTT2016]
1.3       takayama   89: Y.Goto, Y.Tachibana, N.Takayama, 2元分割表に対する差分ホロノミック勾配法の実装,
                     90: 数理研講究録(掲載予定).
1.1       takayama   91: @item [TKT2015]
                     92: N.Takayama, S.Kuriki, A.Takemura,
                     93:          $A$-hypergeometric distributions and Newton polytopes.
                     94:          arxiv:1510.02269
                     95: @end itemize
                     96:
1.3       takayama   97: このマニュアルで説明する関数を用いたプログラム例は
1.1       takayama   98: gtt_ekn/test-t1.rr
1.3       takayama   99: など.
1.1       takayama  100:
1.4     ! takayama  101:
1.3       takayama  102: @node 2元分割表HGMの関数,,, Top
                    103: @chapter 2元分割表HGMの関数
1.1       takayama  104:
1.3       takayama  105: @comment --- section ``実験的関数'' の subsection xyz_abc
                    106: @comment --- subsection xyz_pqr xyz_stu がある.
1.1       takayama  107: @menu
                    108: * gtt_ekn.gmvector::
                    109: * gtt_ekn.nc::
                    110: * gtt_ekn.lognc::
                    111: * gtt_ekn.expectation::
                    112: * gtt_ekn.setup::
                    113: * gtt_ekn.upAlpha::
                    114: @end menu
                    115:
1.3       takayama  116: @node 超幾何関数E(k,n),,, 2元分割表HGMの関数
                    117: @section 超幾何関数E(k,n)
1.1       takayama  118:
                    119: @comment **********************************************************
1.3       takayama  120: @comment --- ◯◯◯◯  の説明
                    121: @comment --- 個々の関数の説明の開始 ---
                    122: @comment --- section 名を正確に ---
                    123: @node gtt_ekn.gmvector,,, 超幾何関数E(k,n)
1.1       takayama  124: @subsection @code{gtt_ekn.gmvector}
1.3       takayama  125: @comment --- 索引用キーワード
1.1       takayama  126: @findex gtt_ekn.gmvector
                    127:
                    128: @table @t
                    129: @item gtt_ekn.gmvector(@var{beta},@var{p})
1.3       takayama  130: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表に付随する超幾何関数
                    131: E(k,n) の値およびその微分の値を戻す.
1.1       takayama  132: @item gtt_ekn.ekn_cBasis_2(@var{beta},@var{p})
1.3       takayama  133: の別名である.
1.1       takayama  134: @end table
                    135:
1.3       takayama  136: @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
1.1       takayama  137: @table @var
                    138: @item return
1.3       takayama  139: ベクトル, 超幾何関数の値とその微分. 詳しくは下記.
1.1       takayama  140: @item beta
1.3       takayama  141: 行和, 列和のリスト. 成分はすべて正であること.
1.1       takayama  142: @item p
1.3       takayama  143: 二元分割表のセルの確率のリスト
1.1       takayama  144: @end table
                    145:
1.3       takayama  146: @comment --- ここで関数の詳しい説明 ---
                    147: @comment --- @itemize〜@end itemize は箇条書き ---
                    148: @comment --- @bullet は黒点付き ---
1.1       takayama  149: @itemize @bullet
                    150: @item
1.3       takayama  151: gmvector は Gauss-Manin vector の略である [GM2016].
1.1       takayama  152: @item
1.3       takayama  153: gmvector の戻り値は
                    154: [GM2016] の 6章 p.23 のベクトル Sである.
                    155: これは
                    156: [GM2016] の4章で定義されているベクトル F の定数倍であり,
                    157: その定数は
                    158: 第一成分が [GM2016] の6章で定義されている級数 S の値と等しく
                    159: なるように決められている.
1.1       takayama  160: @item
1.3       takayama  161:  r1 x r2 分割表を考える.
                    162:  m+1=r1, n+1=r2 とおく.
                    163:  正規化定数 Z は分割表 u を (m+1) × (n+1) 行列とするとき p^u/u! の和である.
                    164:  ここで和は行和列和が @var{beta} であるような u 全体でとる
1.1       takayama  165:  [TKT2015], [GM2016].
1.3       takayama  166:  S はこの多項式 Z の p を
1.1       takayama  167: @verbatim
                    168:   [[1,y11,...,y1n],
                    169:    [1,y21,...,y2n],...,
                    170:    [1,ym1, ...,ymn],
                    171:    [1,1, ..., 1]]
                    172: @end verbatim
1.3       takayama  173:  (1 が L 字型に並ぶ),
                    174: と正規化した級数である.
1.1       takayama  175: @item
1.3       takayama  176: 2x(n+1)分割表で, gmvector の戻り値を Lauricella  F_D で書くことが
                    177: 以下のようにできる
                    178: (b[2][1]-b[1][1] >= 0 の場合).
                    179: ここで b[1][1], b[1][2] は, それぞれ 1 行目の行和, 2 行目の行和,
                    180: b[2][i] は i 列目の列和である.
1.1       takayama  181: @comment ekn/Talks/2015-12-3-goto.tex
                    182: @verbatim
                    183: S=F_D(-b[1,1], [-b[2,2],...,-b[2,n+1]], b[2,1]-b[1,1]+1 ; y)/C,
                    184: @end verbatim
                    185: C=b[1,1]! b[2,2]! ... b[2][n+1]! (b[2,1]-b[1,1])!
1.3       takayama  186: とおく.
                    187: 1/C は L 字型の分割表
1.1       takayama  188: @verbatim
                    189: [[b[1,1],       0,      ..., 0       ],
                    190:  [b[2,1]-b[1,1],b[2,2], ..., b[2,n+1]]]
                    191: @end verbatim
1.3       takayama  192: に対応.
                    193: gmvector は
1.1       takayama  194: @verbatim
                    195: [S,(y11/a2) d_11 S,(y12/a3) d_12 S, ..., (y1n/a_(n+1)) d_1n S]
                    196: @end verbatim
1.3       takayama  197: である.
                    198: ここで d_ij は yij についての微分,
1.1       takayama  199: @verbatim
                    200:   [a0,     a1, ...                      ,a_(n+2)]
                    201: = [-b[1,2],-b[1,1],b[2,2], ..., b[2,n+1],b[2,1]]
                    202: @end verbatim
1.3       takayama  203: である.
1.1       takayama  204: @item
1.3       takayama  205: 周辺和 @var{beta}の時の正規化定数のセル確率 @var{p} に対する値は 多項式に退化した E(k,n) の値で表現できる. 文献 [TKT2015], [GM2016] 参照.
1.1       takayama  206: @item
1.3       takayama  207: option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう
1.1       takayama  208: [T2016].
1.3       takayama  209: 分散計算用の各種パラメータの設定は
                    210: gtt_ekn.setup で行なう.
1.1       takayama  211: @end itemize
                    212:
1.3       takayama  213: @comment --- @example〜@end example は実行例の表示 ---
                    214: 例: 次は2 x 2 分割表で行和が [5,1],  列和が [3,3], 各セルの確率が
                    215: [[1/2,1/3],[1/7,1/5]] の場合の gmvector の値である.
1.1       takayama  216: @example
                    217: [3000] load("gtt_ekn.rr");
                    218: [3001] ekn_gtt.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])
                    219: [775/27783]
                    220: [200/9261]
                    221: @end example
                    222:
1.3       takayama  223: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
                    224: 計算ができる.
                    225: 守備範囲の異なるプログラム同士の比較, debug 用参考.
1.1       takayama  226: @example
                    227: [3080] import("tk_fd.rr");
                    228: [3081] A=tk_fd.marginal2abc([4,5],[2,4,3]);
1.3       takayama  229: [-4,[-4,-3],-1]  // 2変数 FD のパラメータ. a,[b1,b2],c
1.1       takayama  230: [3082] tk_fd.fd_hessian2(A[0],A[1],A[2],[1/2,1/3]);
                    231: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    232: [4483/124416,[ 1961/15552 185/1728 ],
                    233:  [ 79/288 259/864 ]
                    234:  [ 259/864 47/288 ]]
1.3       takayama  235: // 戻値は [F=F_D, gradient(F), Hessian(F)]
1.1       takayama  236:
1.3       takayama  237: // ekn_gt での例と同じパラメータ.
1.1       takayama  238: [3543] A=tk_fd.marginal2abc([5,1],[3,3]);
                    239: [-5,[-3],-1]
                    240: [3544] tk_fd.fd_hessian2(A[0],A[1],A[2],[(1/3)*(1/7)/((1/2)*(1/5))]);
                    241: Computing Dmat(ca) for parameters B=[-3],X=[ 10/21 ]
                    242: [775/27783,[ 20/147 ],[ 17/42 ]]
                    243: @end example
                    244:
1.3       takayama  245: 参考: 一般の A 分布の正規化定数についての Hessian の計算は実験的 package ot_hessian_ahg.rr
                    246: で実装のテストがされている. (これはまだ未完成のテスト版なので出力形式等も将来的には変更される.)
1.1       takayama  247: @example
                    248: import("ot_hgm_ahg.rr");
                    249: import("ot_hessian_ahg.rr");
                    250: def  htest4() @{
                    251:   extern C11_A;
                    252:   extern C11_Beta;
                    253:   Hess=newmat(7,7);
                    254:   A =C11_A;
                    255:   Beta0= [b0,b1,b2,b3];
                    256:   BaseIdx=[4,5,6];
                    257:   X=[x0,x1,x2,x3,x4,x5,x6];
                    258:   for (I=0; I<7; I++) for (J=0; J<7; J++) @{
                    259:     Idx = [I,J];
                    260:     H=hessian_simplify(A,Beta0,X,BaseIdx,Idx);
                    261:     Hess[I][J]=H;
                    262:     printf("[I,J]=%a, Hessian_ij=%a\n",Idx,H);
                    263:   @}
                    264:   return(Hess);
                    265: @}
                    266: [2917] C11_A;
                    267: [[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]]
                    268: [2918] C11_Beta;
                    269: [166,36,290,214]
                    270: [2919] Ans=htest4$
                    271: [2920] Ans[0][0];
                    272: [[((b1-b0-1)*x4)/(x0^2),[4]],[((b1-b0-1)*x6)/(x0^2),[6]],
                    273:  [(b1^2+(-2*b0-1)*b1+b0^2+b0)/(x0^2),[]],[(x6)/(x0),[6,0]],[(x4)/(x0),[4,0]]]
                    274: @end example
                    275:
1.3       takayama  276: @comment --- 参照(リンク)を書く ---
1.1       takayama  277: @table @t
1.3       takayama  278: @item 参照
1.1       takayama  279: @ref{gtt_ekn.setup}
                    280: @ref{gtt_ekn.pfaffian_basis}
                    281: @end table
                    282:
1.3       takayama  283: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1       takayama  284: @noindent
                    285: ChangeLog
                    286: @itemize @bullet
                    287: @item
1.3       takayama  288:  この関数は
                    289: [GM2016] のアルゴリズムおよび
                    290: [T2016] による modular method を用いた高速化を実装したものである.
1.1       takayama  291: @item
1.3       takayama  292:  変更を受けたファイルは
1.1       takayama  293:  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_pfaffian_8.rr
                    294: @end itemize
                    295:
                    296:
                    297: @comment **********************************************************
1.3       takayama  298: @node gtt_ekn.nc,,, 超幾何関数E(k,n)
1.1       takayama  299: @subsection @code{gtt_ekn.nc}
1.3       takayama  300: @comment --- 索引用キーワード
1.1       takayama  301: @findex gtt_ekn.nc
                    302:
                    303: @table @t
                    304: @item gtt_ekn.nc(@var{beta},@var{p})
1.3       takayama  305: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z
                    306: およびその微分の値を戻す.
1.1       takayama  307: @end table
                    308:
1.3       takayama  309: @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
1.1       takayama  310: @table @var
                    311: @item return
1.3       takayama  312: ベクトル [Z,[[d_11 Z, d_12 Z, ...], ..., [d_m1 Z, d_m2 Z, ...., d_mn Z]]]
1.1       takayama  313: @item beta
1.3       takayama  314: 行和, 列和のリスト. 成分はすべて正であること.
1.1       takayama  315: @item p
1.3       takayama  316: 二元分割表のセルの確率のリスト
1.1       takayama  317: @end table
                    318:
1.3       takayama  319: @comment --- ここで関数の詳しい説明 ---
                    320: @comment --- @itemize〜@end itemize は箇条書き ---
                    321: @comment --- @bullet は黒点付き ---
1.1       takayama  322: @itemize @bullet
                    323: @item
1.3       takayama  324:  r1 x r2 分割表を考える.
                    325:  m=r1, n=r2 とおく.
                    326:  正規化定数 Z は分割表 u を m × n 行列とするとき p^u/u! の和である.
                    327:  ここで和は行和列和が @var{beta} であるような u 全体でとる
1.1       takayama  328:  [TKT2015], [GM2016].
1.3       takayama  329:  p^u は p_ij^u_ij の積, u! は u_ij! の積である.
                    330:  d_ij Z で Z の変数 p_ij についての偏微分を表す.
1.1       takayama  331: @item
1.3       takayama  332: nc は gmvector の値を元に, [GM2016] の Prop
                    333:  7.1 に基づいて Z の値を計算する.
1.1       takayama  334: @item
1.3       takayama  335: option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
                    336: 分散計算用の各種パラメータの設定は
                    337: gtt_ekn.setup で行なう.
1.1       takayama  338: @end itemize
                    339:
1.3       takayama  340: @comment --- @example〜@end example は実行例の表示 ---
                    341: 例: 2x3 分割表での Z とその微分の計算.
1.1       takayama  342: @example
                    343: [2237] gtt_ekn.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
                    344: [4483/124416,[ 353/7776 1961/15552 185/1728 ]
                    345: [ 553/20736 1261/15552 1001/13824 ]]
                    346: @end example
                    347:
1.3       takayama  348: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
                    349: 計算ができる.
1.1       takayama  350: @example
                    351: [3076] import("tk_fd.rr");
                    352: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
                    353: [-4,[-4,-3],-1]
                    354: [3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
                    355: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
                    356: [ 1 1 1 ]
                    357: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    358: [4483/124416,[[353/7776,1961/15552,185/1728],
                    359:               [553/20736,1261/15552,1001/13824]]]
1.3       takayama  360: // 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],
                    361: //             [d_21 Z, d_22 Z, d_23 Z]]] の値.
                    362: //           ここで d_ij は i,j 成分についての微分を表す.
1.1       takayama  363: @end example
                    364:
1.3       takayama  365: @comment --- 参照(リンク)を書く ---
1.1       takayama  366: @table @t
1.3       takayama  367: @item 参照
1.1       takayama  368: @ref{gtt_ekn.setup}
                    369: @ref{gtt_ekn.lognc}
                    370: @end table
                    371:
1.3       takayama  372: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1       takayama  373: @noindent
                    374: ChangeLog
                    375: @itemize @bullet
                    376: @item
1.3       takayama  377:  変更を受けたファイルは
1.1       takayama  378:  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_eval.rr
                    379: @end itemize
                    380:
                    381:
                    382: @comment **********************************************************
1.3       takayama  383: @node gtt_ekn.lognc,,, 超幾何関数E(k,n)
1.1       takayama  384: @subsection @code{gtt_ekn.lognc}
1.3       takayama  385: @comment --- 索引用キーワード
1.1       takayama  386: @findex gtt_ekn.lognc
                    387:
                    388: @table @t
                    389: @item gtt_ekn.lognc(@var{beta},@var{p})
1.3       takayama  390: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z
                    391: の log の近似値およびその微分の近似値を戻す.
1.1       takayama  392: @end table
                    393:
1.3       takayama  394: @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
1.1       takayama  395: @table @var
                    396: @item return
1.3       takayama  397: ベクトル [log(Z), [[d_11 log(Z), d_12 log(Z), ...], [d_21 log(Z),...], ... ]
1.1       takayama  398: @item beta
1.3       takayama  399: 行和, 列和のリスト. 成分はすべて正であること.
1.1       takayama  400: @item p
1.3       takayama  401: 二元分割表のセルの確率のリスト
1.1       takayama  402: @end table
                    403:
1.3       takayama  404: @comment --- ここで関数の詳しい説明 ---
                    405: @comment --- @itemize〜@end itemize は箇条書き ---
                    406: @comment --- @bullet は黒点付き ---
1.1       takayama  407: @itemize @bullet
                    408: @item
1.3       takayama  409: 条件付き最尤推定に利用する [TKT2015].
                    410: @item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
                    411: 分散計算用の各種パラメータの設定は
                    412: gtt_ekn.setup で行なう.
1.1       takayama  413: @end itemize
                    414:
1.3       takayama  415: @comment --- @example〜@end example は実行例の表示 ---
                    416: 例: 2 × 3 分割表での例. 第一成分のみ近似値.
1.1       takayama  417: @example
                    418: [2238] gtt_ekn.lognc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
                    419: [-3.32333832422461674630,[ 5648/4483 15688/4483 13320/4483 ]
                    420: [ 3318/4483 10088/4483 9009/4483 ]]
                    421: @end example
                    422:
1.3       takayama  423: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
                    424: 計算ができる.
1.1       takayama  425: @example
                    426: [3076] import("tk_fd.rr");
                    427: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
                    428: [-4,[-4,-3],-1]
                    429: [3078] tk_fd.log_ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
                    430: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
                    431: [ 1 1 1 ]
                    432: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    433: [-3.32333832422461674639485797719209322217260539267246045320,
                    434:  [[1.2598706, 3.499442, 2.971224],
                    435:   [0.7401293, 2.250278, 2.009591]]]
1.3       takayama  436: // 戻値は [log(Z),
1.1       takayama  437: //          [[d_11 log(Z), d_12 log(Z), d_13 log(Z)],
                    438: //           [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]]
1.3       takayama  439: // の近似値.
1.1       takayama  440: @end example
                    441:
1.3       takayama  442: @comment --- 参照(リンク)を書く ---
1.1       takayama  443: @table @t
1.3       takayama  444: @item 参照
1.1       takayama  445: @ref{gtt_ekn.setup}
                    446: @ref{gtt_ekn.nc}
                    447: @end table
                    448:
1.3       takayama  449: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1       takayama  450: @noindent
                    451: ChangeLog
                    452: @itemize @bullet
                    453: @item
1.3       takayama  454:  変更を受けたファイルは
1.1       takayama  455:  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.
                    456: @end itemize
                    457:
                    458: @comment **********************************************************
1.3       takayama  459: @node gtt_ekn.expectation,,, 超幾何関数E(k,n)
1.1       takayama  460: @subsection @code{gtt_ekn.expectation}
1.3       takayama  461: @comment --- 索引用キーワード
1.1       takayama  462: @findex gtt_ekn.expectation
                    463:
                    464: @table @t
                    465: @item gtt_ekn.expectation(@var{beta},@var{p})
1.3       takayama  466: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の期待値を計算する.
1.1       takayama  467: @end table
                    468:
1.3       takayama  469: @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
1.1       takayama  470: @table @var
                    471: @item return
1.3       takayama  472: 二元分割表の各セルの期待値のリスト.
1.1       takayama  473: @item beta
1.3       takayama  474: 行和, 列和のリスト. 成分はすべて正であること.
1.1       takayama  475: @item p
1.3       takayama  476: 二元分割表のセルの確率のリスト
1.1       takayama  477: @end table
                    478:
1.3       takayama  479: @comment --- ここで関数の詳しい説明 ---
                    480: @comment --- @itemize〜@end itemize は箇条書き ---
                    481: @comment --- @bullet は黒点付き ---
1.1       takayama  482: @itemize @bullet
                    483: @item
1.3       takayama  484: [GM2016] の Algorithm 7.8 の実装.
                    485: @item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
                    486: 分散計算用の各種パラメータの設定は
                    487: gtt_ekn.setup で行なう.
                    488: @item option index を与えると, 指定された成分の期待値のみ計算する.
                    489: たとえば 2 x 2 分割表で index=[[0,0],[1,1]] と指定すると, 1 のある成分の期待値のみ計算する.
1.1       takayama  490: @end itemize
                    491:
1.3       takayama  492: @comment --- @example〜@end example は実行例の表示 ---
1.1       takayama  493:
1.3       takayama  494: 2×2, 3×3 の分割表の期待値計算例.
1.1       takayama  495: @example
                    496: [2235] gtt_ekn.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);
                    497: [ 2/3 1/3 ]
                    498: [ 4/3 8/3 ]
                    499: [2236] gtt_ekn.expectation([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
                    500: [ 5648/4483 7844/4483 4440/4483 ]
                    501: [ 3318/4483 10088/4483 9009/4483 ]
                    502:
                    503: [2442] gtt_ekn.expectation([[4,14,9],[11,6,10]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]]);
                    504: [ 207017568232262040/147000422096729819 163140751505489940/147000422096729819
                    505:                                         217843368649167296/147000422096729819 ]
                    506: [ 1185482401011137878/147000422096729819 358095302885438604/147000422096729819
                    507:                                          514428205457640984/147000422096729819 ]
                    508: [ 224504673820628091/147000422096729819 360766478189450370/147000422096729819
                    509:                                         737732646860489910/147000422096729819 ]
                    510: @end example
                    511:
1.3       takayama  512: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
                    513: 計算ができる.
1.1       takayama  514: @example
                    515: [3076] import("tk_fd.rr");
                    516: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
                    517: [-4,[-4,-3],-1]
                    518: [3078] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
                    519: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
                    520: [ 1 1 1 ]
                    521: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    522: [[5648/4483,7844/4483,4440/4483],
                    523:  [3318/4483,10088/4483,9009/4483]]
1.3       takayama  524: // 各セルの期待値.
1.1       takayama  525: @end example
                    526:
1.3       takayama  527: 参考: 一般の A 分布の計算は ot_hgm_ahg.rr. まだ実験的なため, module 化されていない.
                    528: ot_hgm_ahg.rr についての参考文献:
1.1       takayama  529: K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Systems II --- Holonomic Gradient Method, arxiv:1505.02947
                    530: @example
                    531: [3237] import("ot_hgm_ahg.rr");
1.3       takayama  532: // 2 x 2 分割表.
1.1       takayama  533: [3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],
                    534:         [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
                    535: oohg_native=0, oohg_curl=1
                    536: [1376777025/625400597,1750225960/625400597,
                    537:  2375626557/625400597,3252978816/625400597]
1.3       takayama  538: // 2 x 2 分割表の期待値.
1.1       takayama  539:
1.3       takayama  540: // 2 x 3 分割表.
1.1       takayama  541: [3238] hgm_ahg_expected_values_contiguity(
                    542:  [[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
                    543:  [5,2,4,3],[1,1/2,1/3,1,1,1],[x1,x2,x3,x4,x5,x6]|geometric=1);
                    544: [5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
1.3       takayama  545: // 2 x 3 分割表の期待値. 上と同じ問題.
1.1       takayama  546: @end example
                    547:
1.3       takayama  548: 3 x 3 分割表. 構造的0が一つ.
1.1       takayama  549: @example
                    550: /*
1.3       takayama  551:   dojo, p.221 のデータ.  成績3以下の生徒は集めてひとつに.
1.1       takayama  552:   2 1 1
                    553:   8 3 3
                    554:   0 2 6
                    555:
                    556:   row sum: 4,14,8
                    557:   column sum: 10,6,10
1.3       takayama  558:   0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
1.1       takayama  559: */
                    560:
                    561: A=[[0,0,0,1,1,1, 0,0],
                    562:    [0,0,0,0,0,0, 1,1],
                    563:    [1,0,0,1,0,0, 0,0],
                    564:    [0,1,0,0,1,0, 1,0],
                    565:    [0,0,1,0,0,1, 0,1]];
                    566: B=[14,8,10,6,10];
                    567: hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1],
1.3       takayama  568:                 [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);
1.1       takayama  569:
1.3       takayama  570: // 答.
1.1       takayama  571: [14449864949304/9556267369631,
                    572:  10262588586540/9556267369631, 13512615942680/9556267369631,
                    573:  81112808747006/9556267369631,
                    574:  21816297744346/9556267369631, 30858636683482/9556267369631,
                    575:
                    576:  25258717886900/9556267369631,51191421070148/9556267369631]
                    577: @end example
                    578:
1.3       takayama  579: 3 x 3 分割表.
1.1       takayama  580: @example
                    581: /*
1.3       takayama  582:  上のデータで 0 を 1 に変更.
1.1       takayama  583:   2 1 1
                    584:   8 3 3
                    585:   1 2 6
                    586:
                    587:   row sum: 4,14,9
                    588:   column sum: 11,6,10
                    589: */
                    590: A=[[0,0,0,1,1,1,0,0,0],
                    591:    [0,0,0,0,0,0,1,1,1],
                    592:    [1,0,0,1,0,0,1,0,0],
                    593:    [0,1,0,0,1,0,0,1,0],
                    594:    [0,0,1,0,0,1,0,0,1]];
                    595: B=[14,9,11,6,10];
                    596: hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1,1],
                    597:                               [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);
                    598:
1.3       takayama  599: // 期待値, 答.   x9 を指定していないので, 9番目の期待値は出力してない.
1.1       takayama  600: [207017568232262040/147000422096729819,
                    601:  163140751505489940/147000422096729819,217843368649167296/147000422096729819,
                    602:  1185482401011137878/147000422096729819,
                    603:  358095302885438604/147000422096729819,514428205457640984/147000422096729819,
                    604:  224504673820628091/147000422096729819,360766478189450370/147000422096729819]
                    605:
1.3       takayama  606: // Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
                    607: // まだ書いてない.
1.1       takayama  608: @end example
                    609:
                    610:
                    611:
1.3       takayama  612: @comment --- 参照(リンク)を書く ---
1.1       takayama  613: @table @t
1.3       takayama  614: @item 参照
1.1       takayama  615: @ref{gtt_ekn.setup}
                    616: @ref{gtt_ekn.nc}
                    617: @end table
                    618:
1.3       takayama  619: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1       takayama  620: @noindent
                    621: ChangeLog
                    622: @itemize @bullet
                    623: @item
1.3       takayama  624:  変更を受けたファイルは
1.1       takayama  625:  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.
                    626: @end itemize
                    627:
                    628:
                    629: @comment **********************************************************
1.3       takayama  630: @comment --- ◯◯◯◯  の説明
                    631: @comment --- 個々の関数の説明の開始 ---
                    632: @comment --- section 名を正確に ---
                    633: @node gtt_ekn.setup,,, 超幾何関数E(k,n)
1.1       takayama  634: @subsection @code{gtt_ekn.setup}
1.3       takayama  635: @comment --- 索引用キーワード
1.1       takayama  636: @findex gtt_ekn.setup
                    637:
                    638: @table @t
                    639: @item gtt_ekn.setup()
1.3       takayama  640: :: 分散計算用の環境設定をおこなう. 現在の環境を報告する.
1.1       takayama  641: @end table
                    642:
1.3       takayama  643: @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
1.1       takayama  644: @table @var
                    645: @item return
                    646:
                    647: @end table
                    648:
1.3       takayama  649: @comment --- ここで関数の詳しい説明 ---
                    650: @comment --- @itemize〜@end itemize は箇条書き ---
                    651: @comment --- @bullet は黒点付き ---
                    652: @itemize @bullet
                    653: @item 使用するプロセスと素数の個数, 最小の素数を表示する. 準備されていない場合はその旨を表示.
                    654: @item このパッケージでの分散計算は複数のcpuを搭載した計算機で実行されることを想定している.
                    655: @item option nps (または number_of_processes)を与えると指定した数だけプロセスを用意する.
                    656: @item option nprm (または number_of_primes)を与えるとnprmが文字列の場合指定された素数リストのファイルを読み込む. nprmが自然数の場合さらにoption minp (minp =MINimum Prime)を与えるとminpより大きな素数をnprm個生成する. その際option fgp (または file_of_generated_primes)を与えると生成した素数リストをファイル名をfgpとして保存する.
                    657: @item 上記のoption を指定しなかった場合次のデフォルト値が用いられる. nps=1. nprm=10. fgp=0.
                    658: @item option report=1を与えると現在の環境の報告のみを行う. setup(|report=1)の別名としてreport関数を使用することもできる.
1.4     ! takayama  659: @item option subprogs=[file1,file2,...] により分散計算の子供プロセスにロードすべきファイル file1, file2, ... を指定する. default は subprogs=["gtt_ekn/childprocess.rr"] である.
1.1       takayama  660: @end itemize
                    661:
1.3       takayama  662: @comment --- @example〜@end example は実行例の表示 ---
                    663: 例: 素数のリストを生成してファイル p.txt へ書き出す.
1.1       takayama  664: @example
                    665: gtt_ekn.setup(|nps=2,nprm=20,minp=10^10,fgp="p.txt")$
                    666: @end example
                    667:
                    668:
1.3       takayama  669: @comment --- 参照(リンク)を書く ---
1.1       takayama  670: @table @t
1.3       takayama  671: @item 参照
1.1       takayama  672: @ref{gtt_ekn.nc}
                    673: @ref{gtt_ekn.gmvector}
                    674: @end table
                    675:
1.3       takayama  676: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1       takayama  677: @noindent
                    678: ChangeLog
                    679: @itemize @bullet
                    680: @item
1.3       takayama  681:  変更を受けたファイルは
1.1       takayama  682:  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1,
                    683:  gtt_ekn/g_mat_fac.rr
                    684:
                    685: @end itemize
                    686:
                    687: @comment **********************************************************
1.3       takayama  688: @comment --- ◯◯◯◯  の説明
                    689: @comment --- 個々の関数の説明の開始 ---
                    690: @comment --- section 名を正確に ---
                    691: @node gtt_ekn.upAlpha,,, 超幾何関数E(k,n)
1.1       takayama  692: @subsection @code{gtt_ekn.upAlpha}
1.3       takayama  693: @comment --- 索引用キーワード
1.1       takayama  694: @findex gtt_ekn.upAlpha
                    695:
                    696: @table @t
                    697: @item gtt_ekn.upAlpha(@var{i},@var{k},@var{n})
                    698: ::
                    699: @end table
                    700:
1.3       takayama  701: @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
1.1       takayama  702: @table @var
1.3       takayama  703: @item i  a_i を a_i+1 と変化させる contiguity relation.
                    704: @item k  E(k+1,n+k+2)型の超幾何関数の k. 分割表では (k+1)×(n+1).
                    705: @item n  E(k+1,n+k+2)型の超幾何関数の n. 分割表では (k+1)×(n+1).
                    706: @item return  contiguity relation の pfaffian_basis についての行列表現を戻す. [GM2016] の Cor 6.3.
1.1       takayama  707: @end table
                    708:
1.3       takayama  709: @comment --- ここで関数の詳しい説明 ---
                    710: @comment --- @itemize〜@end itemize は箇条書き ---
                    711: @comment --- @bullet は黒点付き ---
1.1       takayama  712: @itemize @bullet
                    713: @item
1.3       takayama  714:  upAlpha は [GM2016] の Cor 6.3 の行列 U_i を戻す.
                    715: @item 関連する各関数の簡潔な説明と例も加える.
                    716: @item a_i を a_i-1 と変化させたい場合は関数 downAlpha を用いる.
                    717: @item a_i と分割表の周辺和を見るには, 関数 marginaltoAlpha([行和,列和]) を用いる.
1.1       takayama  718: @item
1.3       takayama  719:    pfaffian_basis は [GM2016] の4章のベクトル F に対応する偏微分を戻す.
1.1       takayama  720: @end itemize
                    721:
1.3       takayama  722: @comment --- @example〜@end example は実行例の表示 ---
                    723: 例: 以下の例は 2×2分割表(E(2,4)), 2×3分割表(E(2,5))の場合である.
                    724: [2225] までは出力を略している.
1.1       takayama  725: @example
                    726: [2221] gtt_ekn.marginaltoAlpha([[1,4],[2,3]]);
                    727: [[a_0,-4],[a_1,-1],[a_2,3],[a_3,2]]
1.3       takayama  728: [2222] gtt_ekn.upAlpha(1,1,1);  // E(2,4) の a_1 方向の
                    729:                                 //     contiguity を表現する行列
                    730: [2223] gtt_ekn.upAlpha(2,1,1);  // E(2,4) の a_2 方向
                    731: [2224] gtt_ekn.upAlpha(3,1,1);  // E(2,4) の a_3 方向
1.1       takayama  732: [2225] function f(x_1_1);
                    733: [2232] gtt_ekn.pfaffian_basis(f(x_1_1),1,1);
                    734: [ f(x_1_1) ]
                    735: [ (f{1}(x_1_1)*x_1_1)/(a_2) ]
                    736: [2233] function f(x_1_1,x_1_2);
                    737: f() redefined.
1.3       takayama  738: [2234] gtt_ekn.pfaffian_basis(f(x_1_1,x_1_2),1,2); // E(2,5), 2*3 分割表
1.1       takayama  739: [ f(x_1_1,x_1_2) ]
                    740: [ (f{1,0}(x_1_1,x_1_2)*x_1_1)/(a_2) ]
                    741: [ (f{0,1}(x_1_1,x_1_2)*x_1_2)/(a_3) ]
                    742: @end example
                    743:
                    744:
1.3       takayama  745: @comment --- 参照(リンク)を書く ---
1.1       takayama  746: @table @t
1.3       takayama  747: @item 参照
1.1       takayama  748: @ref{gtt_ekn.nc}
                    749: @ref{gtt_ekn.gmvector}
                    750: @end table
                    751:
1.3       takayama  752: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1       takayama  753: @noindent
                    754: ChangeLog
                    755: @itemize @bullet
                    756: @item
1.3       takayama  757:  この関数は [GM2016]
                    758: で与えられたアルゴリズムに従い contiguity relation を導出する.
1.1       takayama  759: @item
1.3       takayama  760:  変更を受けたファイルは
1.1       takayama  761:  OpenXM/src/asir-contrib/packages/src/gtt_ekn/ekn_pfaffian_8.rr 1.1.
                    762: @end itemize
                    763:
                    764:
1.4     ! takayama  765: @node modular計算,,, 2元分割表HGMの関数
        !           766: @chapter modular計算
        !           767:
        !           768: @menu
        !           769: * gtt_ekn.chinese_itor::
        !           770: @end menu
        !           771:
        !           772: @node 中国剰余定理とitor,,, modular計算
        !           773: @section 中国剰余定理とitor
        !           774:
        !           775: @comment **********************************************************
        !           776: @comment --- ◯◯◯◯  の説明
        !           777: @comment --- 個々の関数の説明の開始 ---
        !           778: @comment --- section 名を正確に ---
        !           779: @node gtt_ekn.chinese_itor,,,
        !           780: @subsection @code{gtt_ekn.chinese_itor}
        !           781: @comment --- 索引用キーワード
        !           782: @findex gtt_ekn.chinese_itor 中国剰余定理とitor
        !           783:
        !           784: @table @t
        !           785: @item gtt_ekn.chinese_itor(@var{data},@var{idlist})
        !           786: :: mod p で計算した結果(ベクトル)から chinese remainder theorem, itor(integer to rational) で有理数ベクトルを得る.
        !           787: @end table
        !           788:
        !           789: @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
        !           790: @table @var
        !           791: @item return  [val, n]  ここで val は答え. また, n = n1*n2*...
        !           792: @item data   [[val1,n1],[val2,n2], ...], ここで val mod n1 = val1, val mod n2 = val2,...
        !           793: @item idlist  chinese, itor を実行するサーバIDのリスト.
        !           794: @end table
        !           795:
        !           796: @comment --- ここで関数の詳しい説明 ---
        !           797: @comment --- @itemize〜@end itemize は箇条書き ---
        !           798: @comment --- @bullet は黒点付き ---
        !           799: @itemize @bullet
        !           800: @item 中国剰余定理を用いて val0 mod n1 = val1, val0 mod n2 = val2, ... となる val0 を求める. val に algorithm itor を適用する.
        !           801: @item sqrt(n) より val0 が大きい時は itor が適用されて val0 が有理数 val=a/b に変換される. つまり b*x =1 mod n となる逆数 x を考えて, x*a % n = val0 となる数 val を戻す. 見つからないときは failure を戻す.
        !           802: @end itemize
        !           803:
        !           804: @comment --- @example〜@end example は実行例の表示 ---
        !           805: 例: [3!, 5^3*3!]=[6,750] が戻り値.
        !           806: 6 mod 109 =6, 750 mod 109=96 が最初の引数の [[6,96],109]. 以下同様.
        !           807: @example
        !           808: gtt_ekn.setup(|nps=2,nprm=3,minp=101,fgp="p_small.txt");
        !           809: SS=gtt_ekn.get_svalue();
        !           810: SS[0];
        !           811:   [103,107,109]   // list of primes
        !           812: SS[1];
        !           813:   [0,2]           // list of server ID's
        !           814: gtt_ekn.chinese_itor([[[ 6,96 ],109],[[ 6,29 ],103],[[ 6,1 ],107]],SS[1]);
        !           815:   [[ 6 750 ],1201289]
        !           816:
        !           817: // 引数はスカラーでもよい.
        !           818: gtt_ekn.chinese_itor([[96,109],[29,103]],SS[1]);
        !           819:   [[ 750 ],11227]
        !           820: @end example
        !           821:
        !           822:
        !           823: @comment --- @example〜@end example は実行例の表示 ---
        !           824: 例: gtt_ekn/childprocess.rr (server で実行される) の関数 chinese (chinese remainder theorem) と euclid.
        !           825: @example
        !           826: load("gtt_ekn/childprocess.rr");
        !           827: chinese([newvect(2,[6,29]),103],[newvect(2,[6,750]),107*109]);
        !           828:   // mod 103 で [6,29], mod (107*109) で [6,750] となる数を mod 103*(107*109)
        !           829:   // で求めると,
        !           830:   [[ 6 750 ],1201289]
        !           831: euclid(3,103);  // mod 103 での 3 の逆数. つまり 1/3
        !           832:   -34
        !           833: 3*(-34) % 103; // 確かに逆数.
        !           834:    1
        !           835: @end example
        !           836:
        !           837: @comment --- @example〜@end example は実行例の表示 ---
        !           838: 例: gtt_ekn/childprocess.rr (server で実行される) の関数 itor (integer to rational) の例.
        !           839: itor(Y,Q,Q2,Idx) では Y < Q2 なら Y がそのまま戻る.  Idx は 内部用の index で好きな数でよい. 戻り値の第2成分となる.
        !           840: @example
        !           841: load("gtt_ekn/childprocess.rr");
        !           842: for (I=1;I<11; I++) print([I,itor(I,11,3,0)]);
        !           843: [1,[1,0]]
        !           844: [2,[2,0]]
        !           845: [3,[-2/3,0]] //euclid(3,11); ->4,  4*(-2)%11 -> 3 なので確かに -2/3 は元の数の候補
        !           846: [4,[failure,0]]
        !           847: [5,[-1/2,0]]
        !           848: [6,[1/2,0]]
        !           849: [7,[-1/3,0]]
        !           850: [8,[failure,0]]
        !           851: [9,[-2,0]]
        !           852: [10,[-1,0]]
        !           853: @end example
        !           854:
        !           855:
        !           856: @comment --- 参照(リンク)を書く ---
        !           857: @table @t
        !           858: @item 参照
        !           859: @ref{gtt_ekn.setup}
        !           860: @end table
        !           861:
        !           862: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
        !           863: @noindent
        !           864: ChangeLog
        !           865: @itemize @bullet
        !           866: @item
        !           867:  関連ファイルは
        !           868:  gtt_ekn/g_mat_fac.rr
        !           869:  gtt_ekn/childprocess.rr
        !           870: @end itemize
        !           871:
        !           872:
1.1       takayama  873:
1.3       takayama  874: @comment --- おまじない ---
1.1       takayama  875: @node Index,,, Top
                    876: @unnumbered Index
                    877: @printindex fn
                    878: @printindex cp
                    879: @iftex
                    880: @vfill @eject
                    881: @end iftex
                    882: @summarycontents
                    883: @contents
                    884: @bye
1.3       takayama  885: @comment --- おまじない終り ---
1.1       takayama  886:
                    887:
1.3       takayama  888: // 2 x m 分割表において似た機能を有する関数の利用例を参考までに記載する;
                    889: // 正規化定数とその微分関連.
                    890: // その1.
1.1       takayama  891: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
                    892: [-4,[-4,-3],-1]
                    893: [3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
                    894: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
                    895: [ 1 1 1 ]
                    896: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    897: [4483/124416,[[353/7776,1961/15552,185/1728],[553/20736,1261/15552,1001/13824]]]
1.3       takayama  898: // 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],[d_21 Z, d_22 Z, d_23 Z]]] の値.
1.1       takayama  899:
1.3       takayama  900: // その2.
1.1       takayama  901: [3079] tk_fd.log_ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
                    902: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
                    903: [ 1 1 1 ]
                    904: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    905: [-3.32333832422461674639485797719209322217260539267246045320,
                    906:  [[1.25987062235110417131385233102832924380994869507026544724,3.49944233772027660049074280615659156814633058219942003122,2.97122462636627258532232879768012491635065804149007361142],
                    907:   [0.740129377648895828686147668971670756190051304929734552754,2.25027883113986169975462859692170421592683470890028998438,2.00959179121124247155922373410662502788311398616997546285]]]
1.3       takayama  908: // 戻値は [log(Z),
1.1       takayama  909: //          [[d_11 log(Z), d_12 log(Z), d_13 log(Z)],
                    910: //           [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]]
1.3       takayama  911: // の近似値.
1.1       takayama  912:
1.3       takayama  913: // その3.
1.1       takayama  914: [3082] fd_hessian2(A[0],A[1],A[2],[1/2,1/3]);
                    915: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    916: [4483/124416,[ 1961/15552 185/1728 ],
                    917:  [ 79/288 259/864 ]
                    918:  [ 259/864 47/288 ]]
1.3       takayama  919: // 戻値は [F=F_D, gradient(F), Hessian(F)]
1.1       takayama  920:
1.3       takayama  921: // 参考.
                    922: // ygahvec で巾関数分の調整. 独立した関数はないようだ.
1.1       takayama  923:
                    924: //-----------------------------------------------------------------------
1.3       takayama  925: // 2 x m 分割表において似た機能を有する関数の利用例を参考までに記載する;
                    926: // 期待値関連.
1.1       takayama  927: [3079] A=tk_fd.marginal2abc([4,5],[2,4,3]);
                    928: [-4,[-4,-3],-1]
                    929: [3080] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
                    930: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
                    931: [ 1 1 1 ]
                    932: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
                    933: [[5648/4483,7844/4483,4440/4483],
                    934:  [3318/4483,10088/4483,9009/4483]]
1.3       takayama  935: // 各セルの期待値.
1.1       takayama  936:
                    937: //-----------------------------------------------------------------------
1.3       takayama  938: // ot_hgm_ahg.rr の例.  実験的なため module 化されていない.
1.1       takayama  939: [3237] import("ot_hgm_ahg.rr");
1.3       takayama  940: // 2 x 2 分割表.
1.1       takayama  941: [3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],
                    942:         [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
                    943: oohg_native=0, oohg_curl=1
                    944: [1376777025/625400597,1750225960/625400597,2375626557/625400597,3252978816/625400597]
1.3       takayama  945: // 2 x 2 分割表の期待値.
1.1       takayama  946:
1.3       takayama  947: // 2 x 3 分割表.
1.1       takayama  948: [3238] hgm_ahg_expected_values_contiguity(
                    949:  [[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
                    950:  [5,2,4,3],[1,1/2,1/3,1,1,1],[x1,x2,x3,x4,x5,x6]|geometric=1);
                    951: [5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
1.3       takayama  952: // 2 x 3 分割表の期待値. 上と同じ問題.
1.1       takayama  953:
                    954: /*
1.3       takayama  955:   dojo, p.221.  成績3以下の生徒は集めてひとつに.
1.1       takayama  956:   2 1 1
                    957:   8 3 3
                    958:   0 2 6
                    959:
                    960:   row sum: 4,14,8
                    961:   column sum: 10,6,10
1.3       takayama  962:   0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
1.1       takayama  963: */
1.3       takayama  964: // 3 x 3 分割表. 構造的0が一つ.
1.1       takayama  965:
                    966: A=[[0,0,0,1,1,1, 0,0],
                    967:    [0,0,0,0,0,0, 1,1],
                    968:    [1,0,0,1,0,0, 0,0],
                    969:    [0,1,0,0,1,0, 1,0],
                    970:    [0,0,1,0,0,1, 0,1]];
                    971: B=[14,8,10,6,10];
                    972: 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);
                    973:
1.3       takayama  974: // 答.
1.1       takayama  975: [14449864949304/9556267369631,10262588586540/9556267369631,13512615942680/9556267369631,
                    976:  81112808747006/9556267369631,21816297744346/9556267369631,30858636683482/9556267369631,
                    977:                               25258717886900/9556267369631,51191421070148/9556267369631]
                    978:
                    979:
                    980: /*
1.3       takayama  981:  上のデータで 0 を 1 に変更.
1.1       takayama  982:   2 1 1
                    983:   8 3 3
                    984:   1 2 6
                    985:
                    986:   row sum: 4,14,9
                    987:   column sum: 11,6,10
                    988: */
1.3       takayama  989: // 3 x 3 分割表.
1.1       takayama  990: A=[[0,0,0,1,1,1,0,0,0],
                    991:    [0,0,0,0,0,0,1,1,1],
                    992:    [1,0,0,1,0,0,1,0,0],
                    993:    [0,1,0,0,1,0,0,1,0],
                    994:    [0,0,1,0,0,1,0,0,1]];
                    995: B=[14,9,11,6,10];
                    996: 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);
                    997:
1.3       takayama  998: // 期待値, 答.
1.1       takayama  999: [207017568232262040/147000422096729819,163140751505489940/147000422096729819,217843368649167296/147000422096729819,
                   1000:  1185482401011137878/147000422096729819,358095302885438604/147000422096729819,514428205457640984/147000422096729819,
                   1001:  224504673820628091/147000422096729819,360766478189450370/147000422096729819]
                   1002:
1.3       takayama 1003: // Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
                   1004: // まだ書いてない.
1.1       takayama 1005:
                   1006:
1.3       takayama 1007: 4. x_ij は [GM2016] の1章で,
                   1008:  たとえば 3x3 の時 [[1,1,1],[x_11,x_12,1],[x_21,x_22,1]]
                   1009: となっているが, [GM2016] の Prop 7.1 の対応では,
                   1010:    p = [[1,x_11,x_12],[1,x_21,x_22],[1,1,1]] となっているので注意.

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>