Fisherの正確確率検定
nxm分割表の正確確率検定
先日、HWE検定の正確確率検定版について記載した(こちら)。ついでに、いわゆるFisherの正確確率検定のnxm分割表用のソースも載せる。
解説は、群馬大青木先生のこちらのページがよくわかる。
public static double Fishernxm2(int[][] da){
/*
* nxm分割表のFisher正確確率検定
*/
/*
* 関数limitCubeは分割表から、
* 各セルの取りうる整数値の範囲を限定する
* 各セルが取りうる値の範囲は、周辺度数の制約と
* 各セルが0以上である、という制約とから定まる
* 各セルの取りうる範囲は線形計画法的な多次元立体（斜面を持つ）
* で囲まれるが
* ここでは、多次元立方体(斜面なし）の範囲を定め、
* その範囲の点をしらみつぶしにする
*/
int[][][] limit = limitCube(da);
/*
* 各セルの取りうる値の数を格納
*/
int[][] limitnum = new int[limit.length-1][limit[0].length-1];
for(int i=0;i<limitnum.length;i++){
for(int j=0;j<limitnum[i].length;j++){
limitnum[i][j]=limit[i][j][1]-limit[i][j][0]+1;
}
}
/*
* logfact関数の値は繰り返し使うので、初めに
* 一括して出力しておく(logfactTable[][])
*/
double[][][] logfactTable = new double[limit.length][limit[0].length][];
for(int i=0;i<logfactTable.length;i++){
for(int j=0;j<logfactTable[i].length;j++){
logfactTable[i][j] = new double[limit[i][j][1]-limit[i][j][0]+1];
logfactTable[i][j][0] = logfact(limit[i][j][0]);
for(int k=1;k<logfactTable[i][j].length;k++){
logfactTable[i][j][k] = logfactTable[i][j][k-1]+Math.log(limit[i][j][0]+k);
}
}
}
/*
* 周辺度数p と q
* 作成する分割表を得る確率は、分割表のセルの値と
* 周辺度数とで決まる
* 周辺度数による分は繰り返し計算に使うので、
* constantとして、あらかじめ求めておく
*/
int[] p= new int[da.length];
int[] q= new int[da[0].length];
int sum = 0;
for(int i=0;i<da.length;i++){
for(int j=0;j<da[i].length;j++){
p[i]+=da[i][j];
q[j]+=da[i][j];
sum += da[i][j];
}
}
double constant = 0;
for(int i=0;i<p.length;i++){
constant += logfact(p[i]);
}
for(int j=0;j<q.length;j++){
constant += logfact(q[j]);
}
constant -= logfact(sum);
/*
* 観測データを得る確率は、取りうる分割表を得る確率と比較し、
* それ以下の確率の場合を累積してP値を得る。
* value0は、その観測データを得る確率を求めるための値
* constantと併せて最終的に確率値を得る
*/
double value0 = 0;
for(int i=0;i<da.length;i++){
for(int j=0;j<da[i].length;j++){
value0 += logfact(da[i][j]);
}
}
/*
* nxm分割表は(n-1)x(m-1)の自由度を持つ
* (n-1)x(m-1)の値のとりうる範囲を動かすための変数nshin
* (n-1)x(m-1)桁で、各桁が異なるk進数となるような設定
*/
int[][] nshin = new int[da.length-1][da[0].length-1];
for(int i=0;i<nshin.length;i++){
for(int j=0;j<nshin[i].length;j++){
nshin[i][j] = 0;
}
}
/*
* すべての場合を網羅する
*/
boolean returned = true;
double exactP = 0;
while(returned){
/*
* nshinから、分割表を作る
*/
int[][] tmpvar = createFromNshin(limit,nshin,p,q);
/*
* できた分割表の全セルが0以上かをチェック
*/
boolean positiveCheck = deltaChecker0(tmpvar);
/*
* 全セル0以上の場合は確率を求め、
* 観測データのそれ以下なら加算
*/
if(positiveCheck){
int[][] nshinMarginal = new int[p.length][q.length];
for(int i=0;i<p.length-1;i++){
for(int j=0;j<q.length-1;j++){
nshinMarginal[i][j]=nshin[i][j];
}
}
for(int i=0;i<p.length-1;i++){
nshinMarginal[i][q.length-1]=tmpvar[i][q.length-1]-limit[i][q.length-1][0];
}
for(int i=0;i<q.length;i++){
nshinMarginal[p.length-1][i]=tmpvar[p.length-1][i]-limit[p.length-1][i][0];
}
double tmpdenom = 0;
for(int i=0;i<p.length;i++){
for(int j=0;j<q.length;j++){
tmpdenom += logfactTable[i][j][nshinMarginal[i][j]];
}
}
if(tmpdenom>=value0){
double tmplogfact=constant-tmpdenom;
double tmpexp = Math.exp(tmplogfact);
exactP += tmpexp;
}
}
/*
* nshinをひとつ大きくする
*/
nshin = nshinUpper(nshin,limitnum);
/*
* nshinがぐるっと回って、初めに戻ったら
* 網羅ループ終了
*/
boolean zerocheck = deltaCheckerAllZero(nshin);
if(zerocheck){
returned =false;
}
}
return exactP;
}
public static double logfact(int a){
double ret =0;
for(int i=1;i<=a;i++){
ret += Math.log(i);
}
return ret;
}
public static int[][] createFromNshin(int[][][] limit,int[][] n, int[] p, int[] q){
int[][] ret = new int[limit.length][limit[0].length];
int[] psum = new int[n.length];
int[] qsum = new int[n[0].length];
for(int i=0;i<n.length;i++){
for(int j=0;j<n[i].length;j++){
ret[i][j]=limit[i][j][0]+n[i][j];
psum[i] += ret[i][j];
qsum[j] += ret[i][j];
}
}
for(int i=0;i<n.length;i++){
ret[i][ret[i].length-1]=p[i]-psum[i];
}
for(int i=0;i<n[0].length;i++){
ret[ret.length-1][i]=q[i]-qsum[i];
}
ret[ret.length-1][ret[0].length-1]=p[ret.length-1];
for(int i=0;i<q.length-1;i++){
ret[ret.length-1][ret[0].length-1]-=ret[ret.length-1][i];
}
return ret;
}
public static int[][][] limitCube(int[][] d){
int[][][] ret = new int[d.length][d[0].length][2];
int[] p = new int[d.length];
int[] q = new int[d[0].length];
int sum = 0;
for(int i=0;i<p.length;i++){
for(int j=0;j<q.length;j++){
p[i] += d[i][j];
q[j] += d[i][j];
sum += d[i][j];
}
}
for(int i=0;i<p.length;i++){
for(int j=0;j<q.length;j++){
ret[i][j][0] = 0;
ret[i][j][1] = Math.min(p[i], q[j]);
}
}
boolean identical = false;
while(!identical){
int[][][] ret2 = repeatLimitCube(ret,p,q,sum);
identical = identicalInt3(ret2,ret);
for(int i=0;i<ret.length;i++){
ret[i] = MiscUtil.DeepCopyInt2(ret2[i]);
}
/*
if(identical){
for(int i=0;i<ret.length;i++){
ret[i] = MiscUtil.DeepCopyInt2(ret2[i]);
}
}
*/
}
return ret;
}
public static int[][][] repeatLimitCube(int[][][] d,int[] p,int[] q,int sum){
int[][][] ret = new int[d.length][d[0].length][d[0][0].length];
int[] minAccp = new int[p.length];
int[] minAccq = new int[q.length];
int[] maxAccp = new int[p.length];
int[] maxAccq = new int[q.length];
for(int i=0;i<p.length;i++){
for(int j=0;j<q.length;j++){
//System.out.println("in repeat d " + i + " " + j + " " + d[i][j][0] + " " + d[i][j][1]);
minAccp[i]+=d[i][j][0];
minAccq[j]+=d[i][j][0];
maxAccp[i]+=d[i][j][1];
maxAccq[j]+=d[i][j][1];
}
}
for(int i=0;i<p.length;i++){
for(int j=0;j<q.length;j++){
int tmpminp = p[i]-(minAccp[i]-d[i][j][0]);
int tmpminq = q[j]-(minAccq[j]-d[i][j][0]);
int tmpmin = Math.min(tmpminp,tmpminq);
ret[i][j][1] = Math.min(d[i][j][1], tmpmin);
int tmpmaxp = p[i]-(maxAccp[i]-d[i][j][1]);
int tmpmaxq = q[j]-(maxAccq[j]-d[i][j][1]);
int tmpmax = Math.max(tmpmaxp,tmpmaxq);
ret[i][j][0] = Math.max(d[i][j][0], tmpmax);
}
}
return ret;
}
public static boolean identicalInt3(int[][][] a,int[][][] b){
boolean ret = true;
for(int i=0;i<a.length;i++){
for(int j=0;j<a[i].length;j++){
for(int k=0;k<a[i][j].length;k++){
if(a[i][j][k] != b[i][j][k]){
ret = false;
return ret;
}
}
}
}
return ret;
}
public static boolean deltaChecker0(int[][] d){
boolean ret = true;
for(int i=0;i<d.length;i++){
for(int j=0;j<d[i].length;j++){
if(d[i][j]<0){
ret = false;
return ret;
}
}
}
return ret;
}
public static boolean deltaCheckerAllZero(int[][] d){
boolean ret = true;
for(int i=0;i<d.length;i++){
for(int j=0;j<d[i].length;j++){
if(d[i][j]!=0){
ret = false;
return ret;
}
}
}
return ret;
}
