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