OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi - diff

Return to gtt_ekn-ja.texi CVS log

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.6 and 1.19

version 1.6, 2017/03/30 06:16:36

version 1.19, 2019/06/12 05:53:29

Line 1

%% $OpenXM: OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi,v 1.5 2017/03/03 12:02:32 takayama Exp $

%% $OpenXM: OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi,v 1.18 2019/04/09 00:48:14 takayama Exp $

%% ptex -kanji=utf8 gtt_ekn.texi (.texi までつける. platex でなく ptex)

%% xetex gtt_ekn-ja.texi (.texi までつける. )

%% 以下コメントは @comment で始める. \input texinfo 以降は普通の tex 命令は使えない.

\input texinfo

\input texinfo-ja

@iftex

@catcode`@#=6

@def@fref#1{@xrefX[#1,,@code{#1},,,]}

@def@b#1{{@bf@gt #1}}

@def@b#1{{@bf #1}}

@catcode`@#=@other

@end iftex

@overfullrule=0pt

@documentlanguage ja

@c -*-texinfo-*-

@comment %**start of header

@comment --- おまじない終り ---

Line 39

Line 40

@comment --- タイトル, バージョン, 著者名, 著作権表示 ---

@title 2元分割表HGM関数

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

@subtitle 1.1 版

@subtitle 3.0 版

@subtitle 2017 年 3 月 3 日

@subtitle 2019 年 6 月 7 日

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

@page

@vskip 0pt plus 1filll

@end titlepage

@comment --- おまじない ---

Line 65 Copyright @copyright{} Risa/Asir committers

Line 66 Copyright @copyright{} Risa/Asir committers

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

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

* modular計算

* Binary splitting

* Index::

@end menu

Line 77 Copyright @copyright{} Risa/Asir committers

Line 79 Copyright @copyright{} Risa/Asir committers

HGM(holonomic gradient method) を用いた2元分割表の関数について説明する.

ChangeLog の項目は www.openxm.org の cvsweb で

ソースコードを読む時の助けになる情報が書かれている.

このパッケージは下記のようにロードする.

@example

load("gtt_ekn3.rr");

@end example

gtt_ekn3.rr は gtt_ekn.rr を置き換える大きく改良されたパッケージである.

@noindent

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

@example

import("names.rr");

asir_contrib_update(|update=1);

@end example

@noindent

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

@itemize @bullet

@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,

@uref{http://arxiv.org/abs/1602.01637,arxiv:1602.01637 (version 1)}

@item [T2016]

Y.Tachibana, 差分ホロノミック勾配法のモジュラーメソッドによる計算の高速化,

2016, 神戸大学修士論文.

@item [GTT2016]

Y.Goto, Y.Tachibana, N.Takayama, 2元分割表に対する差分ホロノミック勾配法の実装,

数理研講究録(掲載予定).

数理研講究録.

@item [TGKT]

Y.Tachibana, Y.Goto, T.Koyama, N.Takayama,

Holonomic Gradient Method for Two Way Contingency Tables,

@uref{https://arxiv.org/abs/1803.04170, arxiv:1803.04170 (the 2nd version)}

@item [TKT2015]

N.Takayama, S.Kuriki, A.Takemura,

$A$-hypergeometric distributions and Newton polytopes.

arxiv:1510.02269

@uref{http://arxiv.org/abs/1510.02269,arxiv:1510.02269}

@end itemize

このマニュアルで説明する関数を用いたプログラム例は

gtt_ekn/test-t1.rr

gtt_ekn3/test-t1.rr

など.

Line 105 gtt_ekn/test-t1.rr

Line 123 gtt_ekn/test-t1.rr

@comment --- section ``実験的関数'' の subsection xyz_abc

@comment --- subsection xyz_pqr xyz_stu がある.

@menu

* gtt_ekn.gmvector::

* gtt_ekn3.gmvector::

* gtt_ekn.nc::

* gtt_ekn3.nc::

* gtt_ekn.lognc::

* gtt_ekn3.lognc::

* gtt_ekn.expectation::

* gtt_ekn3.expectation::

* gtt_ekn.setup::

* gtt_ekn3.setup::

* gtt_ekn.upAlpha::

* gtt_ekn3.upAlpha::

* gtt_ekn.cmle::

* gtt_ekn3.cmle::

* gtt_ekn3.set_debug_level::

* gtt_ekn3.contiguity_mat_list_2::

* gtt_ekn3.show_path::

* gtt_ekn3.get_svalue::

* gtt_ekn3.assert1::

* gtt_ekn3.assert2::

* gtt_ekn3.assert3::

* gtt_ekn3.prob1::

@end menu

@node 超幾何関数E(k,n),,, 2元分割表HGMの関数

Line 121 gtt_ekn/test-t1.rr

Line 147 gtt_ekn/test-t1.rr

@comment --- ◯◯◯◯ の説明

@comment --- 個々の関数の説明の開始 ---

@comment --- section 名を正確に ---

@node gtt_ekn.gmvector,,, 超幾何関数E(k,n)

@node gtt_ekn3.gmvector,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn.gmvector}

@subsection @code{gtt_ekn3.gmvector}

@comment --- 索引用キーワード

@findex gtt_ekn.gmvector

@findex gtt_ekn3.gmvector

@table @t

@item gtt_ekn.gmvector(@var{beta},@var{p})

@item gtt_ekn3.gmvector(@var{beta},@var{p})

:: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表に付随する超幾何関数

E(k,n) の値およびその微分の値を戻す.

@item gtt_ekn.ekn_cBasis_2(@var{beta},@var{p})

@item gtt_ekn3.ekn_cBasis_2(@var{beta},@var{p})

の別名である.

@end table

Line 183 b[2][i] は i 列目の列和である.

Line 209 b[2][i] は i 列目の列和である.

@verbatim

S=F_D(-b[1,1], [-b[2,2],...,-b[2,n+1]], b[2,1]-b[1,1]+1 ; y)/C,

@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])!

