exactP値の確率
ある分割表データがある。
その周辺度数がある。
その周辺度数からは、有限個パターンの分割表が作れる。
それらが与えるexactP値は有限個である(分割表の個数以下。異なる分割表から同一のexactP値が与えられることもあるので必ずしも一致しない）。
そのexactP値にも出やすい値と出にくい値がある。
exactP値は小さいものよりも大きいものの方が出やすい。
ただし、小さい値は互いに似通った値(小さいもの同士)がたくさんあり、それぞれの値の出る確率が低いのに対して、大きなexactP値は、そのものずばりの値が出る確率は高いが、それらと近いexactP値をとれないことから、exactP値について、等間隔のヒストグラムを描けば、小さいexactP値と大きいexactP値とで、均衡する。
掲載図は、掲載ソースによって出力したもののプロット。分割表は、２ｘ２で、周辺度数は、行１、行２が1800、列１が2200、列２が1400の場合である。X軸がexactPの値、Y軸がその生起確率である。原点に近い方でプロットが密であり、exactPが大きくなるに連れて、疎になっている。exactPが大きい方が生起確率が高くなっている。
関連記事はこちら。
public static double[][] ProbMoreThanGivenReturnDataAll(int[][] da,double[] logFactList){
double V=2;
double[][] ret = ProbMoreThanGivenReturnData(da,V);
return ret;
}
public static double[][] ProbMoreThanGivenReturnData(int[][] da,double V){
/*
* nxm分割表のFisher正確確率検定
*/
/*
* 関数limitCubeは分割表から、
* 各セルの取りうる整数値の範囲を限定する
* 各セルが取りうる値の範囲は、周辺度数の制約と
* 各セルが0以上である、という制約とから定まる
* 各セルの取りうる範囲は線形計画法的な多次元立体（斜面を持つ）
* で囲まれるが
* ここでは、多次元立方体(斜面なし）の範囲を定め、
* その範囲の点をしらみつぶしにする
*/
double[][] ret =new double[0][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;
double[] prob = {};
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]];
}
}
boolean tobeadded = false;
/*
* 引数で取った値との大小比較
*/
double tmplogfact=constant-tmpdenom;
double logV = Math.log(V);
//double tmpexp = Math.exp(tmplogfact);
if(MiscUtil.IdenticalDouble(tmplogfact,logV,MiscUtil.DEV)){
tobeadded = true;
}else if(tmplogfact<logV){
tobeadded = true;
}
/*
if(MiscUtil.IdenticalDouble(tmpdenom,value0,MiscUtil.DEV)){
tobeadded = true;
}else if(tmpdenom>value0){
tobeadded = true;
}
*/
//if(tmpdenom>=value0){
if(tobeadded){
prob = MiscUtil.AddElemDouble1(prob,tmplogfact);
//System.out.println("prol len " + prob.length);
/*
String sss="";
for(int i=0;i<tmpvar.length;i++){
for(int j=0;j<tmpvar[i].length;j++){
sss += tmpvar[i][j] + "\t";
}
}
*/
//System.out.println(sss);
//double tmplogfact=constant-tmpdenom;
//double tmpexp = Math.exp(tmplogfact);
//System.out.println("tmpexp\t" + tmpexp);
//exactP += tmpexp;
//}
}
}
/*
* nshinをひとつ大きくする
*/
nshin = nshinUpper(nshin,limitnum);
/*
* nshinがぐるっと回って、初めに戻ったら
* 網羅ループ終了
*/
//累積で足して行って、基準以下なら足す
boolean zerocheck = deltaCheckerAllZero(nshin);
if(zerocheck){
returned =false;
}
}
MiscUtil.bublSortDouble(prob);
exactP = 0;
double exactPstock = 0;
for(int i=0;i<prob.length;i++){
double tmpexp = Math.exp(prob[i]);
//System.out.println("%%tmpexp\t" + tmpexp);
exactPstock += tmpexp;
boolean identical = true;
if(i>0 && i<prob.length-1){
if(MiscUtil.IdenticalDouble(prob[i-1],prob[i],MiscUtil.DEV)){
identical = false;
//System.out.println("identical!");
}
}
if(exactPstock>V){
break;
}else{
if(identical){
//System.out.println("to add");
double dif = exactPstock-exactP;
exactP=exactPstock;
double[] tmp = {exactP,dif};
ret = MiscUtil.AddElemDouble2(ret,tmp);
//System.out.println("exactP\tdif\t"+exactP+"\t"+dif);
}
}
//System.out.println("exactPstock exactP\t"+exactPstock + "\t" + exactP);
}
ret[ret.length-1][0]=1;
return ret;
//return exactP;
}
