J言語でF分布とt分布とカイ二乗分布の累積分布関数 (CDF of χ² distribution, Student's t-distribution and F-distribution in J programming language)

cdf.ijs

NB. ガンマ関数。gamma x = Γ(x)
gamma =: ! @ <:

NB. ベータ関数。 x B y = B(x, y)
beta =: * & gamma % gamma @ +

NB. 不完全ベータ関数。a b incomplete_beta x = B(x;a,b)
NB. H.で一般化超幾何級数を求められるらしい cf. http://www.jsoftware.com/jwiki/Vocabulary/hcapdot
incomplete_beta =: dyad define
  'm n' =. x
  (y ^ m) * (((m , 1 - n) H. >: m) y) % m
)

NB. 正規化不完全ベータ関数。a b regularized_incomplete_beta x = I(x;a,b)
regularized_incomplete_beta =: dyad define
  (x incomplete_beta y) % beta/ x
)

NB. 不完全ガンマ関数。s incomplete_gamma x = γ(s,x)
incomplete_gamma =: dyad define
  (y ^ x) * ((x H. >: x) -y) % x
)

NB. 正規化不完全ガンマ関数。s incomplete_gamma x = Γ(s,x)
regularized_incomplete_gamma =: dyad define
  (x incomplete_gamma y) % gamma x
)

NB. χ²分布の累積分布関数。 k chi2_dist_cdf x = chi2_CDF(x;k)
chi2_dist_cdf =: dyad define
  (-: x) regularized_incomplete_gamma -: y
)

NB. t分布の累積分布関数。 ν t_dist_cdf x = t_CDF(x;ν)
t_dist_cdf =: dyad define
  1 - -: (1r2 , -: x) regularized_incomplete_beta x % x + *: y
)

NB. F分布の累積分布関数。d1 d2 f_dist_cdf x = F_CDF(x;d1,d2)
f_dist_cdf =: dyad define
  'd1 d2' =. x
  d1x =. d1 * y
  (-: x) regularized_incomplete_beta d1x % d1x + d2
)

error.ijs

   NB. 拾った信頼率95%のF分布表。行は群内の自由度で10,20...,90、列は群間の自由度で1,2,...,5
   tf =. 9 5 $	4.96	4.10	3.71	3.48	3.33 	4.35	3.49	3.10	2.87	2.71 	4.17	3.32	2.92	2.69	2.53 	4.08	3.23	2.84	2.61	2.45 	4.03	3.18	2.79	2.56	2.40 	4.00	3.15	2.76	2.53	2.37 	3.98	3.13	2.74	2.50	2.35 	3.96	3.11	2.72	2.49	2.33 	3.95	3.10	2.71	2.47	2.32

   NB. テスト用に生成した引数リスト
   (45 2 $ , (10 * >: i. 9) ,~"0/ (>: i. 5)) ,"1 ,. , tf
1 10 4.96
2 10  4.1
3 10 3.71
4 10 3.48
5 10 3.33
1 20 4.35
2 20 3.49
3 20  3.1
4 20 2.87
5 20 2.71
1 30 4.17
2 30 3.32
3 30 2.92
4 30 2.69
5 30 2.53
1 40 4.08
2 40 3.23
3 40 2.84
4 40 2.61
5 40 2.45
1 50 4.03
2 50 3.18
3 50 2.79
4 50 2.56
5 50  2.4
1 60    4
2 60 3.15
3 60 2.76
4 60 2.53
5 60 2.37
1 70 3.98
2 70 3.13
3 70 2.74
4 70  2.5
5 70 2.35
1 80 3.96
2 80 3.11
3 80 2.72
4 80 2.49
5 80 2.33
1 90 3.95
2 90  3.1
3 90 2.71
4 90 2.47
5 90 2.32
   NB. 表を全部CDFに通して誤差を出してみた
   0j6 ": 9 5 $ 0.95 - , (45 2 $ , (10 * >: i. 9) ,~"0/ (>: i. 5)) f_dist_cdf"1 ,. , tf
 0.000088  0.000078 _0.000057 _0.000073 _0.000169
 0.000030  0.000105 _0.000075 _0.000214  0.000055
 0.000023 _0.000170  0.000120 _0.000023  0.000253
 0.000127  0.000074 _0.000070 _0.000269 _0.000041
 0.000118  0.000116  0.000000 _0.000197  0.000033
 0.000033  0.000019 _0.000114 _0.000343 _0.000143
_0.000062 _0.000107 _0.000269  0.000195 _0.000371
 0.000010  0.000036 _0.000075 _0.000305 _0.000110
_0.000089 _0.000108 _0.000257  0.000221 _0.000374

後で探したら H. について日本語の解説と既にCDFの実装があった。日本語で書かれたリファレンスまで。すごく有り難い。