とおく.

1/C は L 字型の分割表

@verbatim

Line 205 gmvector は

Line 231 gmvector は

@item

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

@item

以下の option は expectation その他でも使える.

@item

option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう

[T2016].

分散計算用の各種パラメータの設定は

gtt_ekn.setup で行なう.

gtt_ekn3.setup で行なう.

@item

option bs=1. binary splitting method で matrix factorial を計算.

一般に 3x3 では効果あり(assert2(15|bs=1)), 5x5 (test5x5(20|bs=1))では遅くなる.

デフォールトは bs=0.

@item

option path. contiguity を適用する path をきめるアルゴリズムを指定.

path=2 (後藤, 松本の論文 [GM2016] の path). path=3 (論文 [TGKT] の path).

デフォールトは path=3.

@item

option interval. 通常の matrix factorial の計算では, 分母と分子をそれぞれ整数計算で計算し最後に約分をする. しかしながら数の中間膨張が一般的に発生しその中間膨張を解消するため

約分を一定間隔で行うと計算効率がよくなる.

interval に整数値を設定することにより行列による一次変換を interval 回するたびに約分を行う.

interval の最適値は問題毎に異なるためシステムがデフォールト値を設定することはない.

@item

option x=1. subprocess 毎に window を開く.

@end itemize

@comment --- @example〜@end example は実行例の表示 ---

例: 次は2 x 2 分割表で行和が [5,1], 列和が [3,3], 各セルの確率が

[[1/2,1/3],[1/7,1/5]] の場合の gmvector の値である.

@example

[3000] load("gtt_ekn.rr");

[3000] load("gtt_ekn3.rr");

[3001] ekn_gtt.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])

[3001] gtt_ekn3.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])

[775/27783]

[200/9261]

@end example

例: N を2以上の自然数とする時, Gauss の超幾何関数(この場合は多項式となる)

F(-36N,-11N,2N,(1-1/N)/56) の値は T3 に代入される ( [TGKT] ).

@comment ekn/Prog2/2x2.rr

@example

N=2;

T2=gtt_ekn3.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

例: interval option

@example

[4009] P=gtt_ekn3.prob1(5,5,100);

[[[100,200,300,400,500],[100,200,300,400,500]],[[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1/23,1/29,1/31,1/37],[1,1/41,1/43,1/47,1/53],[1,1,1,1,1]]]

[4010] util_timing(quote(gtt_ekn3.gmvector([[100,200,300,400,500],[100,200,300,400,500]], [[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1/23,1/29,1/31,1/37],[1,1/41,1/43,1/47,1/53],[1,1,1,1,1]])))[1];

