ディプロタイプを観測するということ - ryamadaの遺伝学・遺伝統計学メモ

今、全部で $N$ 人いるとする。染色体ハプロイドは $2N$ セットある。１つの２アレル型SNPについて着目すると、２種類のアレル数 $N_A,N_B$ は $N_A+N_B=2N$ を満足する。他方、２アレルが作る３ジェノタイプの人数 $N_{AA},N_{AB},N_{BB}$ は $N_{AA}+N_{AB}+N_{BB}=N$ を満足する。
この分配方式を一般的に定義しなおすと次のようになる。
N個の観測単位(個人)がある。
観測単位は、同数 $k$ の最小単位の実体(アレル)からなる(SNPの場合は、 $k=2$ )。実体総数は $kN$ 個ある。
観測単位は、 $T$ 種類にわけれれ(SNPの場合は $T=2$ )、 $\sum_{t=1}^TN_t=kN$ を満足する。
観測個体が最小単位実体を持つ、そのパターン(SNPの場合はジェノタイプ)は、 $X=\frac{(k+T-1)!}{k!(T-1)!}$ ある(SNPの場合は、 $k=2,T=2$ で、ジェノタイプは $\frac{3!}{2!1!}=3$ タイプ)。パターンごとの人数を $N_{p,x}$ , $x=1,2,...,X$ とすると、 $\sum_{x=1}^XN_{p,x}=N$ を満足する
最小単位実体を中心にNxTテーブル(最小単位実体テーブル)を作ると以下のようになる。

アレル１	アレル２	...	アレルT	合算
k	0	...	0	k
k	0	...	0	k
k-1	1	...	0	k
k-1	0	...	1	k
k-2	2	...	0	k
k-2	1	...	1	k
1	1	...	k-2	k
0	2	...	0	k
...	...	...	...	...
N_1	N_2	...	N_T	kN

ここで、各個人に対応する行を構成する $T$ 個の数値の組は、観測個体が持つ最小単位実体の持ち方のパターンに相当するから、この表において、 $N_{p,x}$ ずつ、同一の行がある。このパターンごとの $T$ 個の値を、 $n_{x,t}$ , $x=1,2,...,X$ , $t=1,2,...T$ とする。それらについて、観測個体数に着目してテーブル(観測個体テーブル)を作成すると、1xXの表になる。

パターン１	パターン２	...	パターンX	合算
N_{p,1}	N_{p,2}	...	N_{p,X}	N

最小単位実体テーブルについて、その観測確率をその周辺度数から求めるには、
$P=\frac{k!^N\prod_{t=1}^TN_t!}{(kN)!\prod_{x=1}^X(\prod_{t=1}^Tn_{x,t}!)^{N_{p,x}})}$

ここで、観測個体のパターンが同じ個体について、区別しないこととすると、観測個体テーブルを観測するような確率は
$Q=P\times \frac{N!}{\prod_{i=x}^XN_{p,x}!}$
と表せる。
この式は、HWEの正確検定で用いる、あるアレルの観測本数を前提としたときの、ジェノタイプ数の観測確率に相当している。

もし、観測個体について、[texC]種類のカテゴリにわけるとする(ケース・コントロールの場合は $C=2$ )と、その分割表が

パターン１	パターン２	...	パターンX	合算
m_{1,1}	m_{1,2}	...	m_{1,X}	m_{1,*}
m_{2,1}	m_{2,2}	...	m_{2,X}	m_{2,*}
...	...	...	...	...
m_{C,1}	m_{C,2}	...	m_{C,X}	m_{C,*}
N_{p,1}	N_{p,2}	...	N_{p,X}	N

