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