Last active
August 29, 2015 14:02
-
-
Save snipsnipsnip/051c277dbed9ff2ee12e to your computer and use it in GitHub Desktop.
J言語でF分布とt分布とカイ二乗分布の累積分布関数 (CDF of χ² distribution, Student's t-distribution and F-distribution in J programming language)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
後で探したら H. について日本語の解説と既にCDFの実装があった。日本語で書かれたリファレンスまで。すごく有り難い。