となって、このようなテーブルを観測する確率は
$R=Q\times \frac{\prod_{c=1}^CM_{c,*}!\prod_{x=1}^XN_{p,x}!}{N!\prod_{c=1}^C(\prod_{x=1}^Xm_{c,x}!)}$
これは、
$R=P\times \frac{\prod_{c=1}^CM_{c,*}!}{\prod_{c=1}^C(\prod_{x=1}^Xm_{c,x}!)}$
とも書ける。
これについての検算用の雑なjavaソースはこちら

public class generalExact {
	/*
	 * kNの実体がT種類に分けられており、
	 * それがk個ずつに分配され、
	 * さらに、それが、
	 * Cカテゴリに分かれているような分配形式である
	 * SNPでケースコントロールの場合は
	 * T=2 C=2である
	 * kはディプロイドなら2
	 * dはTｘNの表である
	 */
	public int k;//1人あたりの染色体本数
	public int N;//人数
	public int T;//アレル種類数
	public int[] numT;//アレル別本数
	public int C;
	public int[] numC;//カテゴリ別人数
	
	public int GenType;//genotype数
	public int[] numGenType;//genotype別人数
	public int[][] table;//2x3
	public int kshinmax;
	public int[][] tableT;//2x2
	//public double[] serialLog;
	public generalExact(int k_,int[][] d,int[] phen,int C_,double[] serialLog){
		k=k_;
		N=d.length;
		//serialLog = serialLog_;
		T=d[0].length;
		numT = new int[T];
		
		GenType=Fisher.combinseriallog(k+T-1,T-1,serialLog);
		C=C_;
		numC = new int[C];
		kshinmax=(int)Math.pow((double)(k+1),(double)(T));
		//System.out.print("kshinmax="+kshinmax+"\n");
		table = new int[C][kshinmax];
		tableT = new int[C][T];
		numGenType = new int[kshinmax];
		for(int i=0;i<d.length;i++){
			int kshinnum=0;
			for(int j=0;j<d[i].length;j++){
				numT[j]+=d[i][j];
				tableT[phen[i]][j]+=d[i][j];
				int tmp=(int)Math.pow((double)(k+1),(double)(j));
				//System.out.println("keta="+tmp);
				//System.out.println("d[i][j]="+d[i][j]);
				kshinnum+=d[i][j]*tmp;				
			}
			//System.out.println("kshinnum\t"+kshinnum);
			table[phen[i]][kshinnum]++;
			numGenType[kshinnum]++;
			numC[phen[i]]++;
		}
	}
	public generalExact(int k_,int[][] d,int[] phen,int C_){
		k=k_;
		N=d.length;
		double[] serialLog = Fisher.serialLogfact(N*k);
		generalExact gE = new generalExact(k,d,phen,C_,serialLog);
		T=gE.T;
		numT = new int[T];
		
		//GenType=Fisher.combinseriallog(k+T-1,T-1,serialLog);
		C=C_;
		numC = gE.numC;
		//new int[C];
		kshinmax=gE.kshinmax;
		//kshinmax=(int)Math.pow((double)(k+1),(double)(T));
		//System.out.print("kshinmax="+kshinmax+"\n");
		//table = new int[C][kshinmax];
		table = gE.table;
		tableT = gE.tableT;
		//tableT = new int[C][T];
		numGenType = gE.numGenType;
		//numGenType = new int[kshinmax];
		
		/*
		for(int i=0;i<d.length;i++){
			int kshinnum=0;
			for(int j=0;j<d[i].length;j++){
				numT[j]+=d[i][j];
				tableT[phen[i]][j]+=d[i][j];
				int tmp=(int)Math.pow((double)(k+1),(double)(j));
				//System.out.println("keta="+tmp);
				//System.out.println("d[i][j]="+d[i][j]);
				kshinnum+=d[i][j]*tmp;				
			}
			//System.out.println("kshinnum\t"+kshinnum);
			table[phen[i]][kshinnum]++;
			numGenType[kshinnum]++;
			numC[phen[i]]++;
		}
		*/
		
	}
	public double[] probForK2(double[] s){
		double[] ret = new double[3];
		double judgeTableT=0;
		for(int i=0;i<tableT.length;i++){
			for(int j=0;j<tableT[i].length;j++){
				judgeTableT+=s[tableT[i][j]];
			}
		}
		System.out.println("judgeTableT=\t"+judgeTableT);
		/*
		 * k=2,T=2,C=2である
		 */
		if(k!=2 || T!=2 || C!=2){
			return ret;
		}
		
		/*
		 * ExactHWEはとりうるすべての３ジェノタイプ数を網羅してくれる
		 * numGenType[2],[4],[6]に値が格納されている
		 */
		boolean evenodd = true;
		if(numGenType[4]%2==1)evenodd=false;
		int minhetero=0;
		int maxhetero=Math.min(numT[0], numT[1]);
		int numPattern=maxhetero-minhetero+1;
		double[] Qs = new double[numPattern];
		double sumQ=0;
		if(evenodd){
			minhetero=0;
			
		}else if(!evenodd){
			minhetero=1;
			
		}
		int here = N+1;
		for(int i=minhetero;i<=maxhetero;i=i+2){
			//System.out.println("i="+i);
			int g1 = (numT[0]-i)/2;
			int g2 = i;
			int g3 = N-(numT[0]+i)/2;
			
			//System.out.println("3types="+g1 +" "+g2+" "+g3);
			//int hetero = allele-i*2;
			//if(hetero<0){
				//continue;
			//}
			//int neghomo = sum-i-hetero;
			//if(neghomo<0){
				//continue;
			//}
			//System.out.println("i hetero neghomo " + i + " " + hetero + " " + neghomo);
			double combin = (s[N]+s[numT[0]]+s[numT[1]])
					-(s[g2] + s[g1] + s[g3] +s[N*2]);
			//System.out.println("combin "+ combin);
			Qs[i]=Math.exp(combin +  i * Math.log(2) );
			sumQ+=Qs[i];
			System.out.print("i=\t"+i +"\tthreeGenotypes\t"+g1+"\t"+g2+"\t"+g3+"\t\tHWE_P_component\t"+Qs[i]+"\n");
			//double P1 = Math.exp(combin +  i * Math.log(2) );
			//System.out.println("P0="+P0+" P1="+P1);
			/*
			 * ３ジェノタイプのとり方ごとに計算する
			 * この３ジェノタイプの分布において、
			 * アレルのケースコントロール分布が同一な２ｘ３表を作れれば作る
			 * 作れなければ
			 */
			int[] margL = numC;
			int[] margC = {g1,g2,g3};
			int N2=N*2;
			
			double[] tmp = extPs(margL,margC,N2,tableT,judgeTableT,s);
			ret[0]+=tmp[0]*Qs[i];
			ret[1]+=tmp[1]*Qs[i];
			System.out.println("extPs=\t"+tmp[0]+"\t"+tmp[1]);
		}	
		System.out.println("Qsum=\t"+sumQ);
		
		
		System.out.print("P=\t"+ret[0]+"\t"+ret[1]+"\n");
		return ret;
	}
	
