Python で Fisher's exact test をする（Mac の場合）

前に書いた Fisher's exact test(Python) では諦めていたのだけど...

Python で Fisher's exact test を行えるようにした話

まず，目的の FORTRAN プログラムソースをダウンロードする
http://netlib.org/toms/ の中にある file: 643.gz
"unordered rxc contingency tables, Fisher's exact test"

展開してできる 643 というファイルを開き，
　　　1965 行以降のテストデータを削除する。
　　　1~31 行のメインプログラムを以下のサブルーチンプログラムで置き換える。
　　　名前を fisher_lib.f90 として保存する。

C      ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488.
c                                   Driving program to read and
c                                   process the data file
      subroutine fisher_sub(nrow, ncol, expect, percnt, emin, table,
     &  prt, pre)
      implicit none

      integer    i, j, ncol, nrow
      double precision emin, expect, percnt, pre, prt, table(*)
      intent(in) nrow, ncol, expect, percnt, emin, table
      intent(out) prt, pre
c
      external   fexact
c                                  output TABLE
c      write (*,99998)
c      do 20  i=1, nrow
c         write (*,99997) (table(i+(j-1)*nrow),j=1,ncol)
c   20 continue

      call fexact (nrow, ncol, table, nrow, expect, percnt, emin, prt,
     &             pre)

c99997 format (1x, 10f7.0)
c99998 format (/, 2x, 'The contingency table for this problem is:')
      return
      end

ターミナルで，以下を実行する（fisher_lib*.so ファイルができる）

$ f2py -m fisher_lib -c fisher_lib.f90
 
使用法は，fisher_lib をインポートして，fisher_lib.fisher_sub 関数を呼び出すだけ。
下の fisher() 関数のようなシュガーコートを介して呼び出すと楽ちん。

import numpy as np
import fisher_lib

def fisher(table, expect=0, percnt=0, emin=0):
    nrow, ncol = np.array(table).shape
    prt, pre = fisher_lib.fisher_sub(nrow, ncol, expect, percnt, emin, table)
    return pre

集計表は，以下のように二重リストで与えればよい（numpy.array() でもよい）。

>>> import numpy as np
>>> import fisher_lib

>>> def fisher(table, expect=0, percnt=0, emin=0):
...     nrow, ncol = np.array(table).shape
...     prt, pre = fisher_lib.fisher_sub(nrow, ncol, expect, percnt, emin, table)
...     return pre
... 
>>> table = [[3, 5, 2, 5], [4, 1, 8, 3], [2, 1, 3, 4]]
>>> p = fisher(table)
>>> print(p)
0.26187809176685495

>>> p = fisher(np.array([[12, 2], [5, 11]]))
>>> print(p)
0.003913481855563419

R の fisher.test() で結果が正しいことを確かめる。

> options(digits=15)
> fisher.test(matrix(c(3, 5, 2, 5,  4, 1, 8, 3,  2, 1, 3, 4), byrow = TRUE, ncol = 4))$p.value
[1] 0.261878091766855
> fisher.test(matrix(c(12, 2, 5, 11), 2))$p.value
[1] 0.00391348185556345

 

Python の変な仕様？

import scipy as sp

x = sp.reshape(sp.arange(36), (6, 6))

で作成される行列 x

array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23],
       [24, 25, 26, 27, 28, 29],
       [30, 31, 32, 33, 34, 35]])

この行列の小行列

array([[ 6,  8, 10],
       [18, 20, 22],
       [30, 32, 34]])

を作成しようとして，

x[[1,3,5], [0,2,4]]

なんてことをやると

array([ 6, 20, 34])

が得られて，なんだこりゃ？となってしまう。

R なら，これで問題ないのになあと（0 オリジンと 1 オリジンの違いはあるものの）。

> x = matrix(0:35, 6, 6, byrow=TRUE)
> x
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    1    2    3    4    5
[2,]    6    7    8    9   10   11
[3,]   12   13   14   15   16   17
[4,]   18   19   20   21   22   23
[5,]   24   25   26   27   28   29
[6,]   30   31   32   33   34   35

> x[c(2,4,6), c(1,3,5)]
     [,1] [,2] [,3]
[1,]    6    8   10
[2,]   18   20   22
[3,]   30   32   34

Python では，

x[[1,3,5], :][:, [0,2,4]]

でようやっと

array([[ 6,  8, 10],
       [18, 20, 22],
       [30, 32, 34]])

x[[1,3,5], [0,2,4]] で得られた array([ 6, 20, 34]) は，目的とした小行列の対角成分だなあ。

R では，x[c(2,4,6), ][, c(1,3,5)] は x[c(2,4,6), c(1,3,5)] と同じなんだよ。