[cpu,72.852,gc,0,memory,4462742364,real,72.856]

[4011] util_timing(quote(gtt_ekn3.gmvector([[100,200,300,400,500],[100,200,300,400,500]], [[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1/23,1/29,1/31,1/37],[1,1/41,1/43,1/47,1/53],[1,1,1,1,1]]|interval=100)))[1];

[cpu,67.484,gc,0,memory,3535280544,real,67.4844]

@end example

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

計算ができる.

守備範囲の異なるプログラム同士の比較, debug 用参考.

Line 277 def htest4() @{

Line 346 def htest4() @{

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.setup}

@ref{gtt_ekn3.setup}

@ref{gtt_ekn.pfaffian_basis}

@ref{gtt_ekn3.pfaffian_basis}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

Line 288 ChangeLog

Line 357 ChangeLog

@item

この関数は

[GM2016] のアルゴリズムおよび

[T2016] による modular method を用いた高速化を実装したものである.

[T2016] による modular method を用いた高速化,

[TGKT] の高速化を実装したものである.

@item

変更を受けたファイルは

OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_pfaffian_8.rr

@item

interval option について変更を受けたファイルは

OpenXM/src/asir-contrib/packages/src/gtt_ekn3/ekn_eval.rr 1.6

@end itemize

@comment **********************************************************

@node gtt_ekn.nc,,, 超幾何関数E(k,n)

@node gtt_ekn3.nc,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn.nc}

@subsection @code{gtt_ekn3.nc}

@comment --- 索引用キーワード

@findex gtt_ekn.nc

@findex gtt_ekn3.nc

@table @t

@item gtt_ekn.nc(@var{beta},@var{p})

@item gtt_ekn3.nc(@var{beta},@var{p})

:: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z

およびその微分の値を戻す.

@end table

Line 335 nc は gmvector の値を元に, [GM2016] の Prop

Line 408 nc は gmvector の値を元に, [GM2016] の Prop

@item

option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.

分散計算用の各種パラメータの設定は

gtt_ekn.setup で行なう.

gtt_ekn3.setup で行なう.

その他の option は gmvector を参照.

@end itemize

@comment --- @example〜@end example は実行例の表示 ---

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

@example

[2237] gtt_ekn.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);

[2237] gtt_ekn3.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);

[4483/124416,[ 353/7776 1961/15552 185/1728 ]

[ 553/20736 1261/15552 1001/13824 ]]

@end example

Line 366 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/

Line 440 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.setup}

@ref{gtt_ekn3.setup}

@ref{gtt_ekn.lognc}

@ref{gtt_ekn3.lognc}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

Line 381 ChangeLog

Line 455 ChangeLog

@comment **********************************************************

@node gtt_ekn.lognc,,, 超幾何関数E(k,n)

@node gtt_ekn3.lognc,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn.lognc}

@subsection @code{gtt_ekn3.lognc}

@comment --- 索引用キーワード

@findex gtt_ekn.lognc

@findex gtt_ekn3.lognc

@table @t

@item gtt_ekn.lognc(@var{beta},@var{p})

@item gtt_ekn3.lognc(@var{beta},@var{p})

:: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z

の log の近似値およびその微分の近似値を戻す.

@end table

Line 410 ChangeLog

Line 484 ChangeLog

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

@item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.

分散計算用の各種パラメータの設定は

gtt_ekn.setup で行なう.

gtt_ekn3.setup で行なう.

@end itemize

@comment --- @example〜@end example は実行例の表示 ---

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

@example

[2238] gtt_ekn.lognc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);

[2238] gtt_ekn3.lognc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);

[-3.32333832422461674630,[ 5648/4483 15688/4483 13320/4483 ]

[ 3318/4483 10088/4483 9009/4483 ]]

@end example

Line 443 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/

Line 517 Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.setup}

@ref{gtt_ekn3.setup}

@ref{gtt_ekn.nc}

@ref{gtt_ekn3.nc}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

Line 457 ChangeLog

Line 531 ChangeLog

@end itemize

@comment **********************************************************

@node gtt_ekn.expectation,,, 超幾何関数E(k,n)

@node gtt_ekn3.expectation,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn.expectation}

@subsection @code{gtt_ekn3.expectation}

