古い記事
ランダムジャンプ
新しい記事
よく使うんだけど、その度にいちいちコード書くのが面倒なのでコピペ用にここにメモ。

適合度の検定


適合度の検定。理論値とのズレを検定。

■ソース (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.0500.0250.0100.005
χ23.841465.023896.634907.87944
(自由度1のχ2乗分布)

- [を] カイ二乗値で単語間の関連の強さを調べる[2007-09-19-1]

更新履歴


追記140108: ソースと実行例を追加。適合度の検定、および自由度1をnに。今までのソースは独立性の検定用。

追記140204: 自由度2以上にも対応した独立性検定のソースと実行例を追加。