	public double[] extPs(int[] margL,int[] margC,int sum,int[][] twoBYtwo,double org,double[] logFact){
		double[] ret = {0,0};
		int[][] da = new int[2][3];
		//{{margL[0],0,0},{margL[1],0,0}};
		da[0][0]=Math.min(margL[0], margC[0]);
		da[1][0]=margC[0]-da[0][0];
		int tmpMargL = margL[0]-da[0][0];
		da[0][1]=Math.min(tmpMargL,margC[1]);
		da[1][1]=margC[1]-da[0][1];
		da[0][2]=tmpMargL-da[0][1];
		da[1][2]=margC[2]-da[0][2];
		
		//System.out.println("dummyGenos\n"+da[0][0]+"\t"+da[0][1]+"\t"+da[0][2]+"\n"+da[1][0]+"\t"+da[1][1]+"\t"+da[1][2]);
		//int[][] twoBYtwo = {{da[0][0]*2+da[0][1],da[0][1]+2*da[0][2]},{da[1][0]*2+da[1][1],da[1][1]+2*da[1][2]}};
		//System.out.println("twoBYtwo\n"+twoBYtwo[0][0]+"\t"+twoBYtwo[0][1]+"\t"+twoBYtwo[1][0]+"\t"+twoBYtwo[1][1]);
		int[] Lmarg = {twoBYtwo[0][0]+twoBYtwo[0][1],twoBYtwo[1][0]+twoBYtwo[1][1]};
		int[] Cmarg = {twoBYtwo[0][0]+twoBYtwo[1][0],twoBYtwo[0][1]+twoBYtwo[1][1]};
		//int sum = Lmarg[0]+Lmarg[1];
		/*
		double[][] exp = {{(double)(Lmarg[0]*Cmarg[0])/sum,(double)(Lmarg[0]*Cmarg[1])/sum},
				{(double)(Lmarg[1]*Cmarg[0])/sum,(double)(Lmarg[1]*Cmarg[1])/sum}};
				*/
		int min=Math.max(0, Cmarg[0]-Lmarg[1]);
		min = Math.max(min, Lmarg[0]-Cmarg[1]);
		int max=Math.min(Lmarg[0], Cmarg[0]);
		//System.out.println("min="+min);
		//System.out.println("max="+max);
		//double[] judge = new double[max-min+1];
		int ceilMinSide=-1;
		int floorMaxSide=0;
		int extORint=0;
		int checkCounter=0;
		//double org = logFact[twoBYtwo[0][0]]+logFact[twoBYtwo[0][1]]+logFact[twoBYtwo[1][0]]+logFact[twoBYtwo[1][1]];
		int[] orgfour = {twoBYtwo[0][0],twoBYtwo[0][1],twoBYtwo[1][0],twoBYtwo[1][1]};
		//int[] orgfour = {tableT[0][0],tableT[0][1],table[1][0],tableT[1][1]};
		StatUtils.MiscUtil.bublSort(orgfour);
		//System.out.println("org="+org);
		//if(twoBYtwo[0][0]>=exp[0][0]){
			ceilMinSide=min;
			for(int i=min;i<=max;i++){
				//System.out.println("i="+i);
				checkCounter++;
				int[] tmp = {i,Lmarg[0]-i,Cmarg[0]-i,Lmarg[1]-Cmarg[0]+i};
				int[] sorttmp = StatUtils.MiscUtil.DeepCopyInt1(tmp);
				StatUtils.MiscUtil.bublSort(sorttmp);
				boolean ident = StatUtils.MiscUtil.IdenticalIntList(orgfour, sorttmp);
				//System.out.println("ident="+ident);
				if(ident){
					//System.out.println("i="+i);
					ceilMinSide=i;
					//System.out.println("$$$$$$"+tmp[0]+" "+tmp[1]+" "+tmp[2]+ " "+tmp[3]);
					double tmpval = Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
					ret[0] += tmpval;
					ret[1] += tmpval;
					//System.out.println("tmpval=\t"+tmpval);
				}else{
					double tmplog = logFact[tmp[0]]+logFact[tmp[1]]+logFact[tmp[2]]+logFact[tmp[3]];
					//System.out.println("tmporg="+tmplog);
					if(tmplog>=org){
						//System.out.println("i="+i);
						ceilMinSide=i;
						double tmpval = Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
						ret[0] += tmpval;
						ret[1] += tmpval;
						//System.out.println("tmpval=\t"+tmpval);
					}else{
						ceilMinSide=i;
						/*
						for(int j=min;j<i;j++){
							System.out.println("j="+j);
						}
						*/
						double tmpval = Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
						ret[1] += tmpval;
						//ret += tmpval;
						//System.out.println("tmpval=\t"+tmpval);
						//break;
					}
					
				}
				
				
			}
			floorMaxSide = max;
			//System.out.println("max="+max+"ceilMinSide="+ceilMinSide);
			for(int i=max;i>ceilMinSide;i--){
				checkCounter++;
				int[] tmp = {i,Lmarg[0]-i,Cmarg[0]-i,Lmarg[1]-Cmarg[0]+i};
				
				int[] sorttmp = StatUtils.MiscUtil.DeepCopyInt1(tmp);
				StatUtils.MiscUtil.bublSort(sorttmp);
				boolean ident = StatUtils.MiscUtil.IdenticalIntList(orgfour, sorttmp);
				//System.out.println("ident="+ident);
				if(ident){
					//System.out.println("i="+i);
					floorMaxSide=i;
					//System.out.println("$$$$$$"+tmp[0]+" "+tmp[1]+" "+tmp[2]+ " "+tmp[3]);
					//ret += Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
					
					double tmpval = Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
					ret[0] += tmpval;
					ret[1] += tmpval;
					//System.out.println("tmpval=\t"+tmpval);
					
				}else{
					double tmplog = logFact[tmp[0]]+logFact[tmp[1]]+logFact[tmp[2]]+logFact[tmp[3]];
					//System.out.println("tmporg="+tmplog);
					if(tmplog>org){
						//System.out.println("i="+i);
						//ret += Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
						
						double tmpval = Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
						ret[0] += tmpval;
						ret[1] += tmpval;
						//System.out.println("tmpval=\t"+tmpval);
					}else{
						floorMaxSide=i;
						/*
						for(int j=max;j>i;j--){
							System.out.println("j="+j);
						}
						*/
						double tmpval = Fisher.Prob2x3sfor2x2X(da,tmp,logFact);
						ret[1] += tmpval;
						//ret += tmpval;
						//System.out.println("tmpval=\t"+tmpval);
						//break;
					}
				}
				
				
			}

		
		return ret;
	}
	public static int sumline(int[] l){
		int ret =0;
		for(int i=0;i<l.length;i++){
			ret += l[i];
		}
		return ret;
	}
	public String stringGeneralExact(String sep1,String sep2){
		String ret = "";
		ret += "No.samples"+sep1 + N+sep2;
		ret += "No.objectsPerSample"+sep1+k+sep2;
		ret += "No.alleles"+sep1+T+sep2;
		ret += "No.categories"+sep1+C+sep2;
		ret += "No.genotypes"+sep1+GenType+sep2;
		ret += "No.objectsPerAllele"+sep2;
		for(int i=0;i<numT.length;i++){
			ret += numT[i]+sep1;
		}
		ret += sep2;
		ret += "No.samplesPerCategory"+sep2;
		for(int i=0;i<numC.length;i++){
			ret += numC[i]+sep1;
		}
		ret += sep2;
		ret += "No.samplesPerGenotype"+sep2;
		for(int i=0;i<kshinmax;i++){
			ret += numGenType[i]+sep1;
		}
		ret += sep2;
		ret += "No.samplesPerCategoryPerGenotype"+sep2;
		for(int j=0;j<table.length;j++){
			for(int i=0;i<table[j].length;i++){
				
				ret +=table[j][i]+sep1;
			}
			ret += sep2;
		}
		ret += "No.objectsPerCategoryPerAllele"+sep2;
		for(int j=0;j<tableT.length;j++){
			for(int i=0;i<tableT[j].length;i++){
				
				ret +=tableT[j][i]+sep1;
			}
			ret += sep2;
		}
		
		return ret;
		
	}
	public void printGeneralExact(String sep1,String sep2){
		String ret = stringGeneralExact(sep1,sep2);
		System.out.print(ret);
	}
	/*
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO 自動生成されたメソッド・スタブ
		int k=2;
		//int[] phenotype = {0,0,1,1};
		int[][] fix = {{10,20,10},{70,10,0}};
		int sumfix=0;
		
		for(int i=0;i<fix.length;i++){
			for(int j=0;j<fix[i].length;j++){
				sumfix+=fix[i][j];
			}
		}
		int[][] data2 = new int[sumfix][2];
		int[] phenotype2 = new int[sumfix];
		int counter=0;
		for(int i=0;i<fix.length;i++){
			for(int j=0;j<fix[i].length;j++){
				for(int l=0;l<fix[i][j];l++){
					phenotype2[counter]=i;
					if(j==0){
						data2[counter][0]=2;
						data2[counter][1]=0;
					}else if(j==1){
						data2[counter][0]=1;
						data2[counter][1]=1;
					}else{
						data2[counter][0]=0;
						data2[counter][1]=2;
					}
					counter++;
				}
				
			}
		}
		int N = 8;
		int T = 2;
		int [][] data = new int[N][T];
		int[] phenotype = new int[N];
		//int[][] data = {{2,0,0,0},{1,0,1,0},{0,0,1,1},{0,2,0,0}};
		int C=2;
		int seed = (int)(Math.random()*100000);
		//System.out.println("seed="+seed);
		Utils.MersenneTwisterFast mz = new Utils.MersenneTwisterFast(seed);
		double af = 0.03;
		for(int i=0;i<N;i++){
			int tmpphen = mz.nextInt(C);
			phenotype[i]=tmpphen;
			//System.out.println("tmpphen\t"+tmpphen);
			for(int j=0;j<k;j++){
				int tmp = mz.nextInt(T);
				double tmpd = mz.nextDouble();
				double tmpr = mz.nextDouble();
				if(tmpr<0.95){
//					System.out.println("tmpd="+tmpd);
					if(tmpd<af){
						tmp=0;
					}else{
						tmp=1;
					}
				}else{
					tmp=0;
				}
				
				//System.out.println("tmp\t"+tmp);
				data[i][tmp]++;
			}
		}
		//int[][] data2={{5,10,5},{10,5,5}};
		//double[] serialLog = Fisher.serialLogfact(N*2);
		double[] serialLog = Fisher.serialLogfact(sumfix*2);
		//int[] phenotype2={0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1};
		//generalExact gE = new generalExact(k,data,phenotype,C,serialLog);
		generalExact gE = new generalExact(k,data2,phenotype2,C,serialLog);
		gE.printGeneralExact("\t","\n");
		gE.probForK2(serialLog);
		int[][] twoBYthree={{gE.table[0][2],gE.table[0][4],gE.table[0][6]},
				{gE.table[1][2],gE.table[1][4],gE.table[1][6]}};
		
		int[] testType=new int[18];
		for(int i=0;i<testType.length;i++){
			testType[i]=i;
		}
		double[] singleSNPtest = BBPPhenotypeVector3GC.SingleSNPtest4.PforSelectedTestsOriginal2(twoBYthree,serialLog,testType);
		System.out.print("2x2Fisher\t"+singleSNPtest[2]+"\t2x2Chi\t"+singleSNPtest[3]+"\tTrendExact\t"+singleSNPtest[10]+"\tTrend\t"+singleSNPtest[9]+"\tHWE\t"+singleSNPtest[17]);
		
		//for(int i=0;i<singleSNPtest.length;i++){
			//System.out.print(singleSNPtest[i]+"\t");
		//}
		/*
		System.out.println("k="+gE.k);
		System.out.println("T="+gE.T);
		System.out.println("gentype="+gE.GenType);
		for(int i=0;i<gE.table.length;i++){
			for(int j=0;j<gE.table[i].length;j++){
				System.out.print(gE.table[i][j]+"\t");
			}
			System.out.print("\n");
		}
		*/
	}

}
>||
public static double Prob2x3sfor2x2X(int[][] da, int[] four,double[] logFact){
		double ret=0;
		//System.out.println("&&&&&&&&");
		//System.out.println("four="+four[0]+" "+four[1]+" "+ four[2]+ " "+four[3]);
		int[] Lmarg={da[0][0]+da[0][1]+da[0][2],da[1][0]+da[1][1]+da[1][2]};
		int[] Cmarg={da[0][0]+da[1][0],da[0][1]+da[1][1],da[0][2]+da[1][2]};
		//System.out.println("Lmarg="+Lmarg[0]+" "+Lmarg[1]);
		//System.out.println("Cmarg="+Cmarg[0]+" "+Cmarg[1]+" "+Cmarg[2]);
		
		double constant = -logFact[Lmarg[0]+Lmarg[1]]+logFact[Lmarg[0]]+logFact[Lmarg[1]]
		                          +logFact[Cmarg[0]]+logFact[Cmarg[1]]+logFact[Cmarg[2]];
		
		int[][] tmpda = new int[2][3];
		int max01 = Cmarg[1]-(Lmarg[1]-Cmarg[0]-Cmarg[2]);
		int min02 = Cmarg[2]-Lmarg[1];
		int min = Math.max(0, (int)((four[0]-max01)/2));
		int max = Math.min((int)(four[0]/2), Cmarg[0]);
		//System.out.println("min="+min);
		//System.out.println("max="+max);
		for(int i=min;i<=max;i++){
			tmpda[0][0]=i;
			tmpda[0][1]=(four[0]-2*tmpda[0][0]);
			tmpda[0][2]=Lmarg[0]-tmpda[0][0]-tmpda[0][1];
			if(tmpda[0][2]>=0){
				tmpda[1][1]=Cmarg[1]-tmpda[0][1];
				if(tmpda[1][1]>=0){
					tmpda[1][0]=Cmarg[0]-tmpda[0][0];
					if(tmpda[1][0]>=0){
						tmpda[1][2]=Cmarg[2]-tmpda[0][2];
						if(tmpda[1][2]>=0){
							String st="!!##\t";
							for(int j=0;j<tmpda.length;j++){
								for(int k=0;k<tmpda[j].length;k++){
									st += tmpda[j][k]+"\t";
								}
							}
							//System.out.print(st);
							double variableLog = constant-(logFact[tmpda[0][0]]+logFact[tmpda[0][1]]+logFact[tmpda[0][2]]+
							 						logFact[tmpda[1][0]]+logFact[tmpda[1][1]]+logFact[tmpda[1][2]]);
							//System.out.println("variablelog="+variableLog);
							double tmpvar=Math.exp(variableLog);
							//System.out.print(tmpvar+"\n");
							ret += tmpvar;
							//System.out.println(st + "\t"+Math.exp(variableLog));
						}
						
					}
					
				}
				
			}
			
		}
		return ret;
	}

アレル１	アレル２	...	アレルT	合算
k	0	...	0	k
k	0	...	0	k
k-1	1	...	0	k
k-1	0	...	1	k
k-2	2	...	0	k
k-2	1	...	1	k
1	1	...	k-2	k
0	2	...	0	k
...	...	...	...	...
N_1	N_2	...	N_T	kN

アレル１	アレル２	...	アレルT	合算
k	0	...	0	k
k	0	...	0	k
k-1	1	...	0	k
k-1	0	...	1	k
k-2	2	...	0	k
k-2	1	...	1	k
1	1	...	k-2	k
0	2	...	0	k
...	...	...	...	...
N_1	N_2	...	N_T	kN

アレル１	アレル２	...	アレルT	合算
k	0	...	0	k
k	0	...	0	k
k-1	1	...	0	k
k-1	0	...	1	k
k-2	2	...	0	k
k-2	1	...	1	k
1	1	...	k-2	k
0	2	...	0	k
...	...	...	...	...
N_1	N_2	...	N_T	kN