@comment --- 索引用キーワード

@findex gtt_ekn.expectation

@findex gtt_ekn3.expectation

@table @t

@item gtt_ekn.expectation(@var{beta},@var{p})

@item gtt_ekn3.expectation(@var{beta},@var{p})

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

@end table

Line 482 ChangeLog

Line 556 ChangeLog

@comment --- @bullet は黒点付き ---

@itemize @bullet

@item

[GM2016] の Algorithm 7.8 の実装.

[GM2016] の Algorithm 7.8 の実装. [TGKT] による高速化版 (path=3) がデフォールト.

@item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.

分散計算用の各種パラメータの設定は

gtt_ekn.setup で行なう.

gtt_ekn3.setup で行なう.

@item option index を与えると, 指定された成分の期待値のみ計算する.

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

@item その他の option は gmvector を参照.

@end itemize

@comment --- @example〜@end example は実行例の表示 ---

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

@example

[2235] gtt_ekn.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);

[2235] gtt_ekn3.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);

[ 2/3 1/3 ]

[ 4/3 8/3 ]

[2236] gtt_ekn.expectation([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);

[2236] gtt_ekn3.expectation([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);

[ 5648/4483 7844/4483 4440/4483 ]

[ 3318/4483 10088/4483 9009/4483 ]

[2442] gtt_ekn.expectation([[4,14,9],[11,6,10]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]]);

[2442] gtt_ekn3.expectation([[4,14,9],[11,6,10]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]]);

[ 207017568232262040/147000422096729819 163140751505489940/147000422096729819

217843368649167296/147000422096729819 ]

[ 1185482401011137878/147000422096729819 358095302885438604/147000422096729819

Line 613 hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/

Line 688 hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.setup}

@ref{gtt_ekn3.setup}

@ref{gtt_ekn.nc}

@ref{gtt_ekn3.nc}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

Line 631 ChangeLog

Line 706 ChangeLog

@comment --- ◯◯◯◯ の説明

@comment --- 個々の関数の説明の開始 ---

@comment --- section 名を正確に ---

@node gtt_ekn.setup,,, 超幾何関数E(k,n)

@node gtt_ekn3.setup,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn.setup}

@subsection @code{gtt_ekn3.setup}

@comment --- 索引用キーワード

@findex gtt_ekn.setup

@findex gtt_ekn3.setup

@table @t

@item gtt_ekn.setup()

@item gtt_ekn3.setup()

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

@end table

Line 657 ChangeLog

Line 732 ChangeLog

@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 option subprogs=[file1,file2,...] により分散計算の子供プロセスにロードすべきファイル file1, file2, ... を指定する. default は subprogs=["gtt_ekn3/childprocess.rr"] である.

@item gtt_ekn.set_debug_mode(Mode) で Ekn_debug の値を設定する.

@item gtt_ekn3.set_debug_level(Mode) で Ekn_debug の値を設定する.

@end itemize

@comment --- @example〜@end example は実行例の表示 ---

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

@example

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

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

@end example

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

@example

[2867] gtt_ekn3.setup(|nprm=20,minp=10^20);

[2868] N=2; T2=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],

[[1,(1-1/N)/56],[1,1]] | crt=1)$

@end example

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.nc}

@ref{gtt_ekn3.nc}

@ref{gtt_ekn.gmvector}

@ref{gtt_ekn3.gmvector}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

Line 690 ChangeLog

Line 772 ChangeLog

@comment --- ◯◯◯◯ の説明

@comment --- 個々の関数の説明の開始 ---

@comment --- section 名を正確に ---

@node gtt_ekn.upAlpha,,, 超幾何関数E(k,n)

@node gtt_ekn3.upAlpha,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn.upAlpha}

@node gtt_ekn3.downAlpha,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn3.upAlpha}, @code{gtt_ekn3.downAlpha}

@comment --- 索引用キーワード

@findex gtt_ekn.upAlpha

@findex gtt_ekn3.upAlpha

@findex gtt_ekn3.downAlpha

@table @t

@item gtt_ekn.upAlpha(@var{i},@var{k},@var{n})

@item gtt_ekn3.upAlpha(@var{i},@var{k},@var{n})

@item gtt_ekn3.downAlpha(@var{i},@var{k},@var{n})

@end table

