Annotation of OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi, Revision 1.20
1.20 ! takayama 1: %% $OpenXM: OpenXM/src/asir-contrib/packages/doc/gtt_ekn/gtt_ekn-ja.texi,v 1.19 2019/06/12 05:53:29 takayama Exp $
1.12 takayama 2: %% xetex gtt_ekn-ja.texi (.texi までつける. )
1.6 takayama 3: %% 以下コメントは @comment で始める. \input texinfo 以降は普通の tex 命令は使えない.
1.7 takayama 4: \input texinfo-ja
1.1 takayama 5: @iftex
6: @catcode`@#=6
7: @def@fref#1{@xrefX[#1,,@code{#1},,,]}
1.7 takayama 8: @def@b#1{{@bf #1}}
1.1 takayama 9: @catcode`@#=@other
10: @end iftex
11: @overfullrule=0pt
1.7 takayama 12: @documentlanguage ja
1.1 takayama 13: @c -*-texinfo-*-
14: @comment %**start of header
1.6 takayama 15: @comment --- おまじない終り ---
1.1 takayama 16:
1.6 takayama 17: @comment --- GNU info ファイルの名前 ---
1.20 ! takayama 18: @setfilename asir-contrib-gtt_ekn
1.1 takayama 19:
1.6 takayama 20: @comment --- タイトル ---
21: @settitle 2元分割表HGM
1.1 takayama 22:
23: @comment %**end of header
24: @comment %@setchapternewpage odd
25:
1.6 takayama 26: @comment --- おまじない ---
1.1 takayama 27: @ifinfo
28: @macro fref{name}
29: @ref{\name\,,@code{\name\}}
30: @end macro
31: @end ifinfo
32:
33: @iftex
34: @comment @finalout
35: @end iftex
36:
37: @titlepage
1.6 takayama 38: @comment --- おまじない終り ---
1.1 takayama 39:
1.6 takayama 40: @comment --- タイトル, バージョン, 著者名, 著作権表示 ---
41: @title 2元分割表HGM関数
42: @subtitle Risa/Asir 2元分割表HGM関数説明書
1.19 takayama 43: @subtitle 3.0 版
44: @subtitle 2019 年 6 月 7 日
1.1 takayama 45:
46: @author by Y.Goto, Y.Tachibana, N.Takayama
47: @page
48: @vskip 0pt plus 1filll
49: Copyright @copyright{} Risa/Asir committers
1.19 takayama 50: 2004--2019. All rights reserved.
1.1 takayama 51: @end titlepage
52:
1.6 takayama 53: @comment --- おまじない ---
1.1 takayama 54: @synindex vr fn
1.6 takayama 55: @comment --- おまじない終り ---
1.1 takayama 56:
1.6 takayama 57: @comment --- @node は GNU info, HTML 用 ---
58: @comment --- @node の引数は node-name, next, previous, up ---
1.1 takayama 59: @node Top,, (dir), (dir)
60:
1.6 takayama 61: @comment --- @menu は GNU info, HTML 用 ---
62: @comment --- chapter 名を正確に並べる ---
63: @comment --- この文書では chapter XYZ, Chapter Index がある.
64: @comment --- Chapter XYZ には section XYZについて, section XYZに関する関数がある.
1.1 takayama 65: @menu
1.6 takayama 66: * 2元分割表HGMの関数説明書について::
67: * 2元分割表HGMの関数::
68: * modular計算
1.19 takayama 69: * Binary splitting
1.1 takayama 70: * Index::
71: @end menu
72:
1.6 takayama 73: @comment --- chapter の開始 ---
74: @comment --- 親 chapter 名を正確に. 親がない場合は Top ---
75: @node 2元分割表HGMの関数説明書について,,, Top
76: @chapter 2元分割表HGMの関数説明書について
77:
78: この説明書では
79: HGM(holonomic gradient method) を用いた2元分割表の関数について説明する.
80: ChangeLog の項目は www.openxm.org の cvsweb で
81: ソースコードを読む時の助けになる情報が書かれている.
1.8 takayama 82: このパッケージは下記のようにロードする.
83: @example
1.13 takayama 84: load("gtt_ekn3.rr");
1.8 takayama 85: @end example
1.13 takayama 86: gtt_ekn3.rr は gtt_ekn.rr を置き換える大きく改良されたパッケージである.
1.8 takayama 87: @noindent
88: 最新版の asir-contrib package を取得するには, 下記のように更新関数を呼び出す.
89: @example
90: import("names.rr");
91: asir_contrib_update(|update=1);
92: @end example
93: @noindent
1.6 takayama 94: 本文中で引用している文献を列挙する.
1.1 takayama 95: @itemize @bullet
96: @item [GM2016]
1.19 takayama 97: Y.Goto, K.Matsumoto, Pfaffian equations and contiguity relations of the hypergeometric function of type (k+1,k+n+2) and their applications,
98: @uref{http://arxiv.org/abs/1602.01637,arxiv:1602.01637 (version 1)}
1.1 takayama 99: @item [T2016]
1.6 takayama 100: Y.Tachibana, 差分ホロノミック勾配法のモジュラーメソッドによる計算の高速化,
101: 2016, 神戸大学修士論文.
1.1 takayama 102: @item [GTT2016]
1.6 takayama 103: Y.Goto, Y.Tachibana, N.Takayama, 2元分割表に対する差分ホロノミック勾配法の実装,
1.8 takayama 104: 数理研講究録.
105: @item [TGKT]
106: Y.Tachibana, Y.Goto, T.Koyama, N.Takayama,
107: Holonomic Gradient Method for Two Way Contingency Tables,
1.19 takayama 108: @uref{https://arxiv.org/abs/1803.04170, arxiv:1803.04170 (the 2nd version)}
1.1 takayama 109: @item [TKT2015]
110: N.Takayama, S.Kuriki, A.Takemura,
111: $A$-hypergeometric distributions and Newton polytopes.
1.19 takayama 112: @uref{http://arxiv.org/abs/1510.02269,arxiv:1510.02269}
1.1 takayama 113: @end itemize
114:
1.6 takayama 115: このマニュアルで説明する関数を用いたプログラム例は
1.19 takayama 116: gtt_ekn3/test-t1.rr
1.6 takayama 117: など.
1.1 takayama 118:
1.4 takayama 119:
1.6 takayama 120: @node 2元分割表HGMの関数,,, Top
121: @chapter 2元分割表HGMの関数
1.1 takayama 122:
1.6 takayama 123: @comment --- section ``実験的関数'' の subsection xyz_abc
124: @comment --- subsection xyz_pqr xyz_stu がある.
1.1 takayama 125: @menu
1.19 takayama 126: * gtt_ekn3.gmvector::
127: * gtt_ekn3.nc::
128: * gtt_ekn3.lognc::
129: * gtt_ekn3.expectation::
130: * gtt_ekn3.setup::
131: * gtt_ekn3.upAlpha::
132: * gtt_ekn3.cmle::
133: * gtt_ekn3.set_debug_level::
134: * gtt_ekn3.contiguity_mat_list_2::
135: * gtt_ekn3.show_path::
136: * gtt_ekn3.get_svalue::
137: * gtt_ekn3.assert1::
138: * gtt_ekn3.assert2::
139: * gtt_ekn3.assert3::
140: * gtt_ekn3.prob1::
1.1 takayama 141: @end menu
142:
1.6 takayama 143: @node 超幾何関数E(k,n),,, 2元分割表HGMの関数
144: @section 超幾何関数E(k,n)
1.1 takayama 145:
146: @comment **********************************************************
1.6 takayama 147: @comment --- ◯◯◯◯ の説明
148: @comment --- 個々の関数の説明の開始 ---
149: @comment --- section 名を正確に ---
1.19 takayama 150: @node gtt_ekn3.gmvector,,, 超幾何関数E(k,n)
151: @subsection @code{gtt_ekn3.gmvector}
1.6 takayama 152: @comment --- 索引用キーワード
1.19 takayama 153: @findex gtt_ekn3.gmvector
1.1 takayama 154:
155: @table @t
1.19 takayama 156: @item gtt_ekn3.gmvector(@var{beta},@var{p})
1.6 takayama 157: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表に付随する超幾何関数
158: E(k,n) の値およびその微分の値を戻す.
1.19 takayama 159: @item gtt_ekn3.ekn_cBasis_2(@var{beta},@var{p})
1.6 takayama 160: の別名である.
1.1 takayama 161: @end table
162:
1.6 takayama 163: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.1 takayama 164: @table @var
165: @item return
1.6 takayama 166: ベクトル, 超幾何関数の値とその微分. 詳しくは下記.
1.1 takayama 167: @item beta
1.6 takayama 168: 行和, 列和のリスト. 成分はすべて正であること.
1.1 takayama 169: @item p
1.6 takayama 170: 二元分割表のセルの確率のリスト
1.1 takayama 171: @end table
172:
1.6 takayama 173: @comment --- ここで関数の詳しい説明 ---
174: @comment --- @itemize〜@end itemize は箇条書き ---
175: @comment --- @bullet は黒点付き ---
1.1 takayama 176: @itemize @bullet
177: @item
1.6 takayama 178: gmvector は Gauss-Manin vector の略である [GM2016].
1.1 takayama 179: @item
1.6 takayama 180: gmvector の戻り値は
181: [GM2016] の 6章 p.23 のベクトル Sである.
182: これは
183: [GM2016] の４章で定義されているベクトル F の定数倍であり,
184: その定数は
185: 第一成分が [GM2016] の６章で定義されている級数 S の値と等しく
186: なるように決められている.
1.1 takayama 187: @item
1.6 takayama 188: r1 x r2 分割表を考える.
189: m+1=r1, n+1=r2 とおく.
190: 正規化定数 Z は分割表 u を (m+1) × (n+1) 行列とするとき p^u/u! の和である.
191: ここで和は行和列和が @var{beta} であるような u 全体でとる
1.1 takayama 192: [TKT2015], [GM2016].
1.6 takayama 193: S はこの多項式 Z の p を
1.1 takayama 194: @verbatim
195: [[1,y11,...,y1n],
196: [1,y21,...,y2n],...,
197: [1,ym1, ...,ymn],
198: [1,1, ..., 1]]
199: @end verbatim
1.6 takayama 200: (1 が L 字型に並ぶ),
201: と正規化した級数である.
1.1 takayama 202: @item
1.6 takayama 203: 2x(n+1)分割表で, gmvector の戻り値を Lauricella F_D で書くことが
204: 以下のようにできる
205: (b[2][1]-b[1][1] >= 0 の場合).
206: ここで b[1][1], b[1][2] は, それぞれ 1 行目の行和, 2 行目の行和,
207: b[2][i] は i 列目の列和である.
1.1 takayama 208: @comment ekn/Talks/2015-12-3-goto.tex
209: @verbatim
210: S=F_D(-b[1,1], [-b[2,2],...,-b[2,n+1]], b[2,1]-b[1,1]+1 ; y)/C,
211: @end verbatim
1.8 takayama 212: C=b[1,1]! b[2,2]! ... b[2,n+1]! (b[2,1]-b[1,1])!
1.6 takayama 213: とおく.
214: 1/C は L 字型の分割表
1.1 takayama 215: @verbatim
216: [[b[1,1], 0, ..., 0 ],
217: [b[2,1]-b[1,1],b[2,2], ..., b[2,n+1]]]
218: @end verbatim
1.6 takayama 219: に対応.
220: gmvector は
1.1 takayama 221: @verbatim
222: [S,(y11/a2) d_11 S,(y12/a3) d_12 S, ..., (y1n/a_(n+1)) d_1n S]
223: @end verbatim
1.6 takayama 224: である.
225: ここで d_ij は yij についての微分,
1.1 takayama 226: @verbatim
227: [a0, a1, ... ,a_(n+2)]
228: = [-b[1,2],-b[1,1],b[2,2], ..., b[2,n+1],b[2,1]]
229: @end verbatim
1.6 takayama 230: である.
1.1 takayama 231: @item
1.6 takayama 232: 周辺和 @var{beta}の時の正規化定数のセル確率 @var{p} に対する値は 多項式に退化した E(k,n) の値で表現できる. 文献 [TKT2015], [GM2016] 参照.
1.1 takayama 233: @item
1.17 takayama 234: 以下の option は expectation その他でも使える.
1.16 takayama 235: @item
1.6 takayama 236: option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう
1.1 takayama 237: [T2016].
1.6 takayama 238: 分散計算用の各種パラメータの設定は
1.19 takayama 239: gtt_ekn3.setup で行なう.
1.16 takayama 240: @item
241: option bs=1. binary splitting method で matrix factorial を計算.
242: 一般に 3x3 では効果あり(assert2(15|bs=1)), 5x5 (test5x5(20|bs=1))では遅くなる.
243: デフォールトは bs=0.
244: @item
245: option path. contiguity を適用する path をきめるアルゴリズムを指定.
1.17 takayama 246: path=2 (後藤, 松本の論文 [GM2016] の path). path=3 (論文 [TGKT] の path).
1.16 takayama 247: デフォールトは path=3.
1.17 takayama 248: @item
249: option interval. 通常の matrix factorial の計算では, 分母と分子をそれぞれ整数計算で計算し最後に約分をする. しかしながら数の中間膨張が一般的に発生しその中間膨張を解消するため
250: 約分を一定間隔で行うと計算効率がよくなる.
251: interval に整数値を設定することにより行列による一次変換を interval 回するたびに約分を行う.
252: interval の最適値は問題毎に異なるためシステムがデフォールト値を設定することはない.
1.18 takayama 253: @item
254: option x=1. subprocess 毎に window を開く.
1.1 takayama 255: @end itemize
256:
1.6 takayama 257: @comment --- @example〜@end example は実行例の表示 ---
258: 例: 次は2 x 2 分割表で行和が [5,1], 列和が [3,3], 各セルの確率が
259: [[1/2,1/3],[1/7,1/5]] の場合の gmvector の値である.
1.1 takayama 260: @example
1.19 takayama 261: [3000] load("gtt_ekn3.rr");
262: [3001] gtt_ekn3.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])
1.1 takayama 263: [775/27783]
264: [200/9261]
265: @end example
266:
1.8 takayama 267: 例: N を2以上の自然数とする時, Gauss の超幾何関数(この場合は多項式となる)
268: F(-36N,-11N,2N,(1-1/N)/56) の値は T3 に代入される ( [TGKT] ).
269: @comment ekn/Prog2/2x2.rr
270: @example
271: N=2;
1.19 takayama 272: T2=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],[[1,(1-1/N)/56],[1,1]])[0][0];
1.8 takayama 273: D=fac(36*N)*fac(11*N)*fac(2*N-1);
274: T3=T2*D;
275: @end example
276: ちなみに同じ値を Mathematica に計算させるには
277: @example
278: n=2; Hypergeometric2F1[-36*n,-11*n,2*n,(1-1/n)/56]
279: @end example
280:
1.17 takayama 281: 例: interval option
282: @example
283: [4009] P=gtt_ekn3.prob1(5,5,100);
284: [[[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]]]
285:
286: [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];
287: [cpu,72.852,gc,0,memory,4462742364,real,72.856]
288:
289: [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];
290: [cpu,67.484,gc,0,memory,3535280544,real,67.4844]
291: @end example
292:
1.6 takayama 293: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
294: 計算ができる.
295: 守備範囲の異なるプログラム同士の比較, debug 用参考.
1.1 takayama 296: @example
297: [3080] import("tk_fd.rr");
298: [3081] A=tk_fd.marginal2abc([4,5],[2,4,3]);
1.6 takayama 299: [-4,[-4,-3],-1] // 2変数 FD のパラメータ. a,[b1,b2],c
1.1 takayama 300: [3082] tk_fd.fd_hessian2(A[0],A[1],A[2],[1/2,1/3]);
301: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
302: [4483/124416,[ 1961/15552 185/1728 ],
303: [ 79/288 259/864 ]
304: [ 259/864 47/288 ]]
1.6 takayama 305: // 戻値は [F=F_D, gradient(F), Hessian(F)]
1.1 takayama 306:
1.6 takayama 307: // ekn_gt での例と同じパラメータ.
1.1 takayama 308: [3543] A=tk_fd.marginal2abc([5,1],[3,3]);
309: [-5,[-3],-1]
310: [3544] tk_fd.fd_hessian2(A[0],A[1],A[2],[(1/3)*(1/7)/((1/2)*(1/5))]);
311: Computing Dmat(ca) for parameters B=[-3],X=[ 10/21 ]
312: [775/27783,[ 20/147 ],[ 17/42 ]]
313: @end example
314:
1.6 takayama 315: 参考: 一般の A 分布の正規化定数についての Hessian の計算は実験的 package ot_hessian_ahg.rr
316: で実装のテストがされている. (これはまだ未完成のテスト版なので出力形式等も将来的には変更される.)
1.1 takayama 317: @example
318: import("ot_hgm_ahg.rr");
319: import("ot_hessian_ahg.rr");
320: def htest4() @{
321: extern C11_A;
322: extern C11_Beta;
323: Hess=newmat(7,7);
324: A =C11_A;
325: Beta0= [b0,b1,b2,b3];
326: BaseIdx=[4,5,6];
327: X=[x0,x1,x2,x3,x4,x5,x6];
328: for (I=0; I<7; I++) for (J=0; J<7; J++) @{
329: Idx = [I,J];
330: H=hessian_simplify(A,Beta0,X,BaseIdx,Idx);
331: Hess[I][J]=H;
332: printf("[I,J]=%a, Hessian_ij=%a\n",Idx,H);
333: @}
334: return(Hess);
335: @}
336: [2917] C11_A;
337: [[0,0,0,1,1,1,1],[1,0,0,1,0,1,0],[0,1,1,0,1,0,1],[1,1,0,1,1,0,0]]
338: [2918] C11_Beta;
339: [166,36,290,214]
340: [2919] Ans=htest4$
341: [2920] Ans[0][0];
342: [[((b1-b0-1)*x4)/(x0^2),[4]],[((b1-b0-1)*x6)/(x0^2),[6]],
343: [(b1^2+(-2*b0-1)*b1+b0^2+b0)/(x0^2),[]],[(x6)/(x0),[6,0]],[(x4)/(x0),[4,0]]]
344: @end example
345:
1.6 takayama 346: @comment --- 参照(リンク)を書く ---
1.1 takayama 347: @table @t
1.6 takayama 348: @item 参照
1.19 takayama 349: @ref{gtt_ekn3.setup}
350: @ref{gtt_ekn3.pfaffian_basis}
1.1 takayama 351: @end table
352:
1.6 takayama 353: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1 takayama 354: @noindent
355: ChangeLog
356: @itemize @bullet
357: @item
1.6 takayama 358: この関数は
359: [GM2016] のアルゴリズムおよび
1.17 takayama 360: [T2016] による modular method を用いた高速化,
361: [TGKT] の高速化を実装したものである.
1.1 takayama 362: @item
1.6 takayama 363: 変更を受けたファイルは
1.1 takayama 364: OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_pfaffian_8.rr
1.17 takayama 365: @item
366: interval option について変更を受けたファイルは
367: OpenXM/src/asir-contrib/packages/src/gtt_ekn3/ekn_eval.rr 1.6
1.1 takayama 368: @end itemize
369:
370:
371: @comment **********************************************************
1.19 takayama 372: @node gtt_ekn3.nc,,, 超幾何関数E(k,n)
373: @subsection @code{gtt_ekn3.nc}
1.6 takayama 374: @comment --- 索引用キーワード
1.19 takayama 375: @findex gtt_ekn3.nc
1.1 takayama 376:
377: @table @t
1.19 takayama 378: @item gtt_ekn3.nc(@var{beta},@var{p})
1.6 takayama 379: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z
380: およびその微分の値を戻す.
1.1 takayama 381: @end table
382:
1.6 takayama 383: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.1 takayama 384: @table @var
385: @item return
1.6 takayama 386: ベクトル [Z,[[d_11 Z, d_12 Z, ...], ..., [d_m1 Z, d_m2 Z, ...., d_mn Z]]]
1.1 takayama 387: @item beta
1.6 takayama 388: 行和, 列和のリスト. 成分はすべて正であること.
1.1 takayama 389: @item p
1.6 takayama 390: 二元分割表のセルの確率のリスト
1.1 takayama 391: @end table
392:
1.6 takayama 393: @comment --- ここで関数の詳しい説明 ---
394: @comment --- @itemize〜@end itemize は箇条書き ---
395: @comment --- @bullet は黒点付き ---
1.1 takayama 396: @itemize @bullet
397: @item
1.6 takayama 398: r1 x r2 分割表を考える.
399: m=r1, n=r2 とおく.
400: 正規化定数 Z は分割表 u を m × n 行列とするとき p^u/u! の和である.
401: ここで和は行和列和が @var{beta} であるような u 全体でとる
1.1 takayama 402: [TKT2015], [GM2016].
1.6 takayama 403: p^u は p_ij^u_ij の積, u! は u_ij! の積である.
404: d_ij Z で Z の変数 p_ij についての偏微分を表す.
1.1 takayama 405: @item
1.6 takayama 406: nc は gmvector の値を元に, [GM2016] の Prop
407: 7.1 に基づいて Z の値を計算する.
1.1 takayama 408: @item
1.6 takayama 409: option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
410: 分散計算用の各種パラメータの設定は
1.19 takayama 411: gtt_ekn3.setup で行なう.
1.17 takayama 412: その他の option は gmvector を参照.
1.1 takayama 413: @end itemize
414:
1.6 takayama 415: @comment --- @example〜@end example は実行例の表示 ---
416: 例: 2x3 分割表での Z とその微分の計算.
1.1 takayama 417: @example
1.19 takayama 418: [2237] gtt_ekn3.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
1.1 takayama 419: [4483/124416,[ 353/7776 1961/15552 185/1728 ]
420: [ 553/20736 1261/15552 1001/13824 ]]
421: @end example
422:
1.6 takayama 423: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
424: 計算ができる.
1.1 takayama 425: @example
426: [3076] import("tk_fd.rr");
427: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
428: [-4,[-4,-3],-1]
429: [3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
430: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
431: [ 1 1 1 ]
432: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
433: [4483/124416,[[353/7776,1961/15552,185/1728],
434: [553/20736,1261/15552,1001/13824]]]
1.6 takayama 435: // 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],
436: // [d_21 Z, d_22 Z, d_23 Z]]] の値.
437: // ここで d_ij は i,j 成分についての微分を表す.
1.1 takayama 438: @end example
439:
1.6 takayama 440: @comment --- 参照(リンク)を書く ---
1.1 takayama 441: @table @t
1.6 takayama 442: @item 参照
1.19 takayama 443: @ref{gtt_ekn3.setup}
444: @ref{gtt_ekn3.lognc}
1.1 takayama 445: @end table
446:
1.6 takayama 447: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1 takayama 448: @noindent
449: ChangeLog
450: @itemize @bullet
451: @item
1.6 takayama 452: 変更を受けたファイルは
1.1 takayama 453: OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1, gtt_ekn/ekn_eval.rr
454: @end itemize
455:
456:
457: @comment **********************************************************
1.19 takayama 458: @node gtt_ekn3.lognc,,, 超幾何関数E(k,n)
459: @subsection @code{gtt_ekn3.lognc}
1.6 takayama 460: @comment --- 索引用キーワード
1.19 takayama 461: @findex gtt_ekn3.lognc
1.1 takayama 462:
463: @table @t
1.19 takayama 464: @item gtt_ekn3.lognc(@var{beta},@var{p})
1.6 takayama 465: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の条件付き確率の正規化定数 Z
466: の log の近似値およびその微分の近似値を戻す.
1.1 takayama 467: @end table
468:
1.6 takayama 469: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.1 takayama 470: @table @var
471: @item return
1.6 takayama 472: ベクトル [log(Z), [[d_11 log(Z), d_12 log(Z), ...], [d_21 log(Z),...], ... ]
1.1 takayama 473: @item beta
1.6 takayama 474: 行和, 列和のリスト. 成分はすべて正であること.
1.1 takayama 475: @item p
1.6 takayama 476: 二元分割表のセルの確率のリスト
1.1 takayama 477: @end table
478:
1.6 takayama 479: @comment --- ここで関数の詳しい説明 ---
480: @comment --- @itemize〜@end itemize は箇条書き ---
481: @comment --- @bullet は黒点付き ---
1.1 takayama 482: @itemize @bullet
483: @item
1.6 takayama 484: 条件付き最尤推定に利用する [TKT2015].
485: @item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
486: 分散計算用の各種パラメータの設定は
1.19 takayama 487: gtt_ekn3.setup で行なう.
1.1 takayama 488: @end itemize
489:
1.6 takayama 490: @comment --- @example〜@end example は実行例の表示 ---
491: 例: 2 × 3 分割表での例. 第一成分のみ近似値.
1.1 takayama 492: @example
1.19 takayama 493: [2238] gtt_ekn3.lognc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
1.1 takayama 494: [-3.32333832422461674630,[ 5648/4483 15688/4483 13320/4483 ]
495: [ 3318/4483 10088/4483 9009/4483 ]]
496: @end example
497:
1.6 takayama 498: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
499: 計算ができる.
1.1 takayama 500: @example
501: [3076] import("tk_fd.rr");
502: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
503: [-4,[-4,-3],-1]
504: [3078] tk_fd.log_ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
505: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
506: [ 1 1 1 ]
507: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
508: [-3.32333832422461674639485797719209322217260539267246045320,
509: [[1.2598706, 3.499442, 2.971224],
510: [0.7401293, 2.250278, 2.009591]]]
1.6 takayama 511: // 戻値は [log(Z),
1.1 takayama 512: // [[d_11 log(Z), d_12 log(Z), d_13 log(Z)],
513: // [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]]
1.6 takayama 514: // の近似値.
1.1 takayama 515: @end example
516:
1.6 takayama 517: @comment --- 参照(リンク)を書く ---
1.1 takayama 518: @table @t
1.6 takayama 519: @item 参照
1.19 takayama 520: @ref{gtt_ekn3.setup}
521: @ref{gtt_ekn3.nc}
1.1 takayama 522: @end table
523:
1.6 takayama 524: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1 takayama 525: @noindent
526: ChangeLog
527: @itemize @bullet
528: @item
1.6 takayama 529: 変更を受けたファイルは
1.1 takayama 530: OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.
531: @end itemize
532:
533: @comment **********************************************************
1.19 takayama 534: @node gtt_ekn3.expectation,,, 超幾何関数E(k,n)
535: @subsection @code{gtt_ekn3.expectation}
1.6 takayama 536: @comment --- 索引用キーワード
1.19 takayama 537: @findex gtt_ekn3.expectation
1.1 takayama 538:
539: @table @t
1.19 takayama 540: @item gtt_ekn3.expectation(@var{beta},@var{p})
1.6 takayama 541: :: 周辺和 @var{beta}, セルの確率 @var{p} の二元分割表の期待値を計算する.
1.1 takayama 542: @end table
543:
1.6 takayama 544: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.1 takayama 545: @table @var
546: @item return
1.6 takayama 547: 二元分割表の各セルの期待値のリスト.
1.1 takayama 548: @item beta
1.6 takayama 549: 行和, 列和のリスト. 成分はすべて正であること.
1.1 takayama 550: @item p
1.6 takayama 551: 二元分割表のセルの確率のリスト
1.1 takayama 552: @end table
553:
1.6 takayama 554: @comment --- ここで関数の詳しい説明 ---
555: @comment --- @itemize〜@end itemize は箇条書き ---
556: @comment --- @bullet は黒点付き ---
1.1 takayama 557: @itemize @bullet
558: @item
1.17 takayama 559: [GM2016] の Algorithm 7.8 の実装. [TGKT] による高速化版 (path=3) がデフォールト.
1.6 takayama 560: @item option crt=1 (crt = Chinese remainder theorem) を与えると, 分散計算をおこなう.
561: 分散計算用の各種パラメータの設定は
1.19 takayama 562: gtt_ekn3.setup で行なう.
1.6 takayama 563: @item option index を与えると, 指定された成分の期待値のみ計算する.
564: たとえば 2 x 2 分割表で index=[[0,0],[1,1]] と指定すると, 1 のある成分の期待値のみ計算する.
1.17 takayama 565: @item その他の option は gmvector を参照.
1.1 takayama 566: @end itemize
567:
1.6 takayama 568: @comment --- @example〜@end example は実行例の表示 ---
1.1 takayama 569:
1.6 takayama 570: 2×2, 3×3 の分割表の期待値計算例.
1.1 takayama 571: @example
1.19 takayama 572: [2235] gtt_ekn3.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);
1.1 takayama 573: [ 2/3 1/3 ]
574: [ 4/3 8/3 ]
1.19 takayama 575: [2236] gtt_ekn3.expectation([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
1.1 takayama 576: [ 5648/4483 7844/4483 4440/4483 ]
577: [ 3318/4483 10088/4483 9009/4483 ]
578:
1.19 takayama 579: [2442] gtt_ekn3.expectation([[4,14,9],[11,6,10]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]]);
1.1 takayama 580: [ 207017568232262040/147000422096729819 163140751505489940/147000422096729819
581: 217843368649167296/147000422096729819 ]
582: [ 1185482401011137878/147000422096729819 358095302885438604/147000422096729819
583: 514428205457640984/147000422096729819 ]
584: [ 224504673820628091/147000422096729819 360766478189450370/147000422096729819
585: 737732646860489910/147000422096729819 ]
586: @end example
587:
1.6 takayama 588: 参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な
589: 計算ができる.
1.1 takayama 590: @example
591: [3076] import("tk_fd.rr");
592: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
593: [-4,[-4,-3],-1]
594: [3078] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
595: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
596: [ 1 1 1 ]
597: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
598: [[5648/4483,7844/4483,4440/4483],
599: [3318/4483,10088/4483,9009/4483]]
1.6 takayama 600: // 各セルの期待値.
1.1 takayama 601: @end example
602:
1.6 takayama 603: 参考: 一般の A 分布の計算は ot_hgm_ahg.rr. まだ実験的なため, module 化されていない.
604: ot_hgm_ahg.rr についての参考文献:
1.1 takayama 605: K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Systems II --- Holonomic Gradient Method, arxiv:1505.02947
606: @example
607: [3237] import("ot_hgm_ahg.rr");
1.6 takayama 608: // 2 x 2 分割表.
1.1 takayama 609: [3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],
610: [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
611: oohg_native=0, oohg_curl=1
612: [1376777025/625400597,1750225960/625400597,
613: 2375626557/625400597,3252978816/625400597]
1.6 takayama 614: // 2 x 2 分割表の期待値.
1.1 takayama 615:
1.6 takayama 616: // 2 x 3 分割表.
1.1 takayama 617: [3238] hgm_ahg_expected_values_contiguity(
618: [[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
619: [5,2,4,3],[1,1/2,1/3,1,1,1],[x1,x2,x3,x4,x5,x6]|geometric=1);
620: [5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
1.6 takayama 621: // 2 x 3 分割表の期待値. 上と同じ問題.
1.1 takayama 622: @end example
623:
1.6 takayama 624: 3 x 3 分割表. 構造的0が一つ.
1.1 takayama 625: @example
626: /*
1.6 takayama 627: dojo, p.221 のデータ. 成績3以下の生徒は集めてひとつに.
1.1 takayama 628: 2 1 1
629: 8 3 3
630: 0 2 6
631:
632: row sum: 4,14,8
633: column sum: 10,6,10
1.6 takayama 634: 0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
1.1 takayama 635: */
636:
637: A=[[0,0,0,1,1,1, 0,0],
638: [0,0,0,0,0,0, 1,1],
639: [1,0,0,1,0,0, 0,0],
640: [0,1,0,0,1,0, 1,0],
641: [0,0,1,0,0,1, 0,1]];
642: B=[14,8,10,6,10];
643: hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1],
1.6 takayama 644: [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);
1.1 takayama 645:
1.6 takayama 646: // 答.
1.1 takayama 647: [14449864949304/9556267369631,
648: 10262588586540/9556267369631, 13512615942680/9556267369631,
649: 81112808747006/9556267369631,
650: 21816297744346/9556267369631, 30858636683482/9556267369631,
651:
652: 25258717886900/9556267369631,51191421070148/9556267369631]
653: @end example
654:
1.6 takayama 655: 3 x 3 分割表.
1.1 takayama 656: @example
657: /*
1.6 takayama 658: 上のデータで 0 を 1 に変更.
1.1 takayama 659: 2 1 1
660: 8 3 3
661: 1 2 6
662:
663: row sum: 4,14,9
664: column sum: 11,6,10
665: */
666: A=[[0,0,0,1,1,1,0,0,0],
667: [0,0,0,0,0,0,1,1,1],
668: [1,0,0,1,0,0,1,0,0],
669: [0,1,0,0,1,0,0,1,0],
670: [0,0,1,0,0,1,0,0,1]];
671: B=[14,9,11,6,10];
672: hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1,1],
673: [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);
674:
1.6 takayama 675: // 期待値, 答. x9 を指定していないので, 9番目の期待値は出力してない.
1.1 takayama 676: [207017568232262040/147000422096729819,
677: 163140751505489940/147000422096729819,217843368649167296/147000422096729819,
678: 1185482401011137878/147000422096729819,
679: 358095302885438604/147000422096729819,514428205457640984/147000422096729819,
680: 224504673820628091/147000422096729819,360766478189450370/147000422096729819]
681:
1.6 takayama 682: // Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
683: // まだ書いてない.
1.1 takayama 684: @end example
685:
686:
687:
1.6 takayama 688: @comment --- 参照(リンク)を書く ---
1.1 takayama 689: @table @t
1.6 takayama 690: @item 参照
1.19 takayama 691: @ref{gtt_ekn3.setup}
692: @ref{gtt_ekn3.nc}
1.1 takayama 693: @end table
694:
1.6 takayama 695: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1 takayama 696: @noindent
697: ChangeLog
698: @itemize @bullet
699: @item
1.6 takayama 700: 変更を受けたファイルは
1.1 takayama 701: OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1.
702: @end itemize
703:
704:
705: @comment **********************************************************
1.6 takayama 706: @comment --- ◯◯◯◯ の説明
707: @comment --- 個々の関数の説明の開始 ---
708: @comment --- section 名を正確に ---
1.19 takayama 709: @node gtt_ekn3.setup,,, 超幾何関数E(k,n)
710: @subsection @code{gtt_ekn3.setup}
1.6 takayama 711: @comment --- 索引用キーワード
1.19 takayama 712: @findex gtt_ekn3.setup
1.1 takayama 713:
714: @table @t
1.19 takayama 715: @item gtt_ekn3.setup()
1.6 takayama 716: :: 分散計算用の環境設定をおこなう. 現在の環境を報告する.
1.1 takayama 717: @end table
718:
1.6 takayama 719: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.1 takayama 720: @table @var
721: @item return
722:
723: @end table
724:
1.6 takayama 725: @comment --- ここで関数の詳しい説明 ---
726: @comment --- @itemize〜@end itemize は箇条書き ---
727: @comment --- @bullet は黒点付き ---
1.3 takayama 728: @itemize @bullet
1.6 takayama 729: @item 使用するプロセスと素数の個数, 最小の素数を表示する. 準備されていない場合はその旨を表示.
730: @item このパッケージでの分散計算は複数のcpuを搭載した計算機で実行されることを想定している.
731: @item option nps (または number_of_processes)を与えると指定した数だけプロセスを用意する.
732: @item option nprm (または number_of_primes)を与えるとnprmが文字列の場合指定された素数リストのファイルを読み込む. nprmが自然数の場合さらにoption minp (minp =MINimum Prime)を与えるとminpより大きな素数をnprm個生成する. その際option fgp (または file_of_generated_primes)を与えると生成した素数リストをファイル名をfgpとして保存する.
733: @item 上記のoption を指定しなかった場合次のデフォルト値が用いられる. nps=1. nprm=10. fgp=0.
734: @item option report=1を与えると現在の環境の報告のみを行う. setup(|report=1)の別名としてreport関数を使用することもできる.
1.19 takayama 735: @item option subprogs=[file1,file2,...] により分散計算の子供プロセスにロードすべきファイル file1, file2, ... を指定する. default は subprogs=["gtt_ekn3/childprocess.rr"] である.
736: @item gtt_ekn3.set_debug_level(Mode) で Ekn_debug の値を設定する.
1.1 takayama 737: @end itemize
738:
1.6 takayama 739: @comment --- @example〜@end example は実行例の表示 ---
740: 例: 素数のリストを生成してファイル p.txt へ書き出す.
1.1 takayama 741: @example
1.19 takayama 742: gtt_ekn3.setup(|nps=2,nprm=20,minp=10^10,fgp="p.txt")$
1.1 takayama 743: @end example
744:
1.8 takayama 745: 例: chinese remainder theorem (crt) を使って gmvector を計算.
746: @example
1.19 takayama 747: [2867] gtt_ekn3.setup(|nprm=20,minp=10^20);
748: [2868] N=2; T2=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
1.8 takayama 749: [[1,(1-1/N)/56],[1,1]] | crt=1)$
750: @end example
751:
1.1 takayama 752:
1.6 takayama 753: @comment --- 参照(リンク)を書く ---
1.1 takayama 754: @table @t
1.6 takayama 755: @item 参照
1.19 takayama 756: @ref{gtt_ekn3.nc}
757: @ref{gtt_ekn3.gmvector}
1.1 takayama 758: @end table
759:
1.6 takayama 760: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1 takayama 761: @noindent
762: ChangeLog
763: @itemize @bullet
764: @item
1.6 takayama 765: 変更を受けたファイルは
1.1 takayama 766: OpenXM/src/asir-contrib/packages/src/gtt_ekn.rr 1.1,
767: gtt_ekn/g_mat_fac.rr
768:
769: @end itemize
770:
771: @comment **********************************************************
1.6 takayama 772: @comment --- ◯◯◯◯ の説明
773: @comment --- 個々の関数の説明の開始 ---
774: @comment --- section 名を正確に ---
1.19 takayama 775: @node gtt_ekn3.upAlpha,,, 超幾何関数E(k,n)
776: @node gtt_ekn3.downAlpha,,, 超幾何関数E(k,n)
777: @subsection @code{gtt_ekn3.upAlpha}, @code{gtt_ekn3.downAlpha}
1.6 takayama 778: @comment --- 索引用キーワード
1.19 takayama 779: @findex gtt_ekn3.upAlpha
780: @findex gtt_ekn3.downAlpha
1.1 takayama 781:
782: @table @t
1.19 takayama 783: @item gtt_ekn3.upAlpha(@var{i},@var{k},@var{n})
784: @item gtt_ekn3.downAlpha(@var{i},@var{k},@var{n})
1.1 takayama 785: ::
786: @end table
787:
1.6 takayama 788: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.1 takayama 789: @table @var
1.12 takayama 790: @item i a_i を a_i+1 (a_i を a_i-1) と変化させる contiguity relation.
1.6 takayama 791: @item k E(k+1,n+k+2)型の超幾何関数の k. 分割表では (k+1)×(n+1).
792: @item n E(k+1,n+k+2)型の超幾何関数の n. 分割表では (k+1)×(n+1).
793: @item return contiguity relation の pfaffian_basis についての行列表現を戻す. [GM2016] の Cor 6.3.
1.1 takayama 794: @end table
795:
1.6 takayama 796: @comment --- ここで関数の詳しい説明 ---
797: @comment --- @itemize〜@end itemize は箇条書き ---
798: @comment --- @bullet は黒点付き ---
1.1 takayama 799: @itemize @bullet
800: @item
1.6 takayama 801: upAlpha は [GM2016] の Cor 6.3 の行列 U_i を戻す.
802: @item 関連する各関数の簡潔な説明と例も加える.
803: @item a_i を a_i-1 と変化させたい場合は関数 downAlpha を用いる.
804: @item a_i と分割表の周辺和を見るには, 関数 marginaltoAlpha([行和,列和]) を用いる.
1.1 takayama 805: @item
1.6 takayama 806: pfaffian_basis は [GM2016] の４章のベクトル F に対応する偏微分を戻す.
1.12 takayama 807: @item optional 引数 arule, xrule で a_i や x_i_j を数にしたものをより効率的に求めることができる. 変化をうけるパラメータを数にしてしまっても特にエラー表示はしない. a_0 で和の条件を調整しているので注意(Todo, double check).
1.1 takayama 808: @end itemize
809:
1.6 takayama 810: @comment --- @example〜@end example は実行例の表示 ---
811: 例: 以下の例は 2×2分割表(E(2,4)), 2×3分割表(E(2,5))の場合である.
812: [2225] までは出力を略している.
1.1 takayama 813: @example
1.19 takayama 814: [2221] gtt_ekn3.marginaltoAlpha([[1,4],[2,3]]);
1.1 takayama 815: [[a_0,-4],[a_1,-1],[a_2,3],[a_3,2]]
1.19 takayama 816: [2222] gtt_ekn3.upAlpha(1,1,1); // E(2,4) の a_1 方向の
1.6 takayama 817: // contiguity を表現する行列
1.19 takayama 818: [2223] gtt_ekn3.upAlpha(2,1,1); // E(2,4) の a_2 方向
819: [2224] gtt_ekn3.upAlpha(3,1,1); // E(2,4) の a_3 方向
1.1 takayama 820: [2225] function f(x_1_1);
1.19 takayama 821: [2232] gtt_ekn3.pfaffian_basis(f(x_1_1),1,1);
1.1 takayama 822: [ f(x_1_1) ]
823: [ (f{1}(x_1_1)*x_1_1)/(a_2) ]
824: [2233] function f(x_1_1,x_1_2);
825: f() redefined.
1.19 takayama 826: [2234] gtt_ekn3.pfaffian_basis(f(x_1_1,x_1_2),1,2); // E(2,5), 2*3 分割表
1.1 takayama 827: [ f(x_1_1,x_1_2) ]
828: [ (f{1,0}(x_1_1,x_1_2)*x_1_1)/(a_2) ]
829: [ (f{0,1}(x_1_1,x_1_2)*x_1_2)/(a_3) ]
1.12 takayama 830:
831: [2235] RuleA=[[a_2,1/3],[a_3,1/2]]$ RuleX=[[x_1_1,1/5]]$
1.19 takayama 832: base_replace(gtt_ekn3.upAlpha(1,1,1),append(RuleA,RuleX))
833: -gtt_ekn3.upAlpha(1,1,1 | arule=RuleA, xrule=RuleX);
1.12 takayama 834:
835: [ 0 0 ]
836: [ 0 0 ]
837:
1.1 takayama 838: @end example
839:
840:
1.6 takayama 841: @comment --- 参照(リンク)を書く ---
1.1 takayama 842: @table @t
1.6 takayama 843: @item 参照
1.19 takayama 844: @ref{gtt_ekn3.nc}
845: @ref{gtt_ekn3.gmvector}
1.1 takayama 846: @end table
847:
1.6 takayama 848: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.1 takayama 849: @noindent
850: ChangeLog
851: @itemize @bullet
852: @item
1.6 takayama 853: この関数は [GM2016]
854: で与えられたアルゴリズムに従い contiguity relation を導出する.
1.1 takayama 855: @item
1.6 takayama 856: 変更を受けたファイルは
1.1 takayama 857: OpenXM/src/asir-contrib/packages/src/gtt_ekn/ekn_pfaffian_8.rr 1.1.
858: @end itemize
859:
860:
1.5 takayama 861: @comment **********************************************************
1.6 takayama 862: @comment --- ◯◯◯◯ の説明
863: @comment --- 個々の関数の説明の開始 ---
864: @comment --- section 名を正確に ---
1.19 takayama 865: @node gtt_ekn3.cmle,,, 超幾何関数E(k,n)
866: @subsection @code{gtt_ekn3.cmle}
1.6 takayama 867: @comment --- 索引用キーワード
1.19 takayama 868: @findex gtt_ekn3.cmle
1.5 takayama 869:
870: @table @t
1.19 takayama 871: @item gtt_ekn3.cmle(@var{u}) u を観測データとするとき, P(U=u | row sum, column sum = these of U) を最大化する, 各セルの確率の近似値を求める.
1.5 takayama 872: ::
873: @end table
874:
1.6 takayama 875: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.5 takayama 876: @table @var
1.6 takayama 877: @item u 観測データ(分割表)
878: @item return セルの確率(分割表形式)
1.5 takayama 879: @end table
880:
1.6 takayama 881: @comment --- ここで関数の詳しい説明 ---
882: @comment --- @itemize〜@end itemize は箇条書き ---
883: @comment --- @bullet は黒点付き ---
1.5 takayama 884: @itemize @bullet
1.6 takayama 885: @item u を観測データとするとき, P(U=u | row sum, column sum = these of U) を最大化する, 各セルの確率の近似値を求める.
886: @item optional parameter で algorithm の振る舞い(たとえば有理数を近似して, 分母分子が小さい有理数にする, gradient descent の step幅)を調整すべきだが, これは作業中. 2017.03.03
1.5 takayama 887: @end itemize
888:
1.6 takayama 889: @comment --- @example〜@end example は実行例の表示 ---
890: 例: 2 x 4 分割表.
1.5 takayama 891: @example
892: U=[[1,1,2,3],[1,3,1,1]];
1.19 takayama 893: gtt_ekn3.cmle(U);
1.5 takayama 894: [[ 1 1 2 3 ]
895: [ 1 3 1 1 ],[[7,6],[2,4,3,4]], // Data, row sum, column sum
896: [ 1 67147/183792 120403/64148 48801/17869 ] // probability obtained.
897: [ 1 1 1 1 ]]
898: @end example
899:
1.6 takayama 900: 例: 上の例は次の関数に.
1.5 takayama 901: @example
1.19 takayama 902: gtt_ekn3.cmle_test3();
1.5 takayama 903: @end example
904:
1.6 takayama 905: @comment --- 参照(リンク)を書く ---
1.5 takayama 906: @table @t
1.6 takayama 907: @item 参照
1.19 takayama 908: @ref{gtt_ekn3.expectation}
1.5 takayama 909: @end table
910:
1.6 takayama 911: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.5 takayama 912: @noindent
913: ChangeLog
914: @itemize @bullet
1.19 takayama 915: @item gtt_ekn3/mle.rr に本体がある.
916: @item gtt_ekn3.rr の cmle 関数は wrapper.
1.5 takayama 917: @end itemize
918: @comment end cmle.
919:
1.8 takayama 920: @comment **********************************************************
921: @comment --- ◯◯◯◯ の説明
922: @comment --- 個々の関数の説明の開始 ---
923: @comment --- section 名を正確に ---
1.19 takayama 924: @node gtt_ekn3.set_debug_level,,, 超幾何関数E(k,n)
925: @node gtt_ekn3.contiguity_mat_list_2,,, 超幾何関数E(k,n)
926: @node gtt_ekn3.show_path,,, 超幾何関数E(k,n)
927: @node gtt_ekn3.get_svalue,,, 超幾何関数E(k,n)
928: @node gtt_ekn3.assert1,,, 超幾何関数E(k,n)
929: @node gtt_ekn3.assert2,,, 超幾何関数E(k,n)
930: @node gtt_ekn3.assert3,,, 超幾何関数E(k,n)
931: @node gtt_ekn3.prob1,,, 超幾何関数E(k,n)
932: @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}
1.8 takayama 933: @comment --- 索引用キーワード
1.19 takayama 934: @findex gtt_ekn3.set_debug_level
935: @findex gtt_ekn3.contiguity_mat_list_2
936: @findex gtt_ekn3.show_path
937: @findex gtt_ekn3.get_svalue
938: @findex gtt_ekn3.assert1
939: @findex gtt_ekn3.assert2
940: @findex gtt_ekn3.assert3
941: @findex gtt_ekn3.prob1
942:
943: @table @t
944: @item gtt_ekn3.set_debug_level(@var{m}) debug メッセージのレベルを設定.
945: @item gtt_ekn3.contiguity_mat_list_2 使用する contiguity を構成.
946: @item gtt_ekn3.show_path() どのように contiguity を適用したかの情報.
947: @item gtt_ekn3.get_svalue() static 変数の値を得る.
948: @item gtt_ekn3.assert1(@var{N}) 2x2 分割表で動作チェック.
949: @item gtt_ekn3.assert2(@var{N}) 3x3 分割表で動作チェック.
950: @item gtt_ekn3.assert3(@var{R1}, @var{R2}, @var{Size}) R1 x R2 分割表で並列動作の動作チェック.
951: @item gtt_ekn3.prob1(@var{R1},@var{R2},@var{Size}) R1 x R2 分割表用のテストデータを作る.
1.8 takayama 952: ::
953: @end table
954:
955: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
956: @table @var
957: @item m レベル.
958: @end table
959:
960: @comment --- ここで関数の詳しい説明 ---
961: @comment --- @itemize〜@end itemize は箇条書き ---
962: @comment --- @bullet は黒点付き ---
963: @itemize @bullet
1.19 takayama 964: @item (@var{m} & 0x1) == 0x1 の時 g_mat_fac_test_plain と g_mat_fac_itor の両方を呼び出し値を比較する (gtt_ekn3.setup した状態で).
1.11 takayama 965: @item (@var{m} & 0x2) == 0x2 の時 g_mat_fac_test への引数を tmp-input-数.ab として保存.
1.8 takayama 966: @item (@var{m} & 0x4) == 0x4 の時 matrix factorial の計算の呼び出し引数を表示.
1.10 takayama 967: @item @var{N} は問題の周辺和のサイズ.
1.12 takayama 968: @item @code{get_svalue} の戻り値は @code{[Ekn_plist,Ekn_IDL,Ekn_debug,Ekn_mesg,XRule,ARule,Verbose,Ekn_Rq]} の値.
1.18 takayama 969: @item assert3 の options: x=1, subprocess の window を表示. nps=m, m 個のプロセスで contiguity を求める (contiguity_mat_list_3). crt, interval などは gmvector などと共通の
970: option. timing data を表示するには load("gtt_ekn3/ekn_eval-timing.rr"); しておく.
1.8 takayama 971: @end itemize
972:
973: @comment --- @example〜@end example は実行例の表示 ---
1.10 takayama 974: 例.
1.8 takayama 975: @example
1.19 takayama 976: [2846] gtt_ekn3.set_debug_level(0x4);
977: [2847] N=2; T2=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
1.8 takayama 978: [[1,(1-1/N)/56],[1,1]])$
979: [2848] level&0x4: g_mat_fac_test([ 113/112 ]
980: [ 1/112 ],[ (t+225/112)/(t^2+4*t+4) (111/112*t+111/112)/(t^2+4*t+4) ]
981: [ (1/112)/(t^2+4*t+4) (111/112*t+111/112)/(t^2+4*t+4) ],0,20,1,t)
1.19 takayama 982: Note: we do not use g_mat_fac_itor. Call gtt_ekn3.setup(); to use the crt option.
1.8 takayama 983: level&0x4: g_mat_fac_test([ 67/62944040755546030080000 ]
984: [ 1/125888081511092060160000 ],[ (t+24)/(t^2+25*t+46) (2442)/(t^2+25*t+46) ]
985: [ (1)/(t^2+25*t+46) (-111*t-111)/(t^2+25*t+46) ],0,73,1,t)
986: level&0x4: g_mat_fac_test ------ snip
987: @end example
988:
1.10 takayama 989: 例.
1.9 takayama 990: @example
1.19 takayama 991: [2659] gtt_ekn3.nc([[4,5,6],[2,4,9]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]])$
992: [2660] L=matrix_transpose(gtt_ekn3.show_path())$
1.9 takayama 993: [2661] L[2];
1.19 takayama 994: [2 1]
1.9 takayama 995: @end example
1.19 takayama 996: [2 1] の index をもつパラメーター alpha の方向の contigity を求めそれを掛けて
1.9 takayama 997: 計算したことがわかる. L[0] は用いた contiguity の行列.
1.10 takayama 998: L[1] はcontiguity を適用する step 数.
999:
1000: 例. 値を計算せずに path のみ求めたい場合.
1001: @example
1.19 takayama 1002: A=gtt_ekn3.marginaltoAlpha_list([[400,410,1011],[910,411,500]])$
1003: [2666] gtt_ekn3.contiguity_mat_list_2(A,2,2)$
1004: [2667] L=matrix_transpose(gtt_ekn3.show_path())$
1.10 takayama 1005: [2668] L[2];
1006: [ 2 1 5 4 3 ]
1.19 takayama 1007: [2669] gtt_ekn3.contiguity_mat_list_3(A,2,2)$ // new alg in [TGKT]
1008: [2670] L=matrix_transpose(gtt_ekn3.show_path())$
1009: [2671] L[2];
1010: [2 1] // shorter
1.10 takayama 1011: @end example
1012:
1.15 takayama 1013: 例. 値を計算せずに path のみ求めたい場合.
1014: gtt_ekn3 による新しいアルゴリズムによる path の表示.
1015: @example
1016: A=gtt_ekn3.marginaltoAlpha_list([[10,20],[15,15]])$
1017: [2666] gtt_ekn3.contiguity_mat_list_3(A,1,1 | xrule=[[x_1_1,1/2]])$
1018: [t,[[ (-t-43/2)/(t-2) (-15/2)/(t-2) ]
1019: [ 1/2 -1/2 ],-9]]
1020: @end example
1021:
1.10 takayama 1022: 例. 0 が戻れば g_mat_fac_plain と指定した計算方法の結果が一致したことがわかる.
1023: option を書かないと g_mat_fac_int との比較となる.
1024: @example
1.19 takayama 1025: [8859] gtt_ekn3.assert2(1);
1.10 takayama 1026: Marginal=[[130,170,353],[90,119,444]]
1027: P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]
1.19 takayama 1028: Try g_mat_fac_test_int: Note: we do not use g_mat_fac_itor. Call gtt_ekn3.setup(); to use the crt option.
1.10 takayama 1029: Timing (int) =0.413916 (CPU) + 0.590723 (GC) = 1.00464 (total), real time=0.990672
1030:
1.19 takayama 1031: Try g_mat_fac_test_plain: Note: we do not use g_mat_fac_itor. Call gtt_ekn3.setup(); to use the crt option.
1.10 takayama 1032: Timing (rational) =4.51349 (CPU) + 6.32174 (GC) = 10.8352 (total)
1033: diff of both method =
1034: [ 0 0 0 ]
1035: [ 0 0 0 ]
1036: [ 0 0 0 ]
1037: [8860]
1038:
1.19 takayama 1039: [8863] gtt_ekn3.setup(|nprm=100,minp=10^50);
1.10 takayama 1040: Number of processes = 1.
1041: Number of primes = 100.
1042: Min of plist = 100000000000000000000000000000000000000000000000151.
1043: 0
1.19 takayama 1044: [8864] gtt_ekn3.assert2(1 | crt=1);
1.10 takayama 1045: Marginal=[[130,170,353],[90,119,444]]
1046: P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]
1047: Try [[crt,1]]
1048: ---- snip
1049: @end example
1050: なお二番目の例の timing (total) [例では省略] は mod 計算を subprocess がやっているので正しい値ではない. real time が計算時間の目安になる.
1.9 takayama 1051:
1.11 takayama 1052: 例.
1053: @example
1.17 takayama 1054: 3x5 分割表. 周辺和は 10 に比例する一定の数(factor option も関係. ソースを参照).
1055: cell 確率は1/素数で生成される.
1056: @comment grep testnxn ekn/Prog2/*.rr ; grep test_nxn ekn/Prog2/*.rr も見よ.
1.19 takayama 1057: [9054] L=gtt_ekn3.prob1(3,5,10 | factor=1, factor_row=3);
1.11 takayama 1058: [[[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]]]
1.19 takayama 1059: [9055] number_eval(gtt_ekn3.expectation(L[0],L[1]));
1060: [ 1.65224223218613 ... snip ]
1.11 takayama 1061: @end example
1062:
1.18 takayama 1063: 例:
1064: @example
1065: [5779] import("gtt_ekn3.rr"); load("gtt_ekn3/ekn_eval-timing.rr");
1066: [5780] gtt_ekn3.assert3(5,5,100 | nps=32, interval=100);
1067: -- snip
1068: Parallel method: Number of process=32, File name tmp-gtt_ekn3/p300.txt is written.
1069: Number of processes = 32.
1070: -- snip
1071: initialPoly of path=3: [ 2.184 0 124341044 2.1831 ] [CPU(s),0,*,real(s)]
1072: contiguity_mat_list_3 of path=3: [ 0.04 0 630644 9.6774 ] [CPU(s),0,*,real(s)]
1073: 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.
1074: g_mat_fac of path=3: [ 21.644 0 1863290168 21.6457 ] [CPU(s),0,*,real(s)]
1075: Done. Saved in 2.ab
1076: Diff (should be 0)=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,..., 0,0,0]
1077: @end example
1078:
1.8 takayama 1079: @comment --- 参照(リンク)を書く ---
1080: @table @t
1081: @item 参照
1.19 takayama 1082: @ref{gtt_ekn3.nc}
1.8 takayama 1083: @end table
1084:
1085: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1086: @noindent
1087: ChangeLog
1088: @itemize @bullet
1.19 takayama 1089: @item gtt_ekn3/ekn_eval.rr で matrix factorial の計算の呼び出し引数を表示する.
1.8 takayama 1090: @item grep 'iand(Ekn_debug,0x1)' *.rr でソースコードの該当の位置をさがす.
1091: @end itemize
1092: @comment end set_debug_level
1093:
1.5 takayama 1094:
1095:
1.19 takayama 1096: @node modular計算,,, Top
1.6 takayama 1097: @chapter modular計算
1.4 takayama 1098:
1099: @menu
1.19 takayama 1100: * gtt_ekn3.chinese_itor::
1.4 takayama 1101: @end menu
1102:
1.6 takayama 1103: @node 中国剰余定理とitor,,, modular計算
1104: @section 中国剰余定理とitor
1.4 takayama 1105:
1106: @comment **********************************************************
1.6 takayama 1107: @comment --- ◯◯◯◯ の説明
1108: @comment --- 個々の関数の説明の開始 ---
1109: @comment --- section 名を正確に ---
1.19 takayama 1110: @node gtt_ekn3.chinese_itor,,,
1111: @subsection @code{gtt_ekn3.chinese_itor}
1.6 takayama 1112: @comment --- 索引用キーワード
1.19 takayama 1113: @findex gtt_ekn3.chinese_itor 中国剰余定理とitor
1.4 takayama 1114:
1115: @table @t
1.19 takayama 1116: @item gtt_ekn3.chinese_itor(@var{data},@var{idlist})
1.6 takayama 1117: :: mod p で計算した結果(ベクトル)から chinese remainder theorem, itor(integer to rational) で有理数ベクトルを得る.
1.4 takayama 1118: @end table
1119:
1.6 takayama 1120: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.4 takayama 1121: @table @var
1.6 takayama 1122: @item return [val, n] ここで val は答え. また, n = n1*n2*...
1123: @item data [[val1,n1],[val2,n2], ...], ここで val mod n1 = val1, val mod n2 = val2,...
1124: @item idlist chinese, itor を実行するサーバIDのリスト.
1.4 takayama 1125: @end table
1126:
1.6 takayama 1127: @comment --- ここで関数の詳しい説明 ---
1128: @comment --- @itemize〜@end itemize は箇条書き ---
1129: @comment --- @bullet は黒点付き ---
1.4 takayama 1130: @itemize @bullet
1.6 takayama 1131: @item 中国剰余定理を用いて val0 mod n1 = val1, val0 mod n2 = val2, ... となる val0 を求める. val に algorithm itor を適用する.
1132: @item sqrt(n) より val0 が大きい時は itor が適用されて val0 が有理数 val=a/b に変換される. つまり b*x =1 mod n となる逆数 x を考えて, x*a % n = val0 となる数 val を戻す. 見つからないときは failure を戻す.
1.4 takayama 1133: @end itemize
1134:
1.6 takayama 1135: @comment --- @example〜@end example は実行例の表示 ---
1136: 例: [3!, 5^3*3!]=[6,750] が戻り値.
1137: 6 mod 109 =6, 750 mod 109=96 が最初の引数の [[6,96],109]. 以下同様.
1.4 takayama 1138: @example
1.19 takayama 1139: gtt_ekn3.setup(|nps=2,nprm=3,minp=101,fgp="p_small.txt");
1140: SS=gtt_ekn3.get_svalue();
1.4 takayama 1141: SS[0];
1142: [103,107,109] // list of primes
1143: SS[1];
1144: [0,2] // list of server ID's
1.19 takayama 1145: gtt_ekn3.chinese_itor([[[ 6,96 ],109],[[ 6,29 ],103],[[ 6,1 ],107]],SS[1]);
1.4 takayama 1146: [[ 6 750 ],1201289]
1147:
1.6 takayama 1148: // 引数はスカラーでもよい.
1.19 takayama 1149: gtt_ekn3.chinese_itor([[96,109],[29,103]],SS[1]);
1.4 takayama 1150: [[ 750 ],11227]
1151: @end example
1152:
1153:
1.6 takayama 1154: @comment --- @example〜@end example は実行例の表示 ---
1.19 takayama 1155: 例: gtt_ekn3/childprocess.rr (server で実行される) の関数 chinese (chinese remainder theorem) と euclid.
1.4 takayama 1156: @example
1.19 takayama 1157: load("gtt_ekn3/childprocess.rr");
1.4 takayama 1158: chinese([newvect(2,[6,29]),103],[newvect(2,[6,750]),107*109]);
1.6 takayama 1159: // mod 103 で [6,29], mod (107*109) で [6,750] となる数を mod 103*(107*109)
1160: // で求めると,
1.4 takayama 1161: [[ 6 750 ],1201289]
1.6 takayama 1162: euclid(3,103); // mod 103 での 3 の逆数. つまり 1/3
1.4 takayama 1163: -34
1.6 takayama 1164: 3*(-34) % 103; // 確かに逆数.
1.4 takayama 1165: 1
1166: @end example
1167:
1.6 takayama 1168: @comment --- @example〜@end example は実行例の表示 ---
1.19 takayama 1169: 例: gtt_ekn3/childprocess.rr (server で実行される) の関数 itor (integer to rational) の例.
1.6 takayama 1170: itor(Y,Q,Q2,Idx) では Y < Q2 なら Y がそのまま戻る. Idx は 内部用の index で好きな数でよい. 戻り値の第２成分となる.
1.4 takayama 1171: @example
1.19 takayama 1172: load("gtt_ekn3/childprocess.rr");
1.4 takayama 1173: for (I=1;I<11; I++) print([I,itor(I,11,3,0)]);
1174: [1,[1,0]]
1175: [2,[2,0]]
1.6 takayama 1176: [3,[-2/3,0]] //euclid(3,11); ->4, 4*(-2)%11 -> 3 なので確かに -2/3 は元の数の候補
1.4 takayama 1177: [4,[failure,0]]
1178: [5,[-1/2,0]]
1179: [6,[1/2,0]]
1180: [7,[-1/3,0]]
1181: [8,[failure,0]]
1182: [9,[-2,0]]
1183: [10,[-1,0]]
1184: @end example
1185:
1186:
1.6 takayama 1187: @comment --- 参照(リンク)を書く ---
1.4 takayama 1188: @table @t
1.6 takayama 1189: @item 参照
1.19 takayama 1190: @ref{gtt_ekn3.setup}
1.4 takayama 1191: @end table
1192:
1.6 takayama 1193: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.4 takayama 1194: @noindent
1195: ChangeLog
1196: @itemize @bullet
1197: @item
1.6 takayama 1198: 関連ファイルは
1.19 takayama 1199: gtt_ekn3/g_mat_fac.rr
1200: gtt_ekn3/childprocess.rr
1.4 takayama 1201: @end itemize
1202:
1.19 takayama 1203: @node Binary splitting,,, Top
1204: @chapter Binary splitting
1.14 takayama 1205:
1206: @menu
1207: * gtt_ekn3.init_dm_bsplit::
1208: * gtt_ekn3.setup_dm_bsplit::
1209: * gtt_ekn3.init_bsplit::
1210: @end menu
1211:
1.19 takayama 1212: @node matrix factorial,,, Binary splitting
1.14 takayama 1213: @section matrix factorial
1214:
1215: @comment **********************************************************
1216: @comment --- ◯◯◯◯ の説明
1217: @comment --- 個々の関数の説明の開始 ---
1218: @comment --- section 名を正確に ---
1219: @node gtt_ekn3.init_bsplit,,,
1220: @node gtt_ekn3.init_dm_bsplit,,,
1221: @node gtt_ekn3.setup_dm_bsplit,,,
1222: @subsection @code{gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit}
1223: @comment --- 索引用キーワード
1224: @findex gtt_ekn3.init_dm_bsplit matrix factorial
1225: @findex gtt_ekn3.setup_dm_bsplit matrix factorial
1226: @findex gtt_ekn3.init_bsplit matrix factorial
1227:
1228: @table @t
1229: @item gtt_ekn3.init_bsplit(|minsize=16,levelmax=1);
1230: :: binary split の実行のためのパラメータを設定する.
1231: @item gtt_ekn3.init_dm_bsplit(|bsplit_x=0, bsplit_reduce=0)
1232: :: binary split の分散実行のためのパラメータを設定する.
1233: @item gtt_ekn3.setup_dm_bsplit(C)
1234: :: binary split の分散実行のために C 個のプロセスを立ち上げる.
1235: @end table
1236:
1237: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1238: @table @var
1.15 takayama 1239: @item C はlevelmax-1 に設定する. 特に levalmax=1 のときは分散計算を行わない.
1240: @item bsplit_x=1 のとき, debug 用に各プロセスを xterm で表示.
1.14 takayama 1241: @end table
1242:
1243: @comment --- ここで関数の詳しい説明 ---
1244: @comment --- @itemize〜@end itemize は箇条書き ---
1245: @comment --- @bullet は黒点付き ---
1246: @itemize @bullet
1247: @item expectation などの関数に bs=1 オプションを与えると matrix factorial を binary
1248: splitting method で計算する.
1249: @end itemize
1250:
1251: @comment --- @example〜@end example は実行例の表示 ---
1.15 takayama 1252: 例: bs=1 と無い場合の比較.
1.14 takayama 1253: @example
1254: [4618] cputime(1)$
1255: [4619] gtt_ekn3.expectation(Marginal=[[1950,2550,5295],[1350,1785,6660]],
1256: P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]|bs=1)$
1257: 4.912sec(4.914sec)
1258: [4621] V2=gtt_ekn3.expectation(Marginal=[[1950,2550,5295],[1350,1785,6660]],
1259: P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]])$
1260: 6.752sec(6.756sec)
1261: @end example
1262:
1263:
1264: @comment --- @example〜@end example は実行例の表示 ---
1.15 takayama 1265: 例: 分散計算する場合.
1266: 分散計算はかえって遅くなる場合が多いので注意.
1267: 下記の例での bsplit_x=1 option は
1268: debug windows を開くのでさらに遅くなる.
1269: gtt_ekn3.test_bs_dist(); でもテストできる.
1.14 takayama 1270: @example
1.15 takayama 1271: [3669] C=4$ gtt_ekn3.init_bsplit(|minsize=16,levelmax=C+1)$ gtt_ekn3.init_dm_bsplit(|bsplit_x=1)$
1.14 takayama 1272: [3670] [3671] [3672] gtt_ekn3.setup_dm_bsplit(C);
1273: [0,0]
1274: [3673] gtt_ekn3.assert2(10|bs=1)$
1275: @end example
1276:
1277: @comment --- 参照(リンク)を書く ---
1278: @table @t
1279: @item 参照
1280: @ref{gtt_ekn3.gmvector}
1281: @ref{gtt_ekn3.expectation}
1282: @ref{gtt_ekn3.assert1}
1283: @ref{gtt_ekn3.assert2}
1284: @end table
1285:
1286: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1287: @noindent
1288: ChangeLog
1289: @itemize @bullet
1290: @item
1291: 関連ファイルは
1292: gtt_ekn3/mfac_bs.rr
1293: gtt_ekn3/dm_bsplit.rr
1294: @end itemize
1295:
1.4 takayama 1296:
1.1 takayama 1297:
1.6 takayama 1298: @comment --- おまじない ---
1.1 takayama 1299: @node Index,,, Top
1300: @unnumbered Index
1301: @printindex fn
1302: @printindex cp
1303: @iftex
1304: @vfill @eject
1305: @end iftex
1306: @summarycontents
1307: @contents
1308: @bye
1.6 takayama 1309: @comment --- おまじない終り ---
1.1 takayama 1310:
1311:
1.6 takayama 1312: @comment テンプレート. start_of_template.
1.5 takayama 1313: @comment **********************************************************
1.6 takayama 1314: @comment --- ◯◯◯◯ の説明
1315: @comment --- 個々の関数の説明の開始 ---
1316: @comment --- section 名を正確に ---
1.19 takayama 1317: @node gtt_ekn3.hoge,,, 超幾何関数E(k,n)
1318: @subsection @code{gtt_ekn3.hoge}
1.6 takayama 1319: @comment --- 索引用キーワード
1.19 takayama 1320: @findex gtt_ekn3.hoge
1.5 takayama 1321:
1322: @table @t
1.19 takayama 1323: @item gtt_ekn3.hoge(@var{i})
1.5 takayama 1324: ::
1325: @end table
1326:
1.6 takayama 1327: @comment --- 引数の簡単な説明 --- 以下まだ書いてない.
1.5 takayama 1328: @table @var
1329: @item i hage
1330: @item return
1331: @end table
1332:
1.6 takayama 1333: @comment --- ここで関数の詳しい説明 ---
1334: @comment --- @itemize〜@end itemize は箇条書き ---
1335: @comment --- @bullet は黒点付き ---
1.5 takayama 1336: @itemize @bullet
1.6 takayama 1337: @item 説明.
1.5 takayama 1338: @end itemize
1339:
1.6 takayama 1340: @comment --- @example〜@end example は実行例の表示 ---
1341: 例:
1.5 takayama 1342: @example
1.19 takayama 1343: [2221] gtt_ekn3.hoge([[1,4],[2,3]]);
1.5 takayama 1344: @end example
1345:
1346:
1.6 takayama 1347: @comment --- 参照(リンク)を書く ---
1.5 takayama 1348: @table @t
1.6 takayama 1349: @item 参照
1.19 takayama 1350: @ref{gtt_ekn3.nc}
1351: @ref{gtt_ekn3.gmvector}
1.5 takayama 1352: @end table
1353:
1.6 takayama 1354: @comment --- ChangeLog を書く. ソースコードの位置. 変更日時 など CVSサーバを見るため
1.5 takayama 1355: @noindent
1356: ChangeLog
1357: @itemize @bullet
1358: @item
1359: @end itemize
1360: @comment end_of_template
1361:
1362:
1.6 takayama 1363: // 2 x m 分割表において似た機能を有する関数の利用例を参考までに記載する;
1364: // 正規化定数とその微分関連.
1365: // その1.
1.1 takayama 1366: [3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
1367: [-4,[-4,-3],-1]
1368: [3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
1369: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
1370: [ 1 1 1 ]
1371: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
1372: [4483/124416,[[353/7776,1961/15552,185/1728],[553/20736,1261/15552,1001/13824]]]
1.6 takayama 1373: // 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],[d_21 Z, d_22 Z, d_23 Z]]] の値.
1.1 takayama 1374:
1.6 takayama 1375: // その2.
1.1 takayama 1376: [3079] tk_fd.log_ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
1377: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
1378: [ 1 1 1 ]
1379: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
1380: [-3.32333832422461674639485797719209322217260539267246045320,
1381: [[1.25987062235110417131385233102832924380994869507026544724,3.49944233772027660049074280615659156814633058219942003122,2.97122462636627258532232879768012491635065804149007361142],
1382: [0.740129377648895828686147668971670756190051304929734552754,2.25027883113986169975462859692170421592683470890028998438,2.00959179121124247155922373410662502788311398616997546285]]]
1.6 takayama 1383: // 戻値は [log(Z),
1.1 takayama 1384: // [[d_11 log(Z), d_12 log(Z), d_13 log(Z)],
1385: // [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]]
1.6 takayama 1386: // の近似値.
1.1 takayama 1387:
1.6 takayama 1388: // その3.
1.1 takayama 1389: [3082] fd_hessian2(A[0],A[1],A[2],[1/2,1/3]);
1390: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
1391: [4483/124416,[ 1961/15552 185/1728 ],
1392: [ 79/288 259/864 ]
1393: [ 259/864 47/288 ]]
1.6 takayama 1394: // 戻値は [F=F_D, gradient(F), Hessian(F)]
1.1 takayama 1395:
1.6 takayama 1396: // 参考.
1397: // ygahvec で巾関数分の調整. 独立した関数はないようだ.
1.1 takayama 1398:
1399: //-----------------------------------------------------------------------
1.6 takayama 1400: // 2 x m 分割表において似た機能を有する関数の利用例を参考までに記載する;
1401: // 期待値関連.
1.1 takayama 1402: [3079] A=tk_fd.marginal2abc([4,5],[2,4,3]);
1403: [-4,[-4,-3],-1]
1404: [3080] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
1405: RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
1406: [ 1 1 1 ]
1407: Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
1408: [[5648/4483,7844/4483,4440/4483],
1409: [3318/4483,10088/4483,9009/4483]]
1.6 takayama 1410: // 各セルの期待値.
1.1 takayama 1411:
1412: //-----------------------------------------------------------------------
1.6 takayama 1413: // ot_hgm_ahg.rr の例. 実験的なため module 化されていない.
1.1 takayama 1414: [3237] import("ot_hgm_ahg.rr");
1.6 takayama 1415: // 2 x 2 分割表.
1.1 takayama 1416: [3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],
1417: [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
1418: oohg_native=0, oohg_curl=1
1419: [1376777025/625400597,1750225960/625400597,2375626557/625400597,3252978816/625400597]
1.6 takayama 1420: // 2 x 2 分割表の期待値.
1.1 takayama 1421:
1.6 takayama 1422: // 2 x 3 分割表.
1.1 takayama 1423: [3238] hgm_ahg_expected_values_contiguity(
1424: [[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
1425: [5,2,4,3],[1,1/2,1/3,1,1,1],[x1,x2,x3,x4,x5,x6]|geometric=1);
1426: [5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
1.6 takayama 1427: // 2 x 3 分割表の期待値. 上と同じ問題.
1.1 takayama 1428:
1429: /*
1.6 takayama 1430: dojo, p.221. 成績3以下の生徒は集めてひとつに.
1.1 takayama 1431: 2 1 1
1432: 8 3 3
1433: 0 2 6
1434:
1435: row sum: 4,14,8
1436: column sum: 10,6,10
1.6 takayama 1437: 0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
1.1 takayama 1438: */
1.6 takayama 1439: // 3 x 3 分割表. 構造的0が一つ.
1.1 takayama 1440:
1441: A=[[0,0,0,1,1,1, 0,0],
1442: [0,0,0,0,0,0, 1,1],
1443: [1,0,0,1,0,0, 0,0],
1444: [0,1,0,0,1,0, 1,0],
1445: [0,0,1,0,0,1, 0,1]];
1446: B=[14,8,10,6,10];
1447: 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);
1448:
1.6 takayama 1449: // 答.
1.1 takayama 1450: [14449864949304/9556267369631,10262588586540/9556267369631,13512615942680/9556267369631,
1451: 81112808747006/9556267369631,21816297744346/9556267369631,30858636683482/9556267369631,
1452: 25258717886900/9556267369631,51191421070148/9556267369631]
1453:
1454:
1455: /*
1.6 takayama 1456: 上のデータで 0 を 1 に変更.
1.1 takayama 1457: 2 1 1
1458: 8 3 3
1459: 1 2 6
1460:
1461: row sum: 4,14,9
1462: column sum: 11,6,10
1463: */
1.6 takayama 1464: // 3 x 3 分割表.
1.1 takayama 1465: A=[[0,0,0,1,1,1,0,0,0],
1466: [0,0,0,0,0,0,1,1,1],
1467: [1,0,0,1,0,0,1,0,0],
1468: [0,1,0,0,1,0,0,1,0],
1469: [0,0,1,0,0,1,0,0,1]];
1470: B=[14,9,11,6,10];
1471: 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);
1472:
1.6 takayama 1473: // 期待値, 答.
1.1 takayama 1474: [207017568232262040/147000422096729819,163140751505489940/147000422096729819,217843368649167296/147000422096729819,
1475: 1185482401011137878/147000422096729819,358095302885438604/147000422096729819,514428205457640984/147000422096729819,
1476: 224504673820628091/147000422096729819,360766478189450370/147000422096729819]
1477:
1.6 takayama 1478: // Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
1479: // まだ書いてない.
1.1 takayama 1480:
1481:
1.6 takayama 1482: 4. x_ij は [GM2016] の１章で,
1483: たとえば 3x3 の時 [[1,1,1],[x_11,x_12,1],[x_21,x_22,1]]
1484: となっているが, [GM2016] の Prop 7.1 の対応では,
1485: p = [[1,x_11,x_12],[1,x_21,x_22],[1,1,1]] となっているので注意.
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gtt_ekn-ja.texi CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / packages / doc / gtt_ekn