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