@comment --- 引数の簡単な説明 --- 以下まだ書いてない.

@table @var

@item i a_i を a_i+1 と変化させる contiguity relation.

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

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

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

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

Line 719 ChangeLog

Line 804 ChangeLog

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

@item

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

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

@end itemize

@comment --- @example〜@end example は実行例の表示 ---

例: 以下の例は 2×2分割表(E(2,4)), 2×3分割表(E(2,5))の場合である.

[2225] までは出力を略している.

@example

[2221] gtt_ekn.marginaltoAlpha([[1,4],[2,3]]);

[2221] gtt_ekn3.marginaltoAlpha([[1,4],[2,3]]);

[[a_0,-4],[a_1,-1],[a_2,3],[a_3,2]]

[2222] gtt_ekn.upAlpha(1,1,1); // E(2,4) の a_1 方向の

[2222] gtt_ekn3.upAlpha(1,1,1); // E(2,4) の a_1 方向の

// contiguity を表現する行列

[2223] gtt_ekn.upAlpha(2,1,1); // E(2,4) の a_2 方向

[2223] gtt_ekn3.upAlpha(2,1,1); // E(2,4) の a_2 方向

[2224] gtt_ekn.upAlpha(3,1,1); // E(2,4) の a_3 方向

[2224] gtt_ekn3.upAlpha(3,1,1); // E(2,4) の a_3 方向

[2225] function f(x_1_1);

[2232] gtt_ekn.pfaffian_basis(f(x_1_1),1,1);

[2232] gtt_ekn3.pfaffian_basis(f(x_1_1),1,1);

[ f(x_1_1) ]

[ (f{1}(x_1_1)*x_1_1)/(a_2) ]

[2233] function f(x_1_1,x_1_2);

f() redefined.

[2234] gtt_ekn.pfaffian_basis(f(x_1_1,x_1_2),1,2); // E(2,5), 2*3 分割表

[2234] gtt_ekn3.pfaffian_basis(f(x_1_1,x_1_2),1,2); // E(2,5), 2*3 分割表

[ f(x_1_1,x_1_2) ]

[ (f{1,0}(x_1_1,x_1_2)*x_1_1)/(a_2) ]

[ (f{0,1}(x_1_1,x_1_2)*x_1_2)/(a_3) ]

[2235] RuleA=[[a_2,1/3],[a_3,1/2]]$ RuleX=[[x_1_1,1/5]]$

base_replace(gtt_ekn3.upAlpha(1,1,1),append(RuleA,RuleX))

-gtt_ekn3.upAlpha(1,1,1 | arule=RuleA, xrule=RuleX);

[ 0 0 ]

@end example

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.nc}

@ref{gtt_ekn3.nc}

@ref{gtt_ekn.gmvector}

@ref{gtt_ekn3.gmvector}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

Line 768 ChangeLog

Line 862 ChangeLog

@comment --- ◯◯◯◯ の説明

@comment --- 個々の関数の説明の開始 ---

@comment --- section 名を正確に ---

@node gtt_ekn.cmle,,, 超幾何関数E(k,n)

@node gtt_ekn3.cmle,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn.cmle}

@subsection @code{gtt_ekn3.cmle}

@comment --- 索引用キーワード

@findex gtt_ekn.cmle

@findex gtt_ekn3.cmle

@table @t

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

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

@end table

Line 796 ChangeLog

Line 890 ChangeLog

例: 2 x 4 分割表.

@example

U=[[1,1,2,3],[1,3,1,1]];

gtt_ekn.cmle(U);

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

Line 805 gtt_ekn.cmle(U);

Line 899 gtt_ekn.cmle(U);

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

@example

gtt_ekn.cmle_test3();

gtt_ekn3.cmle_test3();

@end example

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.expectation}

@ref{gtt_ekn3.expectation}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

@noindent

ChangeLog

@itemize @bullet

@item gtt_ekn/mle.rr に本体がある.

@item gtt_ekn3/mle.rr に本体がある.

@item gtt_ekn.rr の cmle 関数は wrapper.

@item gtt_ekn3.rr の cmle 関数は wrapper.

@end itemize

@comment end cmle.

@comment **********************************************************

@comment --- ◯◯◯◯ の説明

@comment --- 個々の関数の説明の開始 ---

