MATLABによるフィッシャーの正確確率検定

MATLAB analysis Bio
  • Fisher's exact test.
  • fisher_exact()関数（霧笛作成）を使用する。



引数

data  : データ値のベクトル(2 * 2)
          data = [a b; c d]
alpha : 有意水準
both  : 0 片側検定
        1 両側検定



返り値

h     : 0 帰無仮説は採択される → 2要因は独立でないとはいえない
        1 帰無仮説は棄却される → 2要因は独立ではない
p     : 生起確率



プログラム例

% ---------------------------------------------------------------
% フィッシャーの正確確率検定
% Fisher's exact test.
%
% Web教材「統計学自習ノート」より、
% フィッシャーの正確確率検定（直接確率）「例題」の計算。
% http://aoki2.si.gunma-u.ac.jp/lecture/Cross/Fisher.html
% ---------------------------------------------------------------
function exact_test()

    % p 生起確率
    % h : 0 帰無仮説は採択される → 2要因は独立でないとはいえない
    %   : 1 帰無仮説は棄却される → 2要因は独立ではない

    data = [13 4; 6 14];
    alpha = 0.01;           % 有意水準

    % 片側検定
    [h,p] = fisher_exact(data,alpha,0)    
    % 両側検定
    [h,p] = fisher_exact(data,alpha,1)
end



プログラム例−実行結果

>> exact_test

p =
   0.0059

h =
    1

p =
   0.0081

h =
    1



fisher_exact()関数の実装

  • fisher_exact.m
% ---------------------------------------------------------------
% フィッシャーの正確確率検定
% Fisher's exact test.
% 引数:
%   data  : データ値のベクトル（2 * 2）
%           data = [a b; c d]
%   alpha : 有意水準
%   both  : 0 片側検定
%           1 両側検定
% 返り値:
%   h     : 0 帰無仮説は採択される → 2要因は独立でないとはいえない
%           1 帰無仮説は棄却される → 2要因は独立ではない
%   allp  : 生起確率
% ---------------------------------------------------------------
function [h,allp] = fisher_exact(data, alpha, both)

    h = 0;
    allp = 0.0;
    % +---------------+--------+
    % |    a     b    | sum_ab |
    % |    c     d    | sum_cd |
    % +---------------+--------+
    % | sum_ac sum_bd |        |
    % +---------------+--------+
    a = data(1,1);
    b = data(1,2);
    c = data(2,1);
    d = data(2,2);
    sum_ab = a + b;
    sum_cd = c + d;
    sum_ac = a + c;
    sum_bd = b + d;

    % 実際の観測値の正確な生起確率を得る
    p0 = get_exact_p(a,b,c,d);
    s0 = a*d - b*c;

    % aは 0〜小さい方の合計値まで推移
    index = 0;
    if(sum_ab <= sum_ac)
        for a=0:sum_ab
            index = index + 1;
            b = sum_ab - a;
            c = sum_ac - a;
            d = sum_cd - c;
            % 正確な生起確率を得る
            p(index) = get_exact_p(a,b,c,d);
            s(index) = a*d - b*c;
        end
    else
        for a=0:sum_ac
            index = index + 1;
            b = sum_ab - a;
            c = sum_ac - a;
            d = sum_cd - c;
            % 正確な生起確率を得る
            p(index) = get_exact_p(a,b,c,d);
            s(index) = a*d - b*c;
        end
    end
    number = index;

    % 実際の観測値よりも極端な場合の生起確率の合計を求める
    for index=1:number
        % 片側検定では（sとs0が同符号でかつ | s | ≧ | s0 |）の生起確率の合計
        if(both == 0)
            sflag = 0;
            % 同符号なら sflag = 1
            if((s0 >= 0 && s(index) >= 0) || (s0 < 0 && s(index) < 0))
                sflag = 1;
            end
            if((sflag == 1) && (abs(s(index)) >= abs(s0)))
                allp = allp + p(index);
            end
        % 両側検定では（| s | ≧ | s0 |）の生起確率の合計
        else
            if(abs(s(index)) >= abs(s0))
                allp = allp + p(index);
            end
        end
    end

    % 生起確率が有意水準以下
    if(allp <= alpha)
        h = 1
    end
end

% ---------------------------------------------------------------
% 正確な生起確率を得る
% 引数:
%   a, b, c, d  : データ値（2 * 2）
% 返り値:
%   p           : 生起確率
% ---------------------------------------------------------------
function p = get_exact_p(a,b,c,d)

    % prod(1:n)は nの階乗の計算
    n = a + b + c + d;
    d1 = prod(1:a+b) * prod(1:c+d) * prod(1:a+c) * prod(1:b+d);
    d2 = prod(1:n) * prod(1:a) * prod(1:b) * prod(1:c) * prod(1:d);
    p = d1 / d2;
end