【Perl】自由度nのカイ二乗値の計算
2013-12-21-1
[Algorithm][Programming]
よく使うんだけど、その度にいちいちコード書くのが面倒なのでコピペ用にここにメモ。
適合度の検定。理論値とのズレを検定。
■ソース (chisq.pl)
■実行結果:
独立性の検定。考えられるテーブルすべての結果を出す。
■ソース (chisqin.pl)
■実行結果:「*」のつく数字は合計を表す。
自由度1のとき専用の独立性の検定プログラム。
■ソース (chisqi.pl)
■実行結果:
- 東京大学教養学部統計学教室 / 統計学入門 (基礎統計学)
- Wikipedia:カイ二乗検定
- カイ2乗分布,カイ2乗検定
http://www.geisya.or.jp/~mwm48961/statistics/kai2.htm
- カイ2乗分布表 chi-square distribution
http://www.biwako.shiga-u.ac.jp/sensei/mnaka/ut/chi2disttab.html
(自由度1のχ2乗分布)
- [を] カイ二乗値で単語間の関連の強さを調べる[2007-09-19-1]
追記140108: ソースと実行例を追加。適合度の検定、および自由度1をnに。今までのソースは独立性の検定用。
追記140204: 自由度2以上にも対応した独立性検定のソースと実行例を追加。
適合度の検定
適合度の検定。理論値とのズレを検定。
■ソース (chisq.pl)
#!/usr/bin/perl use strict; use warnings; my $num = @ARGV; my @x = @ARGV[0..$num/2-1]; # 期待値・理論値 my @y = @ARGV[$num/2..$#ARGV]; # 実測値 my $val = chisq(\@x, \@y); print "$val\n"; sub chisq { my ($x_ref, $y_ref) = @_; my $rv = 0; for (my $i = 0; $i < @$x_ref; $i++) { $rv += ($y_ref->[$i]-$x_ref->[$i])**2/$x_ref->[$i]; } return $rv; }
■実行結果:
% chisq.pl 50 50 45 55 1 % chisq.pl 50 50 47 53 0.36 % chisq.pl 40 20 10 30 37 25 12 26 2.40833333333333
独立性の検定
独立性の検定。考えられるテーブルすべての結果を出す。
■ソース (chisqin.pl)
#!/usr/bin/perl use strict; use warnings; my @nums = @ARGV; for (my $x_max = 2; $x_max < @nums; $x_max++) { next if @nums % $x_max; my $y_max = @nums / $x_max; my ($u2, $m_ref) = chisqin(\@nums, $x_max, $y_max); print "table = $x_max x $y_max, chisquare = $u2\n"; for (my $j = 0; $j <= $y_max; $j++) { for (my $i = 0; $i <= $x_max; $i++) { print $m_ref->[$i][$j].(($i==$x_max or $j==$y_max)?"*":"")."\t"; } print "\n"; } print "\n"; } sub chisqin { my ($m_ref, $x_max, $y_max) = @_; my @mat; for (my $j = 0; $j < $y_max; $j++) { for (my $i = 0; $i < $x_max; $i++) { $mat[$i][$j] = $m_ref->[$x_max * $j + $i]; $mat[$x_max][$j] += $mat[$i][$j]; # x_sum $mat[$i][$y_max] += $mat[$i][$j]; # y_sum $mat[$x_max][$y_max] += $mat[$i][$j]; # all_sum } } my $u2; for (my $j = 0; $j < $y_max; $j++) { for (my $i = 0; $i < $x_max; $i++) { my $z = $mat[$x_max][$j] * $mat[$i][$y_max] / $mat[$x_max][$y_max]; $u2 += ($mat[$i][$j]-$z)**2/$z; } } return $u2, \@mat; }
■実行結果:「*」のつく数字は合計を表す。
% ./chisqin.pl 1 99 5 195 table = 2 x 2, chisquare = 0.76530612244898 1 99 100* 5 195 200* 6* 294* 300* % ./chisqin.pl 50 50 47 53 table = 2 x 2, chisquare = 0.180162145931338 50 50 100* 47 53 100* 97* 103* 200* % ./chisqin.pl 4 2 3 8 4 6 6 3 6 table = 3 x 3, chisquare = 0.186666666666667 4 2 3 9* 8 4 6 18* 6 3 6 15* 18* 9* 15* 42* % ./chisqin.pl 4 2 3 8 4 6 6 3 6 1 1 1 table = 2 x 6, chisquare = 7.74524582560296 4 2 6* 3 8 11* 4 6 10* 6 3 9* 6 1 7* 1 1 2* 24* 21* 45* table = 3 x 4, chisquare = 0.430263157894737 4 2 3 9* 8 4 6 18* 6 3 6 15* 1 1 1 3* 19* 10* 16* 45* table = 4 x 3, chisquare = 12.5141284584009 4 2 3 8 17* 4 6 6 3 19* 6 1 1 1 9* 14* 9* 10* 12* 45* table = 6 x 2, chisquare = 11.0582010582011 4 2 3 8 4 6 27* 6 3 6 1 1 1 18* 10* 5* 9* 9* 5* 7* 45*
自由度1の独立性の検定
自由度1のとき専用の独立性の検定プログラム。
■ソース (chisqi.pl)
#!/usr/bin/perl use strict; use warnings; my ($o11, $o12, $o21, $o22) = @ARGV; # my ($n1, $o11, $n2, $o21) = @ARGV; # my $o12 = $n1 - $o11; # my $o22 = $n2 - $o21; my $val = chisqi($o11, $o12, $o21, $o22); print "$val\n"; sub chisqi { my ($o11, $o12, $o21, $o22) = @_; return ($o11 + $o12 + $o21 + $o22)*(($o11*$o22-$o12*$o21)**2)/ (($o11+$o12)*($o11+$o21)*($o12+$o22)*($o21+$o22)); }
■実行結果:
% ./chisqi.pl 1 99 5 195 0.76530612244898 % ./chisqi.pl 50 50 47 53 0.180162145931338 % ./chisqi.pl 50 50 55 45 0.50125313283208
参考
- 東京大学教養学部統計学教室 / 統計学入門 (基礎統計学)
- Wikipedia:カイ二乗検定
- カイ2乗分布,カイ2乗検定
http://www.geisya.or.jp/~mwm48961/statistics/kai2.htm
- カイ2乗分布表 chi-square distribution
http://www.biwako.shiga-u.ac.jp/sensei/mnaka/ut/chi2disttab.html
上側確率 | 0.050 | 0.025 | 0.010 | 0.005 |
χ2 | 3.84146 | 5.02389 | 6.63490 | 7.87944 |
- [を] カイ二乗値で単語間の関連の強さを調べる[2007-09-19-1]
更新履歴
追記140108: ソースと実行例を追加。適合度の検定、および自由度1をnに。今までのソースは独立性の検定用。
追記140204: 自由度2以上にも対応した独立性検定のソースと実行例を追加。