@comment --- section 名を正確に ---

@node gtt_ekn3.set_debug_level,,, 超幾何関数E(k,n)

@node gtt_ekn3.contiguity_mat_list_2,,, 超幾何関数E(k,n)

@node gtt_ekn3.show_path,,, 超幾何関数E(k,n)

@node gtt_ekn3.get_svalue,,, 超幾何関数E(k,n)

@node gtt_ekn3.assert1,,, 超幾何関数E(k,n)

@node gtt_ekn3.assert2,,, 超幾何関数E(k,n)

@node gtt_ekn3.assert3,,, 超幾何関数E(k,n)

@node gtt_ekn3.prob1,,, 超幾何関数E(k,n)

@subsection @code{gtt_ekn3.set_debug_level}, @code{gtt_ekn3.show_path}, @code{gtt_ekn3.get_svalue}, @code{gtt_ekn3.assert1}, @code{gtt_ekn3.assert2}, @code{gtt_ekn3.assert3}, @code{gtt_ekn3.prob1}

@comment --- 索引用キーワード

@findex gtt_ekn3.set_debug_level

@findex gtt_ekn3.contiguity_mat_list_2

@findex gtt_ekn3.show_path

@findex gtt_ekn3.get_svalue

@findex gtt_ekn3.assert1

@findex gtt_ekn3.assert2

@findex gtt_ekn3.assert3

@findex gtt_ekn3.prob1

@table @t

@item gtt_ekn3.set_debug_level(@var{m}) debug メッセージのレベルを設定.

@item gtt_ekn3.contiguity_mat_list_2 使用する contiguity を構成.

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

@item gtt_ekn3.get_svalue() static 変数の値を得る.

@item gtt_ekn3.assert1(@var{N}) 2x2 分割表で動作チェック.

@item gtt_ekn3.assert2(@var{N}) 3x3 分割表で動作チェック.

@item gtt_ekn3.assert3(@var{R1}, @var{R2}, @var{Size}) R1 x R2 分割表で並列動作の動作チェック.

@item gtt_ekn3.prob1(@var{R1},@var{R2},@var{Size}) R1 x R2 分割表用のテストデータを作る.

@end table

@node modular計算,,, 2元分割表HGMの関数

@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_ekn3.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]} の値.

@item assert3 の options: x=1, subprocess の window を表示. nps=m, m 個のプロセスで contiguity を求める (contiguity_mat_list_3). crt, interval などは gmvector などと共通の

option. timing data を表示するには load("gtt_ekn3/ekn_eval-timing.rr"); しておく.

@end itemize

@comment --- @example〜@end example は実行例の表示 ---

例.

@example

[2846] gtt_ekn3.set_debug_level(0x4);

[2847] N=2; T2=gtt_ekn3.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_ekn3.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_ekn3.nc([[4,5,6],[2,4,9]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]])$

[2660] L=matrix_transpose(gtt_ekn3.show_path())$

[2661] L[2];

[2 1]

@end example

[2 1] の index をもつパラメーター alpha の方向の contigity を求めそれを掛けて

計算したことがわかる. L[0] は用いた contiguity の行列.

L[1] はcontiguity を適用する step 数.

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

@example

A=gtt_ekn3.marginaltoAlpha_list([[400,410,1011],[910,411,500]])$

[2666] gtt_ekn3.contiguity_mat_list_2(A,2,2)$

[2667] L=matrix_transpose(gtt_ekn3.show_path())$

[2668] L[2];

[ 2 1 5 4 3 ]

[2669] gtt_ekn3.contiguity_mat_list_3(A,2,2)$ // new alg in [TGKT]

[2670] L=matrix_transpose(gtt_ekn3.show_path())$

[2671] L[2];

[2 1] // shorter

@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_ekn3.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_ekn3.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_ekn3.setup(); to use the crt option.

Timing (rational) =4.51349 (CPU) + 6.32174 (GC) = 10.8352 (total)

diff of both method =

[ 0 0 0 ]

[8860]

[8863] gtt_ekn3.setup(|nprm=100,minp=10^50);

Number of processes = 1.

Number of primes = 100.

Min of plist = 100000000000000000000000000000000000000000000000151.

[8864] gtt_ekn3.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

3x5 分割表. 周辺和は 10 に比例する一定の数(factor option も関係. ソースを参照).

