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

Diff for /OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi between version 1.1 and 1.15

version 1.1, 2016/03/21 00:16:10 version 1.15, 2019/03/20 02:08:55
Line 1 
Line 1 
 %% $OpenXM$  %% $OpenXM: OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi,v 1.14 2019/03/19 07:36:21 takayama Exp $
 %% ptex gtt_ekn.texi   (.texi $B$^$G$D$1$k(B. platex $B$G$J$/(B ptex)  %% xetex gtt_ekn-ja.texi   (.texi までつける. )
 %% $B0J2<%3%a%s%H$O(B @comment $B$G;O$a$k(B.  \input texinfo $B0J9_$OIaDL$N(B tex $BL?Na$O;H$($J$$(B.  %% 以下コメントは @comment で始める.  \input texinfo 以降は普通の tex 命令は使えない.
 \input texinfo  \input texinfo-ja
 @iftex  @iftex
 @catcode`@#=6  @catcode`@#=6
 @def@fref#1{@xrefX[#1,,@code{#1},,,]}  @def@fref#1{@xrefX[#1,,@code{#1},,,]}
 @def@b#1{{@bf@gt #1}}  @def@b#1{{@bf #1}}
 @catcode`@#=@other  @catcode`@#=@other
 @end iftex  @end iftex
 @overfullrule=0pt  @overfullrule=0pt
   @documentlanguage ja
 @c -*-texinfo-*-  @c -*-texinfo-*-
 @comment %**start of header  @comment %**start of header
 @comment --- $B$*$^$8$J$$=*$j(B ---  @comment --- おまじない終り ---
   
 @comment --- GNU info $B%U%!%$%k$NL>A0(B ---  @comment --- GNU info ファイルの名前 ---
 @setfilename xyzman  @setfilename xyzman
   
 @comment --- $B%?%$%H%k(B ---  @comment --- タイトル ---
 @settitle 2$B85J,3dI=(BHGM  @settitle 2元分割表HGM
   
 @comment %**end of header  @comment %**end of header
 @comment %@setchapternewpage odd  @comment %@setchapternewpage odd
   
 @comment --- $B$*$^$8$J$$(B ---  @comment --- おまじない ---
 @ifinfo  @ifinfo
 @macro fref{name}  @macro fref{name}
 @ref{\name\,,@code{\name\}}  @ref{\name\,,@code{\name\}}
Line 34 
Line 35 
 @end iftex  @end iftex
   
 @titlepage  @titlepage
 @comment --- $B$*$^$8$J$$=*$j(B ---  @comment --- おまじない終り ---
   
 @comment --- $B%?%$%H%k(B, $B%P!<%8%g%s(B, $BCx<TL>(B, $BCx:n8"I=<((B ---  @comment --- タイトル, バージョン, 著者名, 著作権表示 ---
 @title 2$B85J,3dI=(BHGM$B4X?t(B  @title 2元分割表HGM関数
 @subtitle Risa/Asir 2$B85J,3dI=(BHGM$B4X?t@bL@=q(B  @subtitle Risa/Asir 2元分割表HGM関数説明書
 @subtitle 1.0 $BHG(B  @subtitle 1.2 版
 @subtitle 2016 $BG/(B 3 $B7n(B 21 $BF|(B  @subtitle 2019 年 3 月 20 日
   
 @author  by Y.Goto, Y.Tachibana, N.Takayama  @author  by Y.Goto, Y.Tachibana, N.Takayama
 @page  @page
Line 49  Copyright @copyright{} Risa/Asir committers
Line 50  Copyright @copyright{} Risa/Asir committers
 2004--2010. All rights reserved.  2004--2010. All rights reserved.
 @end titlepage  @end titlepage
   
 @comment --- $B$*$^$8$J$$(B ---  @comment --- おまじない ---
 @synindex vr fn  @synindex vr fn
 @comment --- $B$*$^$8$J$$=*$j(B ---  @comment --- おまじない終り ---
   
 @comment --- @node $B$O(B GNU info, HTML $BMQ(B ---  @comment --- @node は GNU info, HTML 用 ---
 @comment --- @node  $B$N0z?t$O(B node-name,  next,  previous,  up ---  @comment --- @node  の引数は node-name,  next,  previous,  up ---
 @node Top,, (dir), (dir)  @node Top,, (dir), (dir)
   
 @comment --- @menu $B$O(B GNU info, HTML $BMQ(B ---  @comment --- @menu は GNU info, HTML 用 ---
 @comment --- chapter $BL>$r@53N$KJB$Y$k(B ---  @comment --- chapter 名を正確に並べる ---
 @comment --- $B$3$NJ8=q$G$O(B chapter XYZ, Chapter Index $B$,$"$k(B.  @comment --- この文書では chapter XYZ, Chapter Index がある.
 @comment ---  Chapter XYZ $B$K$O(B section XYZ$B$K$D$$$F(B, section XYZ$B$K4X$9$k4X?t$,$"$k(B.  @comment ---  Chapter XYZ には section XYZについて, section XYZに関する関数がある.
 @menu  @menu
 * 2$B85J,3dI=(BHGM$B$N4X?t@bL@=q$K$D$$$F(B::  * 2元分割表HGMの関数説明書について::
 * 2$B85J,3dI=(BHGM$B$N4X?t(B::  * 2元分割表HGMの関数::
   * modular計算
 * Index::  * Index::
 @end menu  @end menu
   
 @comment --- chapter $B$N3+;O(B ---  @comment --- chapter の開始 ---
 @comment --- $B?F(B chapter $BL>$r@53N$K(B. $B?F$,$J$$>l9g$O(B Top ---  @comment --- 親 chapter 名を正確に. 親がない場合は Top ---
 @node 2$B85J,3dI=(BHGM$B$N4X?t@bL@=q$K$D$$$F(B,,, Top  @node 2元分割表HGMの関数説明書について,,, Top
 @chapter 2$B85J,3dI=(BHGM$B$N4X?t@bL@=q$K$D$$$F(B  @chapter 2元分割表HGMの関数説明書について
   
 $B$3$N@bL@=q$G$O(B  この説明書では
 HGM(holonomic gradient method) $B$rMQ$$$?(B2$B85J,3dI=$N4X?t$K$D$$$F@bL@$9$k(B.  HGM(holonomic gradient method) を用いた2元分割表の関数について説明する.
 ChangeLog $B$N9`L\$O(B www.openxm.org $B$N(B cvsweb $B$G(B  ChangeLog の項目は www.openxm.org の cvsweb で
 $B%=!<%9%3!<%I$rFI$`;~$N=u$1$K$J$k>pJs$,=q$+$l$F$$$k(B.  ソースコードを読む時の助けになる情報が書かれている.
   このパッケージは下記のようにロードする.
 $BK\J8Cf$G0zMQ$7$F$$$kJ88%$rNs5s$9$k(B.  @example
   load("gtt_ekn3.rr");
   @end example
   gtt_ekn3.rr は gtt_ekn.rr を置き換える大きく改良されたパッケージである.
   以下のモジュール名 gtt_ekn はすべて gtt_ekn3 と読み替えてほしい.
   @noindent
   最新版の asir-contrib package を取得するには, 下記のように更新関数を呼び出す.
   @example
   import("names.rr");
   asir_contrib_update(|update=1);
   @end example
   @noindent
   本文中で引用している文献を列挙する.
 @itemize @bullet  @itemize @bullet
 @item [GM2016]  @item [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)  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)
 @item [T2016]  @item [T2016]
 Y.Tachibana, $B:9J,%[%m%N%_%C%/8{G[K!$N%b%8%e%i!<%a%=%C%I$K$h$k7W;;$N9bB.2=(B,  Y.Tachibana, 差分ホロノミック勾配法のモジュラーメソッドによる計算の高速化,
 2016, $B?@8MBg3X=$;NO@J8(B.  2016, 神戸大学修士論文.
 @item [GTT2016]  @item [GTT2016]
 Y.Goto, Y.Tachibana, N.Takayama, 2$B85J,3dI=$KBP$9$k:9J,%[%m%N%_%C%/8{G[K!$N<BAu(B,  Y.Goto, Y.Tachibana, N.Takayama, 2元分割表に対する差分ホロノミック勾配法の実装,
 $B?tM}8&9V5fO?(B($B7G:\M=Dj(B).  数理研講究録.
   @item [TGKT]
   Y.Tachibana, Y.Goto, T.Koyama, N.Takayama,
   Holonomic Gradient Method for Two Way Contingency Tables,
   arxiv:1803.04170
 @item [TKT2015]  @item [TKT2015]
 N.Takayama, S.Kuriki, A.Takemura,  N.Takayama, S.Kuriki, A.Takemura,
          $A$-hypergeometric distributions and Newton polytopes.           $A$-hypergeometric distributions and Newton polytopes.
          arxiv:1510.02269           arxiv:1510.02269
 @end itemize  @end itemize
   
 $B$3$N%^%K%e%"%k$G@bL@$9$k4X?t$rMQ$$$?%W%m%0%i%`Nc$O(B  このマニュアルで説明する関数を用いたプログラム例は
 gtt_ekn/test-t1.rr  gtt_ekn/test-t1.rr
 $B$J$I(B.  など.
   
 @node 2$B85J,3dI=(BHGM$B$N4X?t(B,,, Top  
 @chapter 2$B85J,3dI=(BHGM$B$N4X?t(B  
   
 @comment --- section ``$B<B83E*4X?t(B'' $B$N(B subsection xyz_abc  @node 2元分割表HGMの関数,,, Top
 @comment --- subsection xyz_pqr xyz_stu $B$,$"$k(B.  @chapter 2元分割表HGMの関数
   
   @comment --- section ``実験的関数'' の subsection xyz_abc
   @comment --- subsection xyz_pqr xyz_stu がある.
 @menu  @menu
 * gtt_ekn.gmvector::  * gtt_ekn.gmvector::
 * gtt_ekn.nc::  * gtt_ekn.nc::
Line 109  gtt_ekn/test-t1.rr 
Line 128  gtt_ekn/test-t1.rr 
 * gtt_ekn.expectation::  * gtt_ekn.expectation::
 * gtt_ekn.setup::  * gtt_ekn.setup::
 * gtt_ekn.upAlpha::  * gtt_ekn.upAlpha::
   * gtt_ekn.cmle::
   * gtt_ekn.set_debug_level::
   * gtt_ekn.contiguity_mat_list_2::
   * gtt_ekn.show_path::
   * gtt_ekn.get_svalue::
   * gtt_ekn.assert1::
   * gtt_ekn.assert2::
   * gtt_ekn.prob2::
 @end menu  @end menu
   
 @node $BD64v2?4X?t(BE(k,n),,, 2$B85J,3dI=(BHGM$B$N4X?t(B  @node 超幾何関数E(k,n),,, 2元分割表HGMの関数
 @section $BD64v2?4X?t(BE(k,n)  @section 超幾何関数E(k,n)
   
 @comment **********************************************************  @comment **********************************************************
 @comment --- $B"~"~"~"~(B  $B$N@bL@(B  @comment --- ◯◯◯◯  の説明
 @comment --- $B8D!9$N4X?t$N@bL@$N3+;O(B ---  @comment --- 個々の関数の説明の開始 ---
 @comment --- section $BL>$r@53N$K(B ---  @comment --- section 名を正確に ---
 @node gtt_ekn.gmvector,,, $BD64v2?4X?t(BE(k,n)  @node gtt_ekn.gmvector,,, 超幾何関数E(k,n)
 @subsection @code{gtt_ekn.gmvector}  @subsection @code{gtt_ekn.gmvector}
 @comment --- $B:w0zMQ%-!<%o!<%I(B  @comment --- 索引用キーワード
 @findex gtt_ekn.gmvector  @findex gtt_ekn.gmvector
   
 @table @t  @table @t
 @item gtt_ekn.gmvector(@var{beta},@var{p})  @item gtt_ekn.gmvector(@var{beta},@var{p})
 :: $B<~JUOB(B @var{beta}, $B%;%k$N3NN((B @var{p} $B$NFs85J,3dI=$KIU?o$9$kD64v2?4X?t(B  :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表に付随する超幾何関数
 E(k,n) $B$NCM$*$h$S$=$NHyJ,$NCM$rLa$9(B.  E(k,n) の値およびその微分の値を戻す.
 @item gtt_ekn.ekn_cBasis_2(@var{beta},@var{p})  @item gtt_ekn.ekn_cBasis_2(@var{beta},@var{p})
 $B$NJLL>$G$"$k(B.  の別名である.
 @end table  @end table
   
 @comment --- $B0z?t$N4JC1$J@bL@(B ---  $B0J2<$^$@=q$$$F$J$$(B.  @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
 @table @var  @table @var
 @item return  @item return
 $B%Y%/%H%k(B, $BD64v2?4X?t$NCM$H$=$NHyJ,(B. $B>\$7$/$O2<5-(B.  ベクトル, 超幾何関数の値とその微分. 詳しくは下記.
 @item beta  @item beta
 $B9TOB(B, $BNsOB$N%j%9%H(B. $B@.J,$O$9$Y$F@5$G$"$k$3$H(B.  行和, 列和のリスト. 成分はすべて正であること.
 @item p  @item p
 $BFs85J,3dI=$N%;%k$N3NN($N%j%9%H(B  二元分割表のセルの確率のリスト
 @end table  @end table
   
 @comment --- $B$3$3$G4X?t$N>\$7$$@bL@(B ---  @comment --- ここで関数の詳しい説明 ---
 @comment --- @itemize$B!A(B@end itemize $B$O2U>r=q$-(B ---  @comment --- @itemize〜@end itemize は箇条書き ---
 @comment --- @bullet $B$O9uE@IU$-(B ---  @comment --- @bullet は黒点付き ---
 @itemize @bullet  @itemize @bullet
 @item  @item
 gmvector $B$O(B Gauss-Manin vector $B$NN,$G$"$k(B [GM2016].  gmvector は Gauss-Manin vector の略である [GM2016].
 @item  @item
 gmvector $B$NLa$jCM$O(B [GM2016] $B$N#4>O$GDj5A$5$l$F$$$k%Y%/%H%k(B F $B$G$"$k(B.  gmvector の戻り値は
 $B$?$@$7Bh0l@.J,$,(B [GM2016] $B$N#6>O$GDj5A$5$l$F$$$k5i?t(B S $B$NCM$HEy$7$/(B  [GM2016] の 6章 p.23 のベクトル Sである.
 $B$J$k$h$&$K%9%+%i!<G\$5$l$F$$$k(B.  これは
   [GM2016] の4章で定義されているベクトル F の定数倍であり,
   その定数は
   第一成分が [GM2016] の6章で定義されている級数 S の値と等しく
   なるように決められている.
 @item  @item
  r1 x r2 $BJ,3dI=$r9M$($k(B.   r1 x r2 分割表を考える.
  m+1=r1, n+1=r2 $B$H$*$/(B.   m+1=r1, n+1=r2 とおく.
  $B@55,2=Dj?t(B Z $B$OJ,3dI=(B u $B$r(B (m+1) $B!_(B (n+1) $B9TNs$H$9$k$H$-(B p^u/u! $B$NOB$G$"$k(B.   正規化定数 Z は分割表 u を (m+1) × (n+1) 行列とするとき p^u/u! の和である.
  $B$3$3$GOB$O9TOBNsOB$,(B @var{beta} $B$G$"$k$h$&$J(B u $BA4BN$G$H$k(B   ここで和は行和列和が @var{beta} であるような u 全体でとる
  [TKT2015], [GM2016].   [TKT2015], [GM2016].
  S $B$O$3$N5i?t$N(B p $B$r(B   S はこの多項式 Z の p を
 @verbatim  @verbatim
   [[1,y11,...,y1n],    [[1,y11,...,y1n],
    [1,y21,...,y2n],...,     [1,y21,...,y2n],...,
    [1,ym1, ...,ymn],     [1,ym1, ...,ymn],
    [1,1, ..., 1]]     [1,1, ..., 1]]
 @end verbatim  @end verbatim
 $B!!(B(1 $B$,(B L $B;z7?$KJB$V(B),   (1 が L 字型に並ぶ),
 $B$H@55,2=$7$?5i?t$G$"$k(B.  と正規化した級数である.
 @item  @item
 2x(n+1)$BJ,3dI=$G(B, gmvector $B$NLa$jCM$r(B Lauricella  F_D $B$G=q$/$3$H$,(B  2x(n+1)分割表で, gmvector の戻り値を Lauricella  F_D で書くことが
 $B0J2<$N$h$&$K$G$-$k(B  以下のようにできる
 (b[2][1]-b[1][1] >= 0 $B$N>l9g(B).  (b[2][1]-b[1][1] >= 0 の場合).
 $B$3$3$G(B b[1][1], b[1][2] $B$O(B, $B$=$l$>$l(B 1 $B9TL\$N9TOB(B, 2 $B9TL\$N9TOB(B,  ここで b[1][1], b[1][2] は, それぞれ 1 行目の行和, 2 行目の行和,
 b[2][i] $B$O(B i $BNsL\$NNsOB$G$"$k(B.  b[2][i] は i 列目の列和である.
 @comment ekn/Talks/2015-12-3-goto.tex  @comment ekn/Talks/2015-12-3-goto.tex
 @verbatim  @verbatim
 S=F_D(-b[1,1], [-b[2,2],...,-b[2,n+1]], b[2,1]-b[1,1]+1 ; y)/C,  S=F_D(-b[1,1], [-b[2,2],...,-b[2,n+1]], b[2,1]-b[1,1]+1 ; y)/C,
 @end verbatim  @end verbatim
 C=b[1,1]! b[2,2]! ... b[2][n+1]! (b[2,1]-b[1,1])!  C=b[1,1]! b[2,2]! ... b[2,n+1]! (b[2,1]-b[1,1])!
 $B$H$*$/(B.  とおく.
 1/C $B$O(B L $B;z7?$NJ,3dI=(B  1/C は L 字型の分割表
 @verbatim  @verbatim
 [[b[1,1],       0,      ..., 0       ],  [[b[1,1],       0,      ..., 0       ],
  [b[2,1]-b[1,1],b[2,2], ..., b[2,n+1]]]   [b[2,1]-b[1,1],b[2,2], ..., b[2,n+1]]]
 @end verbatim  @end verbatim
 $B$KBP1~(B.  に対応.
 gmvector $B$O(B  gmvector は
 @verbatim  @verbatim
 [S,(y11/a2) d_11 S,(y12/a3) d_12 S, ..., (y1n/a_(n+1)) d_1n S]  [S,(y11/a2) d_11 S,(y12/a3) d_12 S, ..., (y1n/a_(n+1)) d_1n S]
 @end verbatim  @end verbatim
 $B$G$"$k(B.  である.
 $B$3$3$G(B d_ij $B$O(B yij $B$K$D$$$F$NHyJ,(B,  ここで d_ij は yij についての微分,
 @verbatim  @verbatim
   [a0,     a1, ...                      ,a_(n+2)]    [a0,     a1, ...                      ,a_(n+2)]
 = [-b[1,2],-b[1,1],b[2,2], ..., b[2,n+1],b[2,1]]  = [-b[1,2],-b[1,1],b[2,2], ..., b[2,n+1],b[2,1]]
 @end verbatim  @end verbatim
 $B$G$"$k(B.  である.
 @item  @item
 $B<~JUOB(B @var{beta}$B$N;~$N@55,2=Dj?t$N%;%k3NN((B @var{p} $B$KBP$9$kCM$O(B $BB?9`<0$KB`2=$7$?(B E(k,n) $B$NCM$GI=8=$G$-$k(B. $BJ88%(B [TKT2015], [GM2016] $B;2>H(B.  周辺和 @var{beta}の時の正規化定数のセル確率 @var{p} に対する値は 多項式に退化した E(k,n) の値で表現できる. 文献 [TKT2015], [GM2016] 参照.
 @item  @item
 option crt=1 (crt = Chinese remainder theorem) $B$rM?$($k$H(B, $BJ,;67W;;$r$*$3$J$&(B  option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう
 [T2016].  [T2016].
 $BJ,;67W;;MQ$N3F<o%Q%i%a!<%?$N@_Dj$O(B  分散計算用の各種パラメータの設定は
 gtt_ekn.setup $B$G9T$J$&(B.  gtt_ekn.setup で行なう.
 @end itemize  @end itemize
   
 @comment --- @example$B!A(B@end example $B$O<B9TNc$NI=<((B ---  @comment --- @example〜@end example は実行例の表示 ---
 $BNc(B: $B<!$O(B2 x 2 $BJ,3dI=$G9TOB$,(B [5,1],  $BNsOB$,(B [3,3], $B3F%;%k$N3NN($,(B  例: 次は2 x 2 分割表で行和が [5,1],  列和が [3,3], 各セルの確率が
 [[1/2,1/3],[1/7,1/5]] $B$N>l9g$N(B gmvector $B$NCM$G$"$k(B.  [[1/2,1/3],[1/7,1/5]] の場合の gmvector の値である.
 @example  @example
 [3000] load("gtt_ekn.rr");  [3000] load("gtt_ekn.rr");
 [3001] ekn_gtt.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])  [3001] ekn_gtt.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])
Line 214  gtt_ekn.setup $B$G9T$J$&(B.
Line 245  gtt_ekn.setup $B$G9T$J$&(B.
 [200/9261]  [200/9261]
 @end example  @end example
   
 $B;29M(B: 2 x m $BJ,3dI=(B(Lauricella FD)$B$K$D$$$F$O%Q%C%1!<%8(B tk_fd $B$G$b2<5-$N$h$&$KF1Ey$J(B  例: N を2以上の自然数とする時, Gauss の超幾何関数(この場合は多項式となる)
 $B7W;;$,$G$-$k(B.  F(-36N,-11N,2N,(1-1/N)/56) の値は T3 に代入される ( [TGKT] ).
 $B<iHwHO0O$N0[$J$k%W%m%0%i%`F1;N$NHf3S(B, debug $BMQ;29M(B.  @comment ekn/Prog2/2x2.rr
 @example  @example
   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;
   @end example
   ちなみに同じ値を Mathematica に計算させるには
   @example
   n=2; Hypergeometric2F1[-36*n,-11*n,2*n,(1-1/n)/56]
   @end example
   
   参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
   計算ができる.
   守備範囲の異なるプログラム同士の比較, debug 用参考.
   @example
 [3080] import("tk_fd.rr");  [3080] import("tk_fd.rr");
 [3081] A=tk_fd.marginal2abc([4,5],[2,4,3]);  [3081] A=tk_fd.marginal2abc([4,5],[2,4,3]);
 [-4,[-4,-3],-1]  // 2$BJQ?t(B FD $B$N%Q%i%a!<%?(B. a,[b1,b2],c  [-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]);  [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 ]  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
 [4483/124416,[ 1961/15552 185/1728 ],  [4483/124416,[ 1961/15552 185/1728 ],
  [ 79/288 259/864 ]   [ 79/288 259/864 ]
  [ 259/864 47/288 ]]   [ 259/864 47/288 ]]
 // $BLaCM$O(B [F=F_D, gradient(F), Hessian(F)]  // 戻値は [F=F_D, gradient(F), Hessian(F)]
   
 // ekn_gt $B$G$NNc$HF1$8%Q%i%a!<%?(B.  // ekn_gt での例と同じパラメータ.
 [3543] A=tk_fd.marginal2abc([5,1],[3,3]);  [3543] A=tk_fd.marginal2abc([5,1],[3,3]);
 [-5,[-3],-1]  [-5,[-3],-1]
 [3544] tk_fd.fd_hessian2(A[0],A[1],A[2],[(1/3)*(1/7)/((1/2)*(1/5))]);  [3544] tk_fd.fd_hessian2(A[0],A[1],A[2],[(1/3)*(1/7)/((1/2)*(1/5))]);
Line 236  Computing Dmat(ca) for parameters B=[-3],X=[ 10/21 ]
Line 281  Computing Dmat(ca) for parameters B=[-3],X=[ 10/21 ]
 [775/27783,[ 20/147 ],[ 17/42 ]]  [775/27783,[ 20/147 ],[ 17/42 ]]
 @end example  @end example
   
 $B;29M(B: $B0lHL$N(B A $BJ,I[$N@55,2=Dj?t$K$D$$$F$N(B Hessian $B$N7W;;$O<B83E*(B package ot_hessian_ahg.rr  参考: 一般の A 分布の正規化定数についての Hessian の計算は実験的 package ot_hessian_ahg.rr
 $B$G<BAu$N%F%9%H$,$5$l$F$$$k(B. ($B$3$l$O$^$@L$40@.$N%F%9%HHG$J$N$G=PNO7A<0Ey$b>-MhE*$K$OJQ99$5$l$k(B.)  で実装のテストがされている. (これはまだ未完成のテスト版なので出力形式等も将来的には変更される.)
 @example  @example
 import("ot_hgm_ahg.rr");  import("ot_hgm_ahg.rr");
 import("ot_hessian_ahg.rr");  import("ot_hessian_ahg.rr");
Line 267  def  htest4() @{
Line 312  def  htest4() @{
  [(b1^2+(-2*b0-1)*b1+b0^2+b0)/(x0^2),[]],[(x6)/(x0),[6,0]],[(x4)/(x0),[4,0]]]   [(b1^2+(-2*b0-1)*b1+b0^2+b0)/(x0^2),[]],[(x6)/(x0),[6,0]],[(x4)/(x0),[4,0]]]
 @end example  @end example
   
 @comment --- $B;2>H(B($B%j%s%/(B)$B$r=q$/(B ---  @comment --- 参照(リンク)を書く ---
 @table @t  @table @t
 @item $B;2>H(B  @item 参照
 @ref{gtt_ekn.setup}  @ref{gtt_ekn.setup}
 @ref{gtt_ekn.pfaffian_basis}  @ref{gtt_ekn.pfaffian_basis}
 @end table  @end table
   
 @comment --- ChangeLog $B$r=q$/(B. $B%=!<%9%3!<%I$N0LCV(B. $BJQ99F|;~(B $B$J$I(B CVS$B%5!<%P$r8+$k$?$a(B  @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
 @noindent  @noindent
 ChangeLog  ChangeLog
 @itemize @bullet  @itemize @bullet
 @item  @item
  $B$3$N4X?t$O(B   この関数は
 [GM2016] $B$N%"%k%4%j%:%`$*$h$S(B  [GM2016] のアルゴリズムおよび
 [T2016] $B$K$h$k(B modular method $B$rMQ$$$?9bB.2=$r<BAu$7$?$b$N$G$"$k(B.  [T2016] による modular method を用いた高速化を実装したものである.
 @item  @item
  $BJQ99$r<u$1$?%U%!%$%k$O(B   変更を受けたファイルは
  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_pfaffian_8.rr   OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_pfaffian_8.rr
 @end itemize  @end itemize
   
   
 @comment **********************************************************  @comment **********************************************************
 @node gtt_ekn.nc,,, $BD64v2?4X?t(BE(k,n)  @node gtt_ekn.nc,,, 超幾何関数E(k,n)
 @subsection @code{gtt_ekn.nc}  @subsection @code{gtt_ekn.nc}
 @comment --- $B:w0zMQ%-!<%o!<%I(B  @comment --- 索引用キーワード
 @findex gtt_ekn.nc  @findex gtt_ekn.nc
   
 @table @t  @table @t
 @item gtt_ekn.nc(@var{beta},@var{p})  @item gtt_ekn.nc(@var{beta},@var{p})
 :: $B<~JUOB(B @var{beta}, $B%;%k$N3NN((B @var{p} $B$NFs85J,3dI=$N>r7oIU$-3NN($N@55,2=Dj?t(B Z  :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z
 $B$*$h$S$=$NHyJ,$NCM$rLa$9(B.  およびその微分の値を戻す.
 @end table  @end table
   
 @comment --- $B0z?t$N4JC1$J@bL@(B ---  $B0J2<$^$@=q$$$F$J$$(B.  @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
 @table @var  @table @var
 @item return  @item return
 $B%Y%/%H%k(B [Z,[[d_11 Z, d_12 Z, ...], ..., [d_m1 Z, d_m2 Z, ...., d_mn Z]]]  ベクトル [Z,[[d_11 Z, d_12 Z, ...], ..., [d_m1 Z, d_m2 Z, ...., d_mn Z]]]
 @item beta  @item beta
 $B9TOB(B, $BNsOB$N%j%9%H(B. $B@.J,$O$9$Y$F@5$G$"$k$3$H(B.  行和, 列和のリスト. 成分はすべて正であること.
 @item p  @item p
 $BFs85J,3dI=$N%;%k$N3NN($N%j%9%H(B  二元分割表のセルの確率のリスト
 @end table  @end table
   
 @comment --- $B$3$3$G4X?t$N>\$7$$@bL@(B ---  @comment --- ここで関数の詳しい説明 ---
 @comment --- @itemize$B!A(B@end itemize $B$O2U>r=q$-(B ---  @comment --- @itemize〜@end itemize は箇条書き ---
 @comment --- @bullet $B$O9uE@IU$-(B ---  @comment --- @bullet は黒点付き ---
 @itemize @bullet  @itemize @bullet
 @item  @item
  r1 x r2 $BJ,3dI=$r9M$($k(B.   r1 x r2 分割表を考える.
  m=r1, n=r2 $B$H$*$/(B.   m=r1, n=r2 とおく.
  $B@55,2=Dj?t(B Z $B$OJ,3dI=(B u $B$r(B m $B!_(B n $B9TNs$H$9$k$H$-(B p^u/u! $B$NOB$G$"$k(B.   正規化定数 Z は分割表 u を m × n 行列とするとき p^u/u! の和である.
  $B$3$3$GOB$O9TOBNsOB$,(B @var{beta} $B$G$"$k$h$&$J(B u $BA4BN$G$H$k(B   ここで和は行和列和が @var{beta} であるような u 全体でとる
  [TKT2015], [GM2016].   [TKT2015], [GM2016].
  p^u $B$O(B p_ij^u_ij $B$N@Q(B, u! $B$O(B u_ij! $B$N@Q$G$"$k(B.   p^u は p_ij^u_ij の積, u! は u_ij! の積である.
  d_ij Z $B$G(B Z $B$NJQ?t(B p_ij $B$K$D$$$F$NJPHyJ,$rI=$9(B.   d_ij Z で Z の変数 p_ij についての偏微分を表す.
 @item  @item
 nc $B$O(B gmvector $B$NCM$r85$K(B, [GM2016] $B$N(B Prop  nc は gmvector の値を元に, [GM2016] の Prop
  7.1 $B$K4p$E$$$F(B Z $B$NCM$r7W;;$9$k(B.   7.1 に基づいて Z の値を計算する.
 @item  @item
 option crt=1 (crt = Chinese remainder theorem) $B$rM?$($k$H(B, $BJ,;67W;;$r$*$3$J$&(B.  option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
 $BJ,;67W;;MQ$N3F<o%Q%i%a!<%?$N@_Dj$O(B  分散計算用の各種パラメータの設定は
 gtt_ekn.setup $B$G9T$J$&(B.  gtt_ekn.setup で行なう.
 @end itemize  @end itemize
   
 @comment --- @example$B!A(B@end example $B$O<B9TNc$NI=<((B ---  @comment --- @example〜@end example は実行例の表示 ---
 $BNc(B: 2x3 $BJ,3dI=$G$N(B Z $B$H$=$NHyJ,$N7W;;(B.  例: 2x3 分割表での Z とその微分の計算.
 @example  @example
 [2237] gtt_ekn.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);  [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 ]  [4483/124416,[ 353/7776 1961/15552 185/1728 ]
 [ 553/20736 1261/15552 1001/13824 ]]  [ 553/20736 1261/15552 1001/13824 ]]
 @end example  @end example
   
 $B;29M(B: 2 x m $BJ,3dI=(B(Lauricella FD)$B$K$D$$$F$O%Q%C%1!<%8(B tk_fd $B$G$b2<5-$N$h$&$KF1Ey$J(B  参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
 $B7W;;$,$G$-$k(B.  計算ができる.
 @example  @example
 [3076] import("tk_fd.rr");  [3076] import("tk_fd.rr");
 [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);  [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
Line 351  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
Line 396  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
 [4483/124416,[[353/7776,1961/15552,185/1728],  [4483/124416,[[353/7776,1961/15552,185/1728],
               [553/20736,1261/15552,1001/13824]]]                [553/20736,1261/15552,1001/13824]]]
 // $BLaCM$O(B [Z, [[d_11 Z, d_12 Z, d_13 Z],  // 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],
 //             [d_21 Z, d_22 Z, d_23 Z]]] $B$NCM(B.  //             [d_21 Z, d_22 Z, d_23 Z]]] の値.
 //           $B$3$3$G(B d_ij $B$O(B i,j $B@.J,$K$D$$$F$NHyJ,$rI=$9(B.  //           ここで d_ij は i,j 成分についての微分を表す.
 @end example  @end example
   
 @comment --- $B;2>H(B($B%j%s%/(B)$B$r=q$/(B ---  @comment --- 参照(リンク)を書く ---
 @table @t  @table @t
 @item $B;2>H(B  @item 参照
 @ref{gtt_ekn.setup}  @ref{gtt_ekn.setup}
 @ref{gtt_ekn.lognc}  @ref{gtt_ekn.lognc}
 @end table  @end table
   
 @comment --- ChangeLog $B$r=q$/(B. $B%=!<%9%3!<%I$N0LCV(B. $BJQ99F|;~(B $B$J$I(B CVS$B%5!<%P$r8+$k$?$a(B  @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
 @noindent  @noindent
 ChangeLog  ChangeLog
 @itemize @bullet  @itemize @bullet
 @item  @item
  $BJQ99$r<u$1$?%U%!%$%k$O(B   変更を受けたファイルは
  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_eval.rr   OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_eval.rr
 @end itemize  @end itemize
   
   
 @comment **********************************************************  @comment **********************************************************
 @node gtt_ekn.lognc,,, $BD64v2?4X?t(BE(k,n)  @node gtt_ekn.lognc,,, 超幾何関数E(k,n)
 @subsection @code{gtt_ekn.lognc}  @subsection @code{gtt_ekn.lognc}
 @comment --- $B:w0zMQ%-!<%o!<%I(B  @comment --- 索引用キーワード
 @findex gtt_ekn.lognc  @findex gtt_ekn.lognc
   
 @table @t  @table @t
 @item gtt_ekn.lognc(@var{beta},@var{p})  @item gtt_ekn.lognc(@var{beta},@var{p})
 :: $B<~JUOB(B @var{beta}, $B%;%k$N3NN((B @var{p} $B$NFs85J,3dI=$N>r7oIU$-3NN($N@55,2=Dj?t(B Z  :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z
 $B$N(B log $B$N6a;wCM$*$h$S$=$NHyJ,$N6a;wCM$rLa$9(B.  の log の近似値およびその微分の近似値を戻す.
 @end table  @end table
   
 @comment --- $B0z?t$N4JC1$J@bL@(B ---  $B0J2<$^$@=q$$$F$J$$(B.  @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
 @table @var  @table @var
 @item return  @item return
 $B%Y%/%H%k(B [log(Z), [[d_11 log(Z), d_12 log(Z), ...], [d_21 log(Z),...], ... ]  ベクトル [log(Z), [[d_11 log(Z), d_12 log(Z), ...], [d_21 log(Z),...], ... ]
 @item beta  @item beta
 $B9TOB(B, $BNsOB$N%j%9%H(B. $B@.J,$O$9$Y$F@5$G$"$k$3$H(B.  行和, 列和のリスト. 成分はすべて正であること.
 @item p  @item p
 $BFs85J,3dI=$N%;%k$N3NN($N%j%9%H(B  二元分割表のセルの確率のリスト
 @end table  @end table
   
 @comment --- $B$3$3$G4X?t$N>\$7$$@bL@(B ---  @comment --- ここで関数の詳しい説明 ---
 @comment --- @itemize$B!A(B@end itemize $B$O2U>r=q$-(B ---  @comment --- @itemize〜@end itemize は箇条書き ---
 @comment --- @bullet $B$O9uE@IU$-(B ---  @comment --- @bullet は黒点付き ---
 @itemize @bullet  @itemize @bullet
 @item  @item
 $B>r7oIU$-:GL`?dDj$KMxMQ$9$k(B [TKT2015].  条件付き最尤推定に利用する [TKT2015].
 @item option crt=1 (crt = Chinese remainder theorem) $B$rM?$($k$H(B, $BJ,;67W;;$r$*$3$J$&(B.  @item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
 $BJ,;67W;;MQ$N3F<o%Q%i%a!<%?$N@_Dj$O(B  分散計算用の各種パラメータの設定は
 gtt_ekn.setup $B$G9T$J$&(B.  gtt_ekn.setup で行なう.
 @end itemize  @end itemize
   
 @comment --- @example$B!A(B@end example $B$O<B9TNc$NI=<((B ---  @comment --- @example〜@end example は実行例の表示 ---
 $BNc(B: 2 $B!_(B 3 $BJ,3dI=$G$NNc(B. $BBh0l@.J,$N$_6a;wCM(B.  例: 2 × 3 分割表での例. 第一成分のみ近似値.
 @example  @example
 [2238] gtt_ekn.lognc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);  [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 ]  [-3.32333832422461674630,[ 5648/4483 15688/4483 13320/4483 ]
 [ 3318/4483 10088/4483 9009/4483 ]]  [ 3318/4483 10088/4483 9009/4483 ]]
 @end example  @end example
   
 $B;29M(B: 2 x m $BJ,3dI=(B(Lauricella FD)$B$K$D$$$F$O%Q%C%1!<%8(B tk_fd $B$G$b2<5-$N$h$&$KF1Ey$J(B  参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
 $B7W;;$,$G$-$k(B.  計算ができる.
 @example  @example
 [3076] import("tk_fd.rr");  [3076] import("tk_fd.rr");
 [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);  [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
Line 427  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/
Line 472  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/
 [-3.32333832422461674639485797719209322217260539267246045320,  [-3.32333832422461674639485797719209322217260539267246045320,
  [[1.2598706, 3.499442, 2.971224],   [[1.2598706, 3.499442, 2.971224],
   [0.7401293, 2.250278, 2.009591]]]    [0.7401293, 2.250278, 2.009591]]]
 // $BLaCM$O(B [log(Z),  // 戻値は [log(Z),
 //          [[d_11 log(Z), d_12 log(Z), d_13 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)]]]  //           [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]]
 // $B$N6a;wCM(B.  // の近似値.
 @end example  @end example
   
 @comment --- $B;2>H(B($B%j%s%/(B)$B$r=q$/(B ---  @comment --- 参照(リンク)を書く ---
 @table @t  @table @t
 @item $B;2>H(B  @item 参照
 @ref{gtt_ekn.setup}  @ref{gtt_ekn.setup}
 @ref{gtt_ekn.nc}  @ref{gtt_ekn.nc}
 @end table  @end table
   
 @comment --- ChangeLog $B$r=q$/(B. $B%=!<%9%3!<%I$N0LCV(B. $BJQ99F|;~(B $B$J$I(B CVS$B%5!<%P$r8+$k$?$a(B  @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
 @noindent  @noindent
 ChangeLog  ChangeLog
 @itemize @bullet  @itemize @bullet
 @item  @item
  $BJQ99$r<u$1$?%U%!%$%k$O(B   変更を受けたファイルは
  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.   OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.
 @end itemize  @end itemize
   
 @comment **********************************************************  @comment **********************************************************
 @node gtt_ekn.expectation,,, $BD64v2?4X?t(BE(k,n)  @node gtt_ekn.expectation,,, 超幾何関数E(k,n)
 @subsection @code{gtt_ekn.expectation}  @subsection @code{gtt_ekn.expectation}
 @comment --- $B:w0zMQ%-!<%o!<%I(B  @comment --- 索引用キーワード
 @findex gtt_ekn.expectation  @findex gtt_ekn.expectation
   
 @table @t  @table @t
 @item gtt_ekn.expectation(@var{beta},@var{p})  @item gtt_ekn.expectation(@var{beta},@var{p})
 :: $B<~JUOB(B @var{beta}, $B%;%k$N3NN((B @var{p} $B$NFs85J,3dI=$N4|BTCM$r7W;;$9$k(B.  :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の期待値を計算する.
 @end table  @end table
   
 @comment --- $B0z?t$N4JC1$J@bL@(B ---  $B0J2<$^$@=q$$$F$J$$(B.  @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
 @table @var  @table @var
 @item return  @item return
 $BFs85J,3dI=$N3F%;%k$N4|BTCM$N%j%9%H(B.  二元分割表の各セルの期待値のリスト.
 @item beta  @item beta
 $B9TOB(B, $BNsOB$N%j%9%H(B. $B@.J,$O$9$Y$F@5$G$"$k$3$H(B.  行和, 列和のリスト. 成分はすべて正であること.
 @item p  @item p
 $BFs85J,3dI=$N%;%k$N3NN($N%j%9%H(B  二元分割表のセルの確率のリスト
 @end table  @end table
   
 @comment --- $B$3$3$G4X?t$N>\$7$$@bL@(B ---  @comment --- ここで関数の詳しい説明 ---
 @comment --- @itemize$B!A(B@end itemize $B$O2U>r=q$-(B ---  @comment --- @itemize〜@end itemize は箇条書き ---
 @comment --- @bullet $B$O9uE@IU$-(B ---  @comment --- @bullet は黒点付き ---
 @itemize @bullet  @itemize @bullet
 @item  @item
 [GM2016] $B$N(B Algorithm 7.8 $B$N<BAu(B.  [GM2016] の Algorithm 7.8 の実装.
 @item option crt=1 (crt = Chinese remainder theorem) $B$rM?$($k$H(B, $BJ,;67W;;$r$*$3$J$&(B.  @item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
 $BJ,;67W;;MQ$N3F<o%Q%i%a!<%?$N@_Dj$O(B  分散計算用の各種パラメータの設定は
 gtt_ekn.setup $B$G9T$J$&(B.  gtt_ekn.setup で行なう.
 @item option index $B$rM?$($k$H(B, $B;XDj$5$l$?@.J,$N4|BTCM$N$_7W;;$9$k(B.  @item option index を与えると, 指定された成分の期待値のみ計算する.
 $B$?$H$($P(B 2 x 2 $BJ,3dI=$G(B index=[[0,0],[1,1]] $B$H;XDj$9$k$H(B, 1 $B$N$"$k@.J,$N4|BTCM$N$_7W;;$9$k(B.  たとえば 2 x 2 分割表で index=[[0,0],[1,1]] と指定すると, 1 のある成分の期待値のみ計算する.
 @end itemize  @end itemize
   
 @comment --- @example$B!A(B@end example $B$O<B9TNc$NI=<((B ---  @comment --- @example〜@end example は実行例の表示 ---
   
 2$B!_(B2, 3$B!_(B3 $B$NJ,3dI=$N4|BTCM7W;;Nc(B.  2×2, 3×3 の分割表の期待値計算例.
 @example  @example
 [2235] gtt_ekn.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);  [2235] gtt_ekn.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);
 [ 2/3 1/3 ]  [ 2/3 1/3 ]
Line 503  gtt_ekn.setup $B$G9T$J$&(B.
Line 548  gtt_ekn.setup $B$G9T$J$&(B.
                                         737732646860489910/147000422096729819 ]                                          737732646860489910/147000422096729819 ]
 @end example  @end example
   
 $B;29M(B: 2 x m $BJ,3dI=(B(Lauricella FD)$B$K$D$$$F$O%Q%C%1!<%8(B tk_fd $B$G$b2<5-$N$h$&$KF1Ey$J(B  参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
 $B7W;;$,$G$-$k(B.  計算ができる.
 @example  @example
 [3076] import("tk_fd.rr");  [3076] import("tk_fd.rr");
 [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);  [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
Line 515  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
Line 560  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
 [[5648/4483,7844/4483,4440/4483],  [[5648/4483,7844/4483,4440/4483],
  [3318/4483,10088/4483,9009/4483]]   [3318/4483,10088/4483,9009/4483]]
 // $B3F%;%k$N4|BTCM(B.  // 各セルの期待値.
 @end example  @end example
   
 $B;29M(B: $B0lHL$N(B A $BJ,I[$N7W;;$O(B ot_hgm_ahg.rr. $B$^$@<B83E*$J$?$a(B, module $B2=$5$l$F$$$J$$(B.  参考: 一般の A 分布の計算は ot_hgm_ahg.rr. まだ実験的なため, module 化されていない.
 ot_hgm_ahg.rr $B$K$D$$$F$N;29MJ88%(B:  ot_hgm_ahg.rr についての参考文献:
 K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Systems II --- Holonomic Gradient Method, arxiv:1505.02947  K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Systems II --- Holonomic Gradient Method, arxiv:1505.02947
 @example  @example
 [3237] import("ot_hgm_ahg.rr");  [3237] import("ot_hgm_ahg.rr");
 // 2 x 2 $BJ,3dI=(B.  // 2 x 2 分割表.
 [3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],  [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);          [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
 oohg_native=0, oohg_curl=1  oohg_native=0, oohg_curl=1
 [1376777025/625400597,1750225960/625400597,  [1376777025/625400597,1750225960/625400597,
  2375626557/625400597,3252978816/625400597]   2375626557/625400597,3252978816/625400597]
 // 2 x 2 $BJ,3dI=$N4|BTCM(B.  // 2 x 2 分割表の期待値.
   
 // 2 x 3 $BJ,3dI=(B.  // 2 x 3 分割表.
 [3238] hgm_ahg_expected_values_contiguity(  [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]],   [[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);   [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]  [5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
 // 2 x 3 $BJ,3dI=$N4|BTCM(B. $B>e$HF1$8LdBj(B.  // 2 x 3 分割表の期待値. 上と同じ問題.
 @end example  @end example
   
 3 x 3 $BJ,3dI=(B. $B9=B$E*(B0$B$,0l$D(B.  3 x 3 分割表. 構造的0が一つ.
 @example  @example
 /*  /*
   dojo, p.221 $B$N%G!<%?(B.  $B@.@S(B3$B0J2<$N@8EL$O=8$a$F$R$H$D$K(B.    dojo, p.221 のデータ.  成績3以下の生徒は集めてひとつに.
   2 1 1    2 1 1
   8 3 3    8 3 3
   0 2 6    0 2 6
   
   row sum: 4,14,8    row sum: 4,14,8
   column sum: 10,6,10    column sum: 10,6,10
   0 $B$r0l$D4^$`$N$G(B, (3,6) $B7?$N(B A $B$+$i(B 7 $BNsL\$rH4$/(B.    0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
 */  */
   
 A=[[0,0,0,1,1,1, 0,0],  A=[[0,0,0,1,1,1, 0,0],
Line 559  A=[[0,0,0,1,1,1, 0,0],
Line 604  A=[[0,0,0,1,1,1, 0,0],
    [0,0,1,0,0,1, 0,1]];     [0,0,1,0,0,1, 0,1]];
 B=[14,8,10,6,10];  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],  hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1],
 $B!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!(B[x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);                  [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);
   
 // $BEz(B.  // 答.
 [14449864949304/9556267369631,  [14449864949304/9556267369631,
  10262588586540/9556267369631, 13512615942680/9556267369631,   10262588586540/9556267369631, 13512615942680/9556267369631,
  81112808747006/9556267369631,   81112808747006/9556267369631,
Line 570  hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/
Line 615  hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/
  25258717886900/9556267369631,51191421070148/9556267369631]   25258717886900/9556267369631,51191421070148/9556267369631]
 @end example  @end example
   
 3 x 3 $BJ,3dI=(B.  3 x 3 分割表.
 @example  @example
 /*  /*
  $B>e$N%G!<%?$G(B 0 $B$r(B 1 $B$KJQ99(B.   上のデータで 0 を 1 に変更.
   2 1 1    2 1 1
   8 3 3    8 3 3
   1 2 6    1 2 6
Line 590  B=[14,9,11,6,10];
Line 635  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],  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);                                [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);
   
 // $B4|BTCM(B, $BEz(B.   x9 $B$r;XDj$7$F$$$J$$$N$G(B, 9$BHVL\$N4|BTCM$O=PNO$7$F$J$$(B.  // 期待値, 答.   x9 を指定していないので, 9番目の期待値は出力してない.
 [207017568232262040/147000422096729819,  [207017568232262040/147000422096729819,
  163140751505489940/147000422096729819,217843368649167296/147000422096729819,   163140751505489940/147000422096729819,217843368649167296/147000422096729819,
  1185482401011137878/147000422096729819,   1185482401011137878/147000422096729819,
  358095302885438604/147000422096729819,514428205457640984/147000422096729819,   358095302885438604/147000422096729819,514428205457640984/147000422096729819,
  224504673820628091/147000422096729819,360766478189450370/147000422096729819]   224504673820628091/147000422096729819,360766478189450370/147000422096729819]
   
 // Z $B$d$=$NHyJ,$N7W;;$O(B hgm_ahg_contiguity $B4X?t$,$*$3$J$&$,(B, $B$3$l$N4J0W%$%s%?!<%U%'!<%9$O(B  // Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
 // $B$^$@=q$$$F$J$$(B.  // まだ書いてない.
 @end example  @end example
   
   
   
 @comment --- $B;2>H(B($B%j%s%/(B)$B$r=q$/(B ---  @comment --- 参照(リンク)を書く ---
 @table @t  @table @t
 @item $B;2>H(B  @item 参照
 @ref{gtt_ekn.setup}  @ref{gtt_ekn.setup}
 @ref{gtt_ekn.nc}  @ref{gtt_ekn.nc}
 @end table  @end table
   
 @comment --- ChangeLog $B$r=q$/(B. $B%=!<%9%3!<%I$N0LCV(B. $BJQ99F|;~(B $B$J$I(B CVS$B%5!<%P$r8+$k$?$a(B  @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
 @noindent  @noindent
 ChangeLog  ChangeLog
 @itemize @bullet  @itemize @bullet
 @item  @item
  $BJQ99$r<u$1$?%U%!%$%k$O(B   変更を受けたファイルは
  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.   OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.
 @end itemize  @end itemize
   
   
 @comment **********************************************************  @comment **********************************************************
 @comment --- $B"~"~"~"~(B  $B$N@bL@(B  @comment --- ◯◯◯◯  の説明
 @comment --- $B8D!9$N4X?t$N@bL@$N3+;O(B ---  @comment --- 個々の関数の説明の開始 ---
 @comment --- section $BL>$r@53N$K(B ---  @comment --- section 名を正確に ---
 @node gtt_ekn.setup,,, $BD64v2?4X?t(BE(k,n)  @node gtt_ekn.setup,,, 超幾何関数E(k,n)
 @subsection @code{gtt_ekn.setup}  @subsection @code{gtt_ekn.setup}
 @comment --- $B:w0zMQ%-!<%o!<%I(B  @comment --- 索引用キーワード
 @findex gtt_ekn.setup  @findex gtt_ekn.setup
   
 @table @t  @table @t
 @item gtt_ekn.setup()  @item gtt_ekn.setup()
 :: $BJ,;67W;;MQ$N4D6-@_Dj$r$*$3$J$&(B. $B8=:_$N4D6-$rJs9p$9$k(B.  :: 分散計算用の環境設定をおこなう. 現在の環境を報告する.
 @end table  @end table
   
 @comment --- $B0z?t$N4JC1$J@bL@(B ---  $B0J2<$^$@=q$$$F$J$$(B.  @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
 @table @var  @table @var
 @item return  @item return
   
 @end table  @end table
   
 @comment --- $B$3$3$G4X?t$N>\$7$$@bL@(B ---  @comment --- ここで関数の詳しい説明 ---
 @comment --- @itemize$B!A(B@end itemize $B$O2U>r=q$-(B ---  @comment --- @itemize〜@end itemize は箇条書き ---
 @comment --- @bullet $B$O9uE@IU$-(B ---  @comment --- @bullet は黒点付き ---
 @itemize @bullet  @itemize @bullet
 @item $B;HMQ$9$k%W%m%;%9$HAG?t$N8D?t(B, $B:G>.$NAG?t$rI=<($9$k(B. $B=`Hw$5$l$F$$$J$$>l9g$O$=$N;]$rI=<((B.  @item 使用するプロセスと素数の個数, 最小の素数を表示する. 準備されていない場合はその旨を表示.
 @item option nid (nid = Number of process ID)$B$rM?$($k$H;XDj$7$??t$@$1%W%m%;%9$rMQ0U$9$k(B.  @item このパッケージでの分散計算は複数のcpuを搭載した計算機で実行されることを想定している.
 @item option npl (npl = Prime List or Number of Prime List)$B$rM?$($k$H(Bnpl$B$,J8;zNs$N>l9g;XDj$5$l$?AG?t%j%9%H$N%U%!%$%k$rFI$_9~$`(B. npl$B$,<+A3?t$N>l9g$5$i$K(Boption minp (minp =MINimum Prime)$B$rM?$($k$H(Bminp$B$h$jBg$-$JAG?t$r(Bnpl$B8D@8@.$9$k(B. $B$=$N:](Boption fname (fname = File NAME)$B$rM?$($k$H@8@.$7$?AG?t%j%9%H$r(Bfname$B$H$7$FJ]B8$9$k(B.  @item option nps (または number_of_processes)を与えると指定した数だけプロセスを用意する.
   @item option nprm (または number_of_primes)を与えるとnprmが文字列の場合指定された素数リストのファイルを読み込む. nprmが自然数の場合さらにoption minp (minp =MINimum Prime)を与えるとminpより大きな素数をnprm個生成する. その際option fgp (または file_of_generated_primes)を与えると生成した素数リストをファイル名をfgpとして保存する.
   @item 上記のoption を指定しなかった場合次のデフォルト値が用いられる. nps=1. nprm=10. fgp=0.
   @item option report=1を与えると現在の環境の報告のみを行う. setup(|report=1)の別名としてreport関数を使用することもできる.
   @item option subprogs=[file1,file2,...] により分散計算の子供プロセスにロードすべきファイル file1, file2, ... を指定する. default は subprogs=["gtt_ekn/childprocess.rr"] である.
   @item gtt_ekn.set_debug_level(Mode) で Ekn_debug の値を設定する.
 @end itemize  @end itemize
   
 @comment --- @example$B!A(B@end example $B$O<B9TNc$NI=<((B ---  @comment --- @example〜@end example は実行例の表示 ---
 $BNc(B: $BAG?t$N%j%9%H$r@8@.$7$F%U%!%$%k(B p.txt $B$X=q$-=P$9(B.  例: 素数のリストを生成してファイル p.txt へ書き出す.
 @example  @example
 gtt_ekn.setup(|nps=2,nprm=20,minp=10^10,fgp="p.txt")$  gtt_ekn.setup(|nps=2,nprm=20,minp=10^10,fgp="p.txt")$
 @end example  @end example
   
   例: chinese remainder theorem (crt) を使って gmvector を計算.
   @example
   [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)$
   @end example
   
 @comment --- $B;2>H(B($B%j%s%/(B)$B$r=q$/(B ---  
   @comment --- 参照(リンク)を書く ---
 @table @t  @table @t
 @item $B;2>H(B  @item 参照
 @ref{gtt_ekn.nc}  @ref{gtt_ekn.nc}
 @ref{gtt_ekn.gmvector}  @ref{gtt_ekn.gmvector}
 @end table  @end table
   
 @comment --- ChangeLog $B$r=q$/(B. $B%=!<%9%3!<%I$N0LCV(B. $BJQ99F|;~(B $B$J$I(B CVS$B%5!<%P$r8+$k$?$a(B  @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
 @noindent  @noindent
 ChangeLog  ChangeLog
 @itemize @bullet  @itemize @bullet
 @item  @item
  $BJQ99$r<u$1$?%U%!%$%k$O(B   変更を受けたファイルは
  OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1,   OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1,
  gtt_ekn/g_mat_fac.rr   gtt_ekn/g_mat_fac.rr
   
 @end itemize  @end itemize
   
 @comment **********************************************************  @comment **********************************************************
 @comment --- $B"~"~"~"~(B  $B$N@bL@(B  @comment --- ◯◯◯◯  の説明
 @comment --- $B8D!9$N4X?t$N@bL@$N3+;O(B ---  @comment --- 個々の関数の説明の開始 ---
 @comment --- section $BL>$r@53N$K(B ---  @comment --- section 名を正確に ---
 @node gtt_ekn.upAlpha,,, $BD64v2?4X?t(BE(k,n)  @node gtt_ekn.upAlpha,,, 超幾何関数E(k,n)
 @subsection @code{gtt_ekn.upAlpha}  @node gtt_ekn.downAlpha,,, 超幾何関数E(k,n)
 @comment --- $B:w0zMQ%-!<%o!<%I(B  @subsection @code{gtt_ekn.upAlpha}, @code{gtt_ekn.downAlpha}
   @comment --- 索引用キーワード
 @findex gtt_ekn.upAlpha  @findex gtt_ekn.upAlpha
   @findex gtt_ekn.downAlpha
   
 @table @t  @table @t
 @item gtt_ekn.upAlpha(@var{i},@var{k},@var{n})  @item gtt_ekn.upAlpha(@var{i},@var{k},@var{n})
   @item gtt_ekn.downAlpha(@var{i},@var{k},@var{n})
 ::  ::
 @end table  @end table
   
 @comment --- $B0z?t$N4JC1$J@bL@(B ---  $B0J2<$^$@=q$$$F$J$$(B.  @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
 @table @var  @table @var
 @item i  a_i $B$r(B a_i+1 $B$HJQ2=$5$;$k(B contiguity relation.  @item i  a_i を a_i+1 (a_i を a_i-1) と変化させる contiguity relation.
 @item k  E(k+1,n+k+2)$B7?$ND64v2?4X?t$N(B k. $BJ,3dI=$G$O(B (k+1)$B!_(B(n+1).  @item k  E(k+1,n+k+2)型の超幾何関数の k. 分割表では (k+1)×(n+1).
 @item n  E(k+1,n+k+2)$B7?$ND64v2?4X?t$N(B n. $BJ,3dI=$G$O(B (k+1)$B!_(B(n+1).  @item n  E(k+1,n+k+2)型の超幾何関数の n. 分割表では (k+1)×(n+1).
 @item return  contiguity relation $B$N(B pfaffian_basis $B$K$D$$$F$N9TNsI=8=$rLa$9(B. [GM2016] $B$N(B Cor 6.3.  @item return  contiguity relation の pfaffian_basis についての行列表現を戻す. [GM2016] の Cor 6.3.
 @end table  @end table
   
 @comment --- $B$3$3$G4X?t$N>\$7$$@bL@(B ---  @comment --- ここで関数の詳しい説明 ---
 @comment --- @itemize$B!A(B@end itemize $B$O2U>r=q$-(B ---  @comment --- @itemize〜@end itemize は箇条書き ---
 @comment --- @bullet $B$O9uE@IU$-(B ---  @comment --- @bullet は黒点付き ---
 @itemize @bullet  @itemize @bullet
 @item  @item
  upAlpha $B$O!!(B[GM2016] $B$N(B Cor 6.3 $B$N9TNs(B U_i $B$rLa$9(B.   upAlpha は [GM2016] の Cor 6.3 の行列 U_i を戻す.
 @item $B4XO"$9$k3F4X?t$N4J7i$J@bL@$HNc$b2C$($k(B.  @item 関連する各関数の簡潔な説明と例も加える.
 @item a_i $B$r(B a_i-1 $B$HJQ2=$5$;$?$$>l9g$O4X?t(B downAlpha $B$rMQ$$$k(B.  @item a_i を a_i-1 と変化させたい場合は関数 downAlpha を用いる.
 @item a_i $B$HJ,3dI=$N<~JUOB$r8+$k$K$O(B, $B4X?t(B marginaltoAlpha([$B9TOB(B,$BNsOB(B]) $B$rMQ$$$k(B.  @item a_i と分割表の周辺和を見るには, 関数 marginaltoAlpha([行和,列和]) を用いる.
 @item  @item
    pfaffian_basis $B$O(B [GM2016] $B$N#4>O$N%Y%/%H%k(B F $B$KBP1~$9$kJPHyJ,$rLa$9(B.     pfaffian_basis は [GM2016] の4章のベクトル F に対応する偏微分を戻す.
   @item optional 引数 arule, xrule で a_i や x_i_j を数にしたものをより効率的に求めることができる. 変化をうけるパラメータを数にしてしまっても特にエラー表示はしない. a_0 で和の条件を調整しているので注意(Todo, double check).
 @end itemize  @end itemize
   
 @comment --- @example$B!A(B@end example $B$O<B9TNc$NI=<((B ---  @comment --- @example〜@end example は実行例の表示 ---
 $BNc(B: $B0J2<$NNc$O(B 2$B!_(B2$BJ,3dI=(B(E(2,4)), 2$B!_(B3$BJ,3dI=(B(E(2,5))$B$N>l9g$G$"$k(B.  例: 以下の例は 2×2分割表(E(2,4)), 2×3分割表(E(2,5))の場合である.
 [2225] $B$^$G$O=PNO$rN,$7$F$$$k(B.  [2225] までは出力を略している.
 @example  @example
 [2221] gtt_ekn.marginaltoAlpha([[1,4],[2,3]]);  [2221] gtt_ekn.marginaltoAlpha([[1,4],[2,3]]);
 [[a_0,-4],[a_1,-1],[a_2,3],[a_3,2]]  [[a_0,-4],[a_1,-1],[a_2,3],[a_3,2]]
 [2222] gtt_ekn.upAlpha(1,1,1);  // E(2,4) $B$N(B a_1 $BJ}8~$N(B  [2222] gtt_ekn.upAlpha(1,1,1);  // E(2,4) の a_1 方向の
                                 //     contiguity $B$rI=8=$9$k9TNs(B                                  //     contiguity を表現する行列
 [2223] gtt_ekn.upAlpha(2,1,1);  // E(2,4) $B$N(B a_2 $BJ}8~(B  [2223] gtt_ekn.upAlpha(2,1,1);  // E(2,4) の a_2 方向
 [2224] gtt_ekn.upAlpha(3,1,1);  // E(2,4) $B$N(B a_3 $BJ}8~(B  [2224] gtt_ekn.upAlpha(3,1,1);  // E(2,4) の a_3 方向
 [2225] function f(x_1_1);  [2225] function f(x_1_1);
 [2232] gtt_ekn.pfaffian_basis(f(x_1_1),1,1);  [2232] gtt_ekn.pfaffian_basis(f(x_1_1),1,1);
 [ f(x_1_1) ]  [ f(x_1_1) ]
 [ (f{1}(x_1_1)*x_1_1)/(a_2) ]  [ (f{1}(x_1_1)*x_1_1)/(a_2) ]
 [2233] function f(x_1_1,x_1_2);  [2233] function f(x_1_1,x_1_2);
 f() redefined.  f() redefined.
 [2234] gtt_ekn.pfaffian_basis(f(x_1_1,x_1_2),1,2); // E(2,5), 2*3 $BJ,3dI=(B  [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(x_1_1,x_1_2) ]
 [ (f{1,0}(x_1_1,x_1_2)*x_1_1)/(a_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) ]  [ (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 ]
   
 @end example  @end example
   
   
 @comment --- $B;2>H(B($B%j%s%/(B)$B$r=q$/(B ---  @comment --- 参照(リンク)を書く ---
 @table @t  @table @t
 @item $B;2>H(B  @item 参照
 @ref{gtt_ekn.nc}  @ref{gtt_ekn.nc}
 @ref{gtt_ekn.gmvector}  @ref{gtt_ekn.gmvector}
 @end table  @end table
   
 @comment --- ChangeLog $B$r=q$/(B. $B%=!<%9%3!<%I$N0LCV(B. $BJQ99F|;~(B $B$J$I(B CVS$B%5!<%P$r8+$k$?$a(B  @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
 @noindent  @noindent
 ChangeLog  ChangeLog
 @itemize @bullet  @itemize @bullet
 @item  @item
  $B$3$N4X?t$O(B [GM2016]   この関数は [GM2016]
 $B$GM?$($i$l$?%"%k%4%j%:%`$K=>$$(B contiguity relation $B$rF3=P$9$k(B.  で与えられたアルゴリズムに従い contiguity relation を導出する.
 @item  @item
  $BJQ99$r<u$1$?%U%!%$%k$O(B   変更を受けたファイルは
  OpenXM/src/asir-contrib/packages/src/gtt_ekn/ekn_pfaffian_8.rr 1.1.   OpenXM/src/asir-contrib/packages/src/gtt_ekn/ekn_pfaffian_8.rr 1.1.
 @end itemize  @end itemize
   
   
   @comment **********************************************************
   @comment --- ◯◯◯◯  の説明
   @comment --- 個々の関数の説明の開始 ---
   @comment --- section 名を正確に ---
   @node gtt_ekn.cmle,,, 超幾何関数E(k,n)
   @subsection @code{gtt_ekn.cmle}
   @comment --- 索引用キーワード
   @findex gtt_ekn.cmle
   
 @comment --- $B$*$^$8$J$$(B ---  @table @t
   @item gtt_ekn.cmle(@var{u}) u を観測データとするとき, P(U=u | row sum, column sum = these of U) を最大化する, 各セルの確率の近似値を求める.
   ::
   @end table
   
   @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
   @table @var
   @item u  観測データ(分割表)
   @item return  セルの確率(分割表形式)
   @end table
   
   @comment --- ここで関数の詳しい説明 ---
   @comment --- @itemize〜@end itemize は箇条書き ---
   @comment --- @bullet は黒点付き ---
   @itemize @bullet
   @item u を観測データとするとき, P(U=u | row sum, column sum = these of U) を最大化する, 各セルの確率の近似値を求める.
   @item optional parameter で algorithm の振る舞い(たとえば有理数を近似して, 分母分子が小さい有理数にする, gradient descent の step幅)を調整すべきだが, これは作業中. 2017.03.03
   @end itemize
   
   @comment --- @example〜@end example は実行例の表示 ---
   例: 2 x 4 分割表.
   @example
   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 ]]
   @end example
   
   例: 上の例は次の関数に.
   @example
   gtt_ekn.cmle_test3();
   @end example
   
   @comment --- 参照(リンク)を書く ---
   @table @t
   @item 参照
   @ref{gtt_ekn.expectation}
   @end table
   
   @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
   @noindent
   ChangeLog
   @itemize @bullet
   @item  gtt_ekn/mle.rr に本体がある.
   @item  gtt_ekn.rr の cmle 関数は wrapper.
   @end itemize
   @comment end cmle.
   
   @comment **********************************************************
   @comment --- ◯◯◯◯  の説明
   @comment --- 個々の関数の説明の開始 ---
   @comment --- section 名を正確に ---
   @node gtt_ekn.set_debug_level,,, 超幾何関数E(k,n)
   @node gtt_ekn.contiguity_mat_list_2,,, 超幾何関数E(k,n)
   @node gtt_ekn.show_path,,, 超幾何関数E(k,n)
   @node gtt_ekn.get_svalue,,, 超幾何関数E(k,n)
   @node gtt_ekn.assert1,,, 超幾何関数E(k,n)
   @node gtt_ekn.assert2,,, 超幾何関数E(k,n)
   @node gtt_ekn.prob1,,, 超幾何関数E(k,n)
   @subsection @code{gtt_ekn.set_debug_level}, @code{gtt_ekn.show_path}, @code{gtt_ekn.get_svalue}, @code{gtt_ekn.assert1}, @code{gtt_ekn.assert2}, @code{gtt_ekn.prob1}
   @comment --- 索引用キーワード
   @findex gtt_ekn.set_debug_level
   @findex gtt_ekn.contiguity_mat_list_2
   @findex gtt_ekn.show_path
   @findex gtt_ekn.get_svalue
   @findex gtt_ekn.assert1
   @findex gtt_ekn.assert2
   @findex gtt_ekn.prob1
   
   @table @t
   @item gtt_ekn.set_debug_level(@var{m}) debug メッセージのレベルを設定.
   @item gtt_ekn.contiguity_mat_list_2  使用する contiguity を構成.
   @item gtt_ekn.show_path()  どのように contiguity を適用したかの情報.
   @item gtt_ekn.get_svalue()  static 変数の値を得る.
   @item gtt_ekn.assert1(@var{N})  2x2 分割表で動作チェック.
   @item gtt_ekn.assert2(@var{N})  3x3 分割表で動作チェック.
   @item gtt_ekn.prob1(@var{R1},@var{R2},@var{Size})  R1 x R2 分割表用のテストデータを作る.
   ::
   @end table
   
   @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
   @table @var
   @item  m  レベル.
   @end table
   
   @comment --- ここで関数の詳しい説明 ---
   @comment --- @itemize〜@end itemize は箇条書き ---
   @comment --- @bullet は黒点付き ---
   @itemize @bullet
   @item (@var{m} & 0x1) == 0x1 の時 g_mat_fac_test_plain と g_mat_fac_itor の両方を呼び出し値を比較する (gtt_ekn.setup した状態で).
   @item (@var{m} & 0x2) == 0x2 の時 g_mat_fac_test への引数を tmp-input-数.ab として保存.
   @item (@var{m} & 0x4) == 0x4 の時 matrix factorial の計算の呼び出し引数を表示.
   @item @var{N} は問題の周辺和のサイズ.
   @item @code{get_svalue} の戻り値は @code{[Ekn_plist,Ekn_IDL,Ekn_debug,Ekn_mesg,XRule,ARule,Verbose,Ekn_Rq]} の値.
   @end itemize
   
   @comment --- @example〜@end example は実行例の表示 ---
   例.
   @example
   [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
   @end example
   
   例.
   @example
   [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]
   @end example
   [1 4 3 2] の index をもつパラメーター alpha の方向の contigity を求めそれを掛けて
   計算したことがわかる.  L[0] は用いた contiguity の行列.
   L[1] はcontiguity を適用する step 数.
   
   例. 値を計算せずに path のみ求めたい場合.
   @example
   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 ]
   @end example
   
   例. 値を計算せずに path のみ求めたい場合.
   gtt_ekn3 による新しいアルゴリズムによる path の表示.
   @example
   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]]
   @end example
   
   例. 0 が戻れば g_mat_fac_plain と指定した計算方法の結果が一致したことがわかる.
   option を書かないと g_mat_fac_int との比較となる.
   @example
   [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
   @end example
   なお二番目の例の timing (total) [例では省略] は mod 計算を subprocess がやっているので正しい値ではない. real time が計算時間の目安になる.
   
   例.
   @example
   [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 ]
   @end example
   
   @comment --- 参照(リンク)を書く ---
   @table @t
   @item 参照
   @ref{gtt_ekn.nc}
   @end table
   
   @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
   @noindent
   ChangeLog
   @itemize @bullet
   @item  gtt_ekn/ekn_eval.rr で matrix factorial の計算の呼び出し引数を表示する.
   @item grep 'iand(Ekn_debug,0x1)' *.rr でソースコードの該当の位置をさがす.
   @end itemize
   @comment end set_debug_level
   
   
   
   @node modular計算,,, 2元分割表HGMの関数
   @chapter modular計算
   
   @menu
   * gtt_ekn.chinese_itor::
   @end menu
   
   @node 中国剰余定理とitor,,, modular計算
   @section 中国剰余定理とitor
   
   @comment **********************************************************
   @comment --- ◯◯◯◯  の説明
   @comment --- 個々の関数の説明の開始 ---
   @comment --- section 名を正確に ---
   @node gtt_ekn.chinese_itor,,,
   @subsection @code{gtt_ekn.chinese_itor}
   @comment --- 索引用キーワード
   @findex gtt_ekn.chinese_itor 中国剰余定理とitor
   
   @table @t
   @item gtt_ekn.chinese_itor(@var{data},@var{idlist})
   :: mod p で計算した結果(ベクトル)から chinese remainder theorem, itor(integer to rational) で有理数ベクトルを得る.
   @end table
   
   @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
   @table @var
   @item return  [val, n]  ここで val は答え. また, n = n1*n2*...
   @item data   [[val1,n1],[val2,n2], ...], ここで val mod n1 = val1, val mod n2 = val2,...
   @item idlist  chinese, itor を実行するサーバIDのリスト.
   @end table
   
   @comment --- ここで関数の詳しい説明 ---
   @comment --- @itemize〜@end itemize は箇条書き ---
   @comment --- @bullet は黒点付き ---
   @itemize @bullet
   @item 中国剰余定理を用いて val0 mod n1 = val1, val0 mod n2 = val2, ... となる val0 を求める. val に algorithm itor を適用する.
   @item sqrt(n) より val0 が大きい時は itor が適用されて val0 が有理数 val=a/b に変換される. つまり b*x =1 mod n となる逆数 x を考えて, x*a % n = val0 となる数 val を戻す. 見つからないときは failure を戻す.
   @end itemize
   
   @comment --- @example〜@end example は実行例の表示 ---
   例: [3!, 5^3*3!]=[6,750] が戻り値.
   6 mod 109 =6, 750 mod 109=96 が最初の引数の [[6,96],109]. 以下同様.
   @example
   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]
   @end example
   
   
   @comment --- @example〜@end example は実行例の表示 ---
   例: gtt_ekn/childprocess.rr (server で実行される) の関数 chinese (chinese remainder theorem) と euclid.
   @example
   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
   @end example
   
   @comment --- @example〜@end example は実行例の表示 ---
   例: gtt_ekn/childprocess.rr (server で実行される) の関数 itor (integer to rational) の例.
   itor(Y,Q,Q2,Idx) では Y < Q2 なら Y がそのまま戻る.  Idx は 内部用の index で好きな数でよい. 戻り値の第2成分となる.
   @example
   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]]
   @end example
   
   
   @comment --- 参照(リンク)を書く ---
   @table @t
   @item 参照
   @ref{gtt_ekn.setup}
   @end table
   
   @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
   @noindent
   ChangeLog
   @itemize @bullet
   @item
    関連ファイルは
    gtt_ekn/g_mat_fac.rr
    gtt_ekn/childprocess.rr
   @end itemize
   
   @node binary splitting,,, 2元分割表HGMの関数
   @chapter binary splitting
   
   @menu
   * gtt_ekn3.init_dm_bsplit::
   * gtt_ekn3.setup_dm_bsplit::
   * gtt_ekn3.init_bsplit::
   @end menu
   
   @node matrix factorial,,, binary splitting
   @section matrix factorial
   
   @comment **********************************************************
   @comment --- ◯◯◯◯  の説明
   @comment --- 個々の関数の説明の開始 ---
   @comment --- section 名を正確に ---
   @node gtt_ekn3.init_bsplit,,,
   @node gtt_ekn3.init_dm_bsplit,,,
   @node gtt_ekn3.setup_dm_bsplit,,,
   @subsection @code{gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit}
   @comment --- 索引用キーワード
   @findex gtt_ekn3.init_dm_bsplit matrix factorial
   @findex gtt_ekn3.setup_dm_bsplit matrix factorial
   @findex gtt_ekn3.init_bsplit matrix factorial
   
   @table @t
   @item gtt_ekn3.init_bsplit(|minsize=16,levelmax=1);
   :: binary split の実行のためのパラメータを設定する.
   @item gtt_ekn3.init_dm_bsplit(|bsplit_x=0, bsplit_reduce=0)
   :: binary split の分散実行のためのパラメータを設定する.
   @item gtt_ekn3.setup_dm_bsplit(C)
   :: binary split の分散実行のために C 個のプロセスを立ち上げる.
   @end table
   
   @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
   @table @var
   @item C はlevelmax-1 に設定する. 特に levalmax=1 のときは分散計算を行わない.
   @item bsplit_x=1 のとき, debug 用に各プロセスを xterm で表示.
   @end table
   
   @comment --- ここで関数の詳しい説明 ---
   @comment --- @itemize〜@end itemize は箇条書き ---
   @comment --- @bullet は黒点付き ---
   @itemize @bullet
   @item expectation などの関数に bs=1 オプションを与えると matrix factorial を binary
   splitting method で計算する.
   @end itemize
   
   @comment --- @example〜@end example は実行例の表示 ---
   例: bs=1 と無い場合の比較.
   @example
   [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)
   @end example
   
   
   @comment --- @example〜@end example は実行例の表示 ---
   例: 分散計算する場合.
   分散計算はかえって遅くなる場合が多いので注意.
   下記の例での bsplit_x=1 option は
   debug windows を開くのでさらに遅くなる.
   gtt_ekn3.test_bs_dist(); でもテストできる.
   @example
   [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)$
   @end example
   
   @comment --- 参照(リンク)を書く ---
   @table @t
   @item 参照
   @ref{gtt_ekn3.gmvector}
   @ref{gtt_ekn3.expectation}
   @ref{gtt_ekn3.assert1}
   @ref{gtt_ekn3.assert2}
   @end table
   
   @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
   @noindent
   ChangeLog
   @itemize @bullet
   @item
    関連ファイルは
    gtt_ekn3/mfac_bs.rr
    gtt_ekn3/dm_bsplit.rr
   @end itemize
   
   
   
   @comment --- おまじない ---
 @node Index,,, Top  @node Index,,, Top
 @unnumbered Index  @unnumbered Index
 @printindex fn  @printindex fn
Line 764  ChangeLog
Line 1241  ChangeLog
 @summarycontents  @summarycontents
 @contents  @contents
 @bye  @bye
 @comment --- $B$*$^$8$J$$=*$j(B ---  @comment --- おまじない終り ---
   
   
 // 2 x m $BJ,3dI=$K$*$$$F;w$?5!G=$rM-$9$k4X?t$NMxMQNc$r;29M$^$G$K5-:\$9$k(B;  @comment テンプレート.  start_of_template.
 // $B@55,2=Dj?t$H$=$NHyJ,4XO"(B.  @comment **********************************************************
 // $B$=$N(B1.  @comment --- ◯◯◯◯  の説明
   @comment --- 個々の関数の説明の開始 ---
   @comment --- section 名を正確に ---
   @node gtt_ekn.hoge,,, 超幾何関数E(k,n)
   @subsection @code{gtt_ekn.hoge}
   @comment --- 索引用キーワード
   @findex gtt_ekn.hoge
   
   @table @t
   @item gtt_ekn.hoge(@var{i})
   ::
   @end table
   
   @comment --- 引数の簡単な説明 ---  以下まだ書いてない.
   @table @var
   @item i  hage
   @item return
   @end table
   
   @comment --- ここで関数の詳しい説明 ---
   @comment --- @itemize〜@end itemize は箇条書き ---
   @comment --- @bullet は黒点付き ---
   @itemize @bullet
   @item 説明.
   @end itemize
   
   @comment --- @example〜@end example は実行例の表示 ---
   例:
   @example
   [2221] gtt_ekn.hoge([[1,4],[2,3]]);
   @end example
   
   
   @comment --- 参照(リンク)を書く ---
   @table @t
   @item 参照
   @ref{gtt_ekn.nc}
   @ref{gtt_ekn.gmvector}
   @end table
   
   @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
   @noindent
   ChangeLog
   @itemize @bullet
   @item
   @end itemize
   @comment end_of_template
   
   
   // 2 x m 分割表において似た機能を有する関数の利用例を参考までに記載する;
   // 正規化定数とその微分関連.
   // その1.
 [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);  [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
 [-4,[-4,-3],-1]  [-4,[-4,-3],-1]
 [3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);  [3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
Line 777  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
Line 1305  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
 [ 1 1 1 ]  [ 1 1 1 ]
 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]  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]]]  [4483/124416,[[353/7776,1961/15552,185/1728],[553/20736,1261/15552,1001/13824]]]
 // $BLaCM$O(B [Z, [[d_11 Z, d_12 Z, d_13 Z],[d_21 Z, d_22 Z, d_23 Z]]] $B$NCM(B.  // 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],[d_21 Z, d_22 Z, d_23 Z]]] の値.
   
 // $B$=$N(B2.  // その2.
 [3079] tk_fd.log_ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);  [3079] 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 ]  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
 [ 1 1 1 ]  [ 1 1 1 ]
Line 787  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/
Line 1315  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/
 [-3.32333832422461674639485797719209322217260539267246045320,  [-3.32333832422461674639485797719209322217260539267246045320,
  [[1.25987062235110417131385233102832924380994869507026544724,3.49944233772027660049074280615659156814633058219942003122,2.97122462636627258532232879768012491635065804149007361142],   [[1.25987062235110417131385233102832924380994869507026544724,3.49944233772027660049074280615659156814633058219942003122,2.97122462636627258532232879768012491635065804149007361142],
   [0.740129377648895828686147668971670756190051304929734552754,2.25027883113986169975462859692170421592683470890028998438,2.00959179121124247155922373410662502788311398616997546285]]]    [0.740129377648895828686147668971670756190051304929734552754,2.25027883113986169975462859692170421592683470890028998438,2.00959179121124247155922373410662502788311398616997546285]]]
 // $BLaCM$O(B [log(Z),  // 戻値は [log(Z),
 //          [[d_11 log(Z), d_12 log(Z), d_13 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)]]]  //           [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]]
 // $B$N6a;wCM(B.  // の近似値.
   
 // $B$=$N(B3.  // その3.
 [3082] fd_hessian2(A[0],A[1],A[2],[1/2,1/3]);  [3082] 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 ]  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
 [4483/124416,[ 1961/15552 185/1728 ],  [4483/124416,[ 1961/15552 185/1728 ],
  [ 79/288 259/864 ]   [ 79/288 259/864 ]
  [ 259/864 47/288 ]]   [ 259/864 47/288 ]]
 // $BLaCM$O(B [F=F_D, gradient(F), Hessian(F)]  // 戻値は [F=F_D, gradient(F), Hessian(F)]
   
 // $B;29M(B.  // 参考.
 // ygahvec $B$G6R4X?tJ,$ND4@0(B. $BFHN)$7$?4X?t$O$J$$$h$&$@(B.  // ygahvec で巾関数分の調整. 独立した関数はないようだ.
   
 //-----------------------------------------------------------------------  //-----------------------------------------------------------------------
 // 2 x m $BJ,3dI=$K$*$$$F;w$?5!G=$rM-$9$k4X?t$NMxMQNc$r;29M$^$G$K5-:\$9$k(B;  // 2 x m 分割表において似た機能を有する関数の利用例を参考までに記載する;
 // $B4|BTCM4XO"(B.  // 期待値関連.
 [3079] A=tk_fd.marginal2abc([4,5],[2,4,3]);  [3079] A=tk_fd.marginal2abc([4,5],[2,4,3]);
 [-4,[-4,-3],-1]  [-4,[-4,-3],-1]
 [3080] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);  [3080] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
Line 814  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
Line 1342  RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]  Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
 [[5648/4483,7844/4483,4440/4483],  [[5648/4483,7844/4483,4440/4483],
  [3318/4483,10088/4483,9009/4483]]   [3318/4483,10088/4483,9009/4483]]
 // $B3F%;%k$N4|BTCM(B.  // 各セルの期待値.
   
 //-----------------------------------------------------------------------  //-----------------------------------------------------------------------
 // ot_hgm_ahg.rr $B$NNc(B.  $B<B83E*$J$?$a(B module $B2=$5$l$F$$$J$$(B.  // ot_hgm_ahg.rr の例.  実験的なため module 化されていない.
 [3237] import("ot_hgm_ahg.rr");  [3237] import("ot_hgm_ahg.rr");
 // 2 x 2 $BJ,3dI=(B.  // 2 x 2 分割表.
 [3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],  [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);          [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
 oohg_native=0, oohg_curl=1  oohg_native=0, oohg_curl=1
 [1376777025/625400597,1750225960/625400597,2375626557/625400597,3252978816/625400597]  [1376777025/625400597,1750225960/625400597,2375626557/625400597,3252978816/625400597]
 // 2 x 2 $BJ,3dI=$N4|BTCM(B.  // 2 x 2 分割表の期待値.
   
 // 2 x 3 $BJ,3dI=(B.  // 2 x 3 分割表.
 [3238] hgm_ahg_expected_values_contiguity(  [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]],   [[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);   [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]  [5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
 // 2 x 3 $BJ,3dI=$N4|BTCM(B. $B>e$HF1$8LdBj(B.  // 2 x 3 分割表の期待値. 上と同じ問題.
   
 /*  /*
   dojo, p.221.  $B@.@S(B3$B0J2<$N@8EL$O=8$a$F$R$H$D$K(B.    dojo, p.221.  成績3以下の生徒は集めてひとつに.
   2 1 1    2 1 1
   8 3 3    8 3 3
   0 2 6    0 2 6
   
   row sum: 4,14,8    row sum: 4,14,8
   column sum: 10,6,10    column sum: 10,6,10
   0 $B$r0l$D4^$`$N$G(B, (3,6) $B7?$N(B A $B$+$i(B 7 $BNsL\$rH4$/(B.    0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
 */  */
 // 3 x 3 $BJ,3dI=(B. $B9=B$E*(B0$B$,0l$D(B.  // 3 x 3 分割表. 構造的0が一つ.
   
 A=[[0,0,0,1,1,1, 0,0],  A=[[0,0,0,1,1,1, 0,0],
    [0,0,0,0,0,0, 1,1],     [0,0,0,0,0,0, 1,1],
Line 853  A=[[0,0,0,1,1,1, 0,0],
Line 1381  A=[[0,0,0,1,1,1, 0,0],
 B=[14,8,10,6,10];  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);  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);
   
 // $BEz(B.  // 答.
 [14449864949304/9556267369631,10262588586540/9556267369631,13512615942680/9556267369631,  [14449864949304/9556267369631,10262588586540/9556267369631,13512615942680/9556267369631,
  81112808747006/9556267369631,21816297744346/9556267369631,30858636683482/9556267369631,   81112808747006/9556267369631,21816297744346/9556267369631,30858636683482/9556267369631,
                               25258717886900/9556267369631,51191421070148/9556267369631]                                25258717886900/9556267369631,51191421070148/9556267369631]
   
   
 /*  /*
  $B>e$N%G!<%?$G(B 0 $B$r(B 1 $B$KJQ99(B.   上のデータで 0 を 1 に変更.
   2 1 1    2 1 1
   8 3 3    8 3 3
   1 2 6    1 2 6
Line 868  hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/
Line 1396  hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/
   row sum: 4,14,9    row sum: 4,14,9
   column sum: 11,6,10    column sum: 11,6,10
 */  */
 // 3 x 3 $BJ,3dI=(B.  // 3 x 3 分割表.
 A=[[0,0,0,1,1,1,0,0,0],  A=[[0,0,0,1,1,1,0,0,0],
    [0,0,0,0,0,0,1,1,1],     [0,0,0,0,0,0,1,1,1],
    [1,0,0,1,0,0,1,0,0],     [1,0,0,1,0,0,1,0,0],
Line 877  A=[[0,0,0,1,1,1,0,0,0],
Line 1405  A=[[0,0,0,1,1,1,0,0,0],
 B=[14,9,11,6,10];  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);  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);
   
 // $B4|BTCM(B, $BEz(B.  // 期待値, 答.
 [207017568232262040/147000422096729819,163140751505489940/147000422096729819,217843368649167296/147000422096729819,  [207017568232262040/147000422096729819,163140751505489940/147000422096729819,217843368649167296/147000422096729819,
  1185482401011137878/147000422096729819,358095302885438604/147000422096729819,514428205457640984/147000422096729819,   1185482401011137878/147000422096729819,358095302885438604/147000422096729819,514428205457640984/147000422096729819,
  224504673820628091/147000422096729819,360766478189450370/147000422096729819]   224504673820628091/147000422096729819,360766478189450370/147000422096729819]
   
 // Z $B$d$=$NHyJ,$N7W;;$O(B hgm_ahg_contiguity $B4X?t$,$*$3$J$&$,(B, $B$3$l$N4J0W%$%s%?!<%U%'!<%9$O(B  // Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
 // $B$^$@=q$$$F$J$$(B.  // まだ書いてない.
   
   
 4. x_ij $B$O(B [GM2016] $B$N#1>O$G(B,  4. x_ij は [GM2016] の1章で,
  $B$?$H$($P(B 3x3 $B$N;~(B [[1,1,1],[x_11,x_12,1],[x_21,x_22,1]]   たとえば 3x3 の時 [[1,1,1],[x_11,x_12,1],[x_21,x_22,1]]
 $B$H$J$C$F$$$k$,(B, [GM2016] $B$N(B Prop 7.1 $B$NBP1~$G$O(B,  となっているが, [GM2016] の Prop 7.1 の対応では,
    p = [[1,x_11,x_12],[1,x_21,x_22],[1,1,1]] $B$H$J$C$F$$$k$N$GCm0U(B.     p = [[1,x_11,x_12],[1,x_21,x_22],[1,1,1]] となっているので注意.

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.15

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