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