cell 確率は1/素数で生成される.

@comment grep testnxn ekn/Prog2/*.rr ; grep test_nxn ekn/Prog2/*.rr も見よ.

[9054] L=gtt_ekn3.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_ekn3.expectation(L[0],L[1]));

[ 1.65224223218613 ... snip ]

@end example

例:

@example

[5779] import("gtt_ekn3.rr"); load("gtt_ekn3/ekn_eval-timing.rr");

[5780] gtt_ekn3.assert3(5,5,100 | nps=32, interval=100);

-- snip

Parallel method: Number of process=32, File name tmp-gtt_ekn3/p300.txt is written.

Number of processes = 32.

-- snip

initialPoly of path=3: [ 2.184 0 124341044 2.1831 ] [CPU(s),0,*,real(s)]

contiguity_mat_list_3 of path=3: [ 0.04 0 630644 9.6774 ] [CPU(s),0,*,real(s)]

Note: interval option will lead faster evaluation. We do not use g_mat_fac_itor (crt). Call gtt_ekn3.setup(); to use the crt option.

g_mat_fac of path=3: [ 21.644 0 1863290168 21.6457 ] [CPU(s),0,*,real(s)]

Done. Saved in 2.ab

Diff (should be 0)=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,..., 0,0,0]

@end example

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn3.nc}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

@noindent

ChangeLog

@itemize @bullet

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

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

@end itemize

@comment end set_debug_level

@node modular計算,,, Top

@chapter modular計算

@menu

* gtt_ekn.chinese_itor::

* gtt_ekn3.chinese_itor::

@end menu

@node 中国剰余定理とitor,,, modular計算

Line 839 ChangeLog

Line 1107 ChangeLog

@comment --- ◯◯◯◯ の説明

@comment --- 個々の関数の説明の開始 ---

@comment --- section 名を正確に ---

@node gtt_ekn.chinese_itor,,,

@node gtt_ekn3.chinese_itor,,,

@subsection @code{gtt_ekn.chinese_itor}

@subsection @code{gtt_ekn3.chinese_itor}

@comment --- 索引用キーワード

@findex gtt_ekn.chinese_itor 中国剰余定理とitor

@findex gtt_ekn3.chinese_itor 中国剰余定理とitor

@table @t

@item gtt_ekn.chinese_itor(@var{data},@var{idlist})

@item gtt_ekn3.chinese_itor(@var{data},@var{idlist})

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

@end table

Line 868 ChangeLog

Line 1136 ChangeLog

例: [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");

gtt_ekn3.setup(|nps=2,nprm=3,minp=101,fgp="p_small.txt");

SS=gtt_ekn.get_svalue();

SS=gtt_ekn3.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]);

gtt_ekn3.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]);

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

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

@example

load("gtt_ekn/childprocess.rr");

load("gtt_ekn3/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)

// で求めると,

Line 898 euclid(3,103); // mod 103 での 3 の逆数. つま�

Line 1166 euclid(3,103); // mod 103 での 3 の逆数. つま�

@end example

@comment --- @example〜@end example は実行例の表示 ---

例: gtt_ekn/childprocess.rr (server で実行される) の関数 itor (integer to rational) の例.

例: gtt_ekn3/childprocess.rr (server で実行される) の関数 itor (integer to rational) の例.

itor(Y,Q,Q2,Idx) では Y < Q2 なら Y がそのまま戻る. Idx は内部用の index で好きな数でよい. 戻り値の第２成分となる.

@example

load("gtt_ekn/childprocess.rr");

load("gtt_ekn3/childprocess.rr");

for (I=1;I<11; I++) print([I,itor(I,11,3,0)]);

[1,[1,0]]

[2,[2,0]]

Line 919 for (I=1;I<11; I++) print([I,itor(I,11,3,0)]);

Line 1187 for (I=1;I<11; I++) print([I,itor(I,11,3,0)]);

@comment --- 参照(リンク)を書く ---

@table @t

@item 参照

@ref{gtt_ekn.setup}

@ref{gtt_ekn3.setup}

@end table

@comment --- ChangeLog を書く. ソースコードの位置. 変更日時など CVSサーバを見るため

Line 928 ChangeLog

Line 1196 ChangeLog

@itemize @bullet

@item