CNV多型をはじめとする多(3以上)アレルのHWE検定と関連検定 2

最大でCなるコピー数アレルを有するCNVを考える。
集団内に実在するか否かは考慮せず、コピー数が0-CC+1種類のアレルが存在するものとする。
アレル数はNA=C+1である。
ディプロタイプを考える。
2アレルを区別して特定するようなジェノタイプ(ジェノタイプx)はNGx=\frac{NA(NA-1)}{2}=\frac{C(C+1)}{2}種類ある。
また、2アレルタイプを区別することができず、ディプロイドが持つコピー数によって、ジェノタイプとする(ジェノタイプy)ときには、最少コピー数が0で、最大コピー数が2Cであるから、NGy=2C+1=2NA-1種類ある。
今、ジェノタイプxのようなデータが得られる場合と、ジェノタイプyのようなデータが得られる場合とがあるとする。

  • ジェノタイプxの場合
    • 観測サンプルにおける、コピー数別染色体本数がジェノタイプxから確定的に算出できる
    • 確定的に算出された染色体本数からモーメントとして求めた、コピー数頻度の最尤推定値を基準として、HWEの場合の期待ジェノタイプx観測数が算出できる。この分布についてカイ自乗統計量を算出することができる。
    • このようにして算出したカイ自乗統計量は、自由度NGx-NAにて評価してHWE帰無仮説の棄却率を求めることができる。
      • ジェノタイプxの観測数を変動させることを考えたとき、ジェノタイプxのタイプ数を表す変数がNGxあり、それを制約する式がNA個できるので、自由度をNGx-NAとして評価する。(ただし、ケースコントロールを通じて観測されないジェノタイプはNGxの勘定に入れない。NAも同様)
    • ケース・コントロール関連検定については、こちらのとおり。
  • ジェノタイプyの場合
    • 観測サンプルにおける、コピー数別染色体本数は、ジェノタイプyから確定的に算出できない。
      • ひとつの方法としては、EMアルゴリズムをまわすことで、コピー数アレルの頻度の最尤推定値を求めることもありである(末尾にjavaソースを掲載)
    • 推定されたコピー数アレル頻度から、ジェノタイプyの観測期待値を算出し、それに対して、カイ自乗検定を行うことが可能。こちらの自由度は、NA-1。
      • ジェノタイプyの観測数を変動させることを考えたとき、ジェノタイプyのタイプ数を表す変数がNGyあり、それを制約する変数がNA個あるので、自由度をNGy-NA=2NA-1-NA=NA-1として評価する。ただし、ケースコントロールを通じて観測されないジェノタイプはNGyの勘定に入れない。NAも同様。

※SNPのように2アレルの場合には、ジェノタイプx型の場合しかないが、仮にジェノタイプy型であるとしても、自由度はいずれも1である。これは、上記の論理で作成した計算ソースにおいて、C=2で実行したときに、ジェノタイプx型での評価とジェノタイプy型での評価が同一になることからも確認できる。

    • ケース・コントロール関連検定については、同一コピー数(ジェノタイプy)をもたらすジェノタイプxを区別していない観測データなので、そのようなテーブルをあるがままに検定すればよい。計算上は、ジェノタイプx型のそれと同じ手続きで算出可能である。

ソースはメモ代わりなので、使えない(嘘)の部分も含み、また、バグ持ちの可能性も大だが、以下のとおり。

package CNV;

import StatUtilsX.MiscUtil;

public class CNVcalculator {
	/*
	 * アレル単位での最大コピー数をK
	 * 最小コピー数をkとしたとき
	 * Na=K-k+1
	 * Ng=Na*(Na+1)/2
	 * とすること!
	 */

	public int Na; //No. allele
	public int Ng; //No. genotype
	public double[][] t;//contingency table
	public int[][] t2;//CopyNumber別ジェノタイプテーブル
	public int[][] t3;//allele本数テーブル
	//public double[] w;//weight vector
	public double N;//No.samples
	public double R;//No.cases
	public double S;//No.controls
	public double[] n;//number vector
	int df=1;//自由度
	public org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf1 = 
		new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df);
	public org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidfX;
	public double[][] hweChi;
	public double[][] hweExact;
	public double[][] assoc;//[0] trend,[1] MWW,[2] CNpool
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO 自動生成されたメソッド・スタブ
		int cnvN=100;
		int caseN=500;
		int contN=500;
		int numallele =4;
		int numgenotype = (int)(numallele*(numallele+1)/2);
		int x=10000000;
		int seed=(int)(Math.random()*x);
		Utils.MersenneTwisterFast mz= new Utils.MersenneTwisterFast(seed);
		double[] sl = StatUtilsX.Fisher.serialLogfact((caseN+contN)*2);
		
		for(int i=0;i<cnvN;i++){
			int[][] table=new int[2][numgenotype];
			double[] weight = new double[numgenotype];
			for(int j=0;j<weight.length;j++){
				weight[j]=j;
			}
			double[] af=new double[numallele];
			for(int j=0;j<af.length;j++){
				af[j]=mz.nextDouble();
				
			}
			StatUtilsX.MiscUtil.standard(af);
			//af=StatUtilsX.MiscUtil.accumstandard(af);
			/*
			for(int j=0;j<af.length;j++){
				System.out.println(af[j]);
			}
			*/
			//System.out.println("##");
			double[] gf = new double[numgenotype];
			//System.out.println("numgenotype="+numgenotype);
			int counter=0;
			for(int j=0;j<af.length;j++){
				for(int k=j;k<af.length;k++){
					double tmp = af[j]*af[k];
					if(j==k){
						gf[counter]=tmp;
					}else{
						gf[counter]=tmp*2;
					}
					
					counter++;
				}
			}
			/*
			for(int j=0;j<gf.length;j++){
				System.out.println(gf[j]+"\t");
			}
			*/
			
			gf=StatUtilsX.MiscUtil.accumstandard(gf);
			/*
			for(int j=0;j<gf.length;j++){
				System.out.println(gf[j]+"\t");
			}
			*/
			for(int j=0;j<caseN;j++){
				double r=mz.nextDouble();
				for(int k=0;k<gf.length;k++){
					if(r<gf[k]){
						table[0][k]++;
						break;
					}
				}
			}
			for(int j=0;j<contN;j++){
				double r=mz.nextDouble();
				for(int k=0;k<gf.length;k++){
					if(r<gf[k]){
						table[1][k]++;
						break;
					}
				}
			}
			double hweexactpCase=StatUtilsX.Fisher.hweFisherAbecasis(table[0]);
			double hweexactpCont=StatUtilsX.Fisher.hweFisherAbecasis(table[1]);
			//double hweexactpCase=StatUtilsX.Fisher.hweFisherAbecasis(table[0]);
			//System.out.print("casehwe\t"+hweexactpCase+"\t");
			//System.out.print("conthwe\t"+hweexactpCont+"\t");
			
			CNVcalculator cnvc=new CNVcalculator(table);
			
			cnvc.assoc[0]= cnvc.testcnvTrend(weight);
			cnvc.assoc[2] = cnvc.testCNpool();
			boolean cor=false;
			cnvc.assoc[1] = cnvc.testMWW(cor);
			cnvc.hweChi = cnvc.hweChicasecontall();
			//cnvc.hweExact=cnvc.hweExactcasecontall(sl);
			//System.out.println();
			/*
			for(int j=0;j<cnvc.t2[0].length;j++){
				System.out.print(cnvc.t2[0][j]+"\t");
			}
			*/
			for(int j=0;j<cnvc.t2.length;j++){
				double[] emaf = emCNV(cnvc.t2[j],cnvc.Na);
				for(int k=0;k<emaf.length;k++){
					System.out.print(emaf[k]+"\t");
				}
				//System.out.println();
				double[] hwepool = hweChiPool(cnvc.t2[j],emaf);
				System.out.print(hwepool[0]+"\t"+hwepool[1]+"\t");
				double[] hwegenFromAf = cnvc.hweChiFromAf(cnvc.t[j],emaf);
				System.out.print(hwegenFromAf[0]+"\t"+hwegenFromAf[1]+"\t");
				
			}
			
			cnvc.Print("\t", "\t");
			
		}
		/*
		int[][] t = {{1,2,3,4,5,6},{2,3,4,5,6,7}};
		//int[][] t = {{5,10,3},{8,9,2}};
		CNVcalculator cnv = new CNVcalculator(t);

		double[] w = {1,2,3,4,5,6};
		cnv.assoc[0]= cnv.testcnvTrend(w);
		cnv.assoc[2] = cnv.testCNpool();
		boolean cor=false;
		cnv.assoc[1] = cnv.testMWW(cor);
		cnv.hweChi = cnv.hweChicasecontall();
		cnv.Print("\t", "\n");
		*/
		/*
		for(int i=0;i<hwe.length;i++){
			for(int j=0;j<hwe[i].length;j++){
				System.out.print(hwe[i][j]+"\t");
			}
			System.out.println();
			
		}
		*/
	}
	
	public CNVcalculator(int[][] t_){
		t = new double[t_.length][t_[0].length];
		for(int i=0;i<t.length;i++){
			for(int j=0;j<t[i].length;j++){
				t[i][j]=(double)t_[i][j];
			}
		}
		//t = StatUtilsX.MiscUtil.DeepCopyInt2(t_);
		Ng=t[0].length;
		
		Na=(int)(Math.round((Math.sqrt(1+8*(double)Ng)-1)/2));
		//System.out.println(""+Na);
		N=0;
		R=0;
		S=0;
		n = new double[t[0].length];
		for(int i=0;i<t.length;i++){
			for(int j=0;j<t[i].length;j++){
				N+=t[i][j];
				if(i==0){
					R+=t[i][j];
				}else if(i==1){
					S+=t[i][j];
				}
				n[j]+=t[i][j];
			}
		}
		int maxCN=(Na-1)*2;
		
		t2= new int[2][maxCN+1];
		t3=new int[2][Na];
		int a1=0;
		int a2=0;
		
		for(int i=0;i<t[0].length;i++){
			int numcp=a1+a2;
			t2[0][numcp]+=t[0][i];
			t2[1][numcp]+=t[1][i];
			t3[0][a1]+=t[0][i];
			t3[0][a2]+=t[0][i];
			t3[1][a1]+=t[1][i];
			t3[1][a2]+=t[1][i];
			a2++;
			if(a2==Na){
				
				a1++;
				a2=a1;
			}
		}
		int dfX=Ng-Na;
		//System.out.print("df="+dfX);
		chidfX = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(dfX);
		System.out.print("dfX\t"+dfX+"\t");
		hweChi = new double[3][2];
		hweExact = new double[3][2];
		assoc = new double[3][2];
		/*
		for(int i=0;i<t2[0].length;i++){
			System.out.println(t2[0][i]+"\t"+t2[1][i]);
		}
		System.out.println("");
		for(int i=0;i<t3[0].length;i++){
			System.out.println(t3[0][i]+"\t"+t3[1][i]);
		}
		*/
		//w = new double[Ng];
	}

	public CNVcalculator(int[][] t_,String st){
		/*
		 * trendテストをするため
		 * t_はアレル数・ジェノタイプ数の基準を満たしていなくてもよい
		 */
		t = new double[t_.length][t_[0].length];
		for(int i=0;i<t.length;i++){
			for(int j=0;j<t[i].length;j++){
				t[i][j]=(double)t_[i][j];
			}
		}
		//t = StatUtilsX.MiscUtil.DeepCopyInt2(t_);
		
		Ng=t[0].length;
		/*
		Na=(int)(Math.round((Math.sqrt(1+8*(double)Ng)-1)/2));
		//System.out.println(""+Na);
		 */
		N=0;
		R=0;
		S=0;
		n = new double[t[0].length];
		for(int i=0;i<t.length;i++){
			for(int j=0;j<t[i].length;j++){
				N+=t[i][j];
				if(i==0){
					R+=t[i][j];
				}else if(i==1){
					S+=t[i][j];
				}
				n[j]+=t[i][j];
			}
		}
		/*
		int maxCN=(Na-1)*2;
		
		t2= new int[2][maxCN+1];
		t3=new int[2][Na];
		int a1=0;
		int a2=0;
		
		for(int i=0;i<t[0].length;i++){
			int numcp=a1+a2;
			t2[0][numcp]+=t[0][i];
			t2[1][numcp]+=t[1][i];
			t3[0][a1]+=t[0][i];
			t3[0][a2]+=t[0][i];
			t3[1][a1]+=t[1][i];
			t3[1][a2]+=t[1][i];
			a2++;
			if(a2==Na){
				
				a1++;
				a2=a1;
			}
		}
		chidfX = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(t3[0].length-1);
		
		hweChi = new double[3][2];
		assoc = new double[3][2];
		*/
		/*
		for(int i=0;i<t2[0].length;i++){
			System.out.println(t2[0][i]+"\t"+t2[1][i]);
		}
		System.out.println("");
		for(int i=0;i<t3[0].length;i++){
			System.out.println(t3[0][i]+"\t"+t3[1][i]);
		}
		*/
		//w = new double[Ng];
	}
	public void Print(String sep1,String sep2){
		System.out.println(StringCNVcalculator(sep1,sep2));
	}
	public String StringCNVcalculator(String sep1,String sep2){
		String ret ="";
		ret += "contTable"+sep2;
		for(int i=0;i<t.length;i++){
			for(int j=0;j<t[i].length;j++){
				ret+=t[i][j]+sep1;
			}
			ret+=sep2;
		}
		for(int i=0;i<n.length;i++){
			ret+=n[i]+sep1;
		}
		ret +=sep2;
		ret+="t2"+sep2;
		for(int i=0;i<t2.length;i++){
			for(int j=0;j<t2[i].length;j++){
				ret+=t2[i][j]+sep1;
			}
			ret+=sep2;
		}
		ret+="t3"+sep2;
		for(int i=0;i<t3.length;i++){
			for(int j=0;j<t3[i].length;j++){
				ret+=t3[i][j]+sep1;
			}
			ret+=sep2;
		}
		ret += "R"+sep1+R+sep1+"S"+sep1+S+sep1+"N"+sep1+N+sep2;
		ret += "Na"+sep1+Na+sep2;
		ret += "Ng"+sep1+Ng+sep2;
		
		ret += "hweChi"+sep2;
		for(int i=0;i<hweChi.length;i++){
			for(int j=0;j<hweChi[i].length;j++){
				ret += hweChi[i][j]+sep1;
			}
			ret +=sep2;
		}
		ret += "hweExact"+sep2;
		for(int i=0;i<hweExact.length;i++){
			for(int j=0;j<hweExact[i].length;j++){
				ret += hweExact[i][j]+sep1;
			}
			ret +=sep2;
		}
		ret+="assoc"+sep2;
		for(int i=0;i<assoc.length;i++){
			for(int j=0;j<assoc[i].length;j++){
				ret+=assoc[i][j]+sep1;
			}
			ret+=sep2;
		}
		return ret;
	}
	public double[] testcnvTrend(double[] w){
		double[] ret = new double[2];
		if(N==0){
			return ret;
		}
		
		
		double[] p = new double[Ng];
		double P = R/N;
		for(int i=0;i<p.length;i++){
			if(n[i]!=0){
				p[i]=t[0][i]/n[i];
			}
			
		}
		double[] wn = new double[Ng];
		double W = 0;
		for(int i=0;i<wn.length;i++){
			wn[i]=n[i]*w[i];
			W+=wn[i];
		}
		double Wav=W/N;
		double[] wcor=new double[Ng];
		for(int i=0;i<wcor.length;i++){
			wcor[i]=w[i]-Wav;
		}
		double[] rsqn=new double[Ng];
		double Rsqn=0;
		for(int i=0;i<rsqn.length;i++){
			rsqn[i]=t[0][i]*t[0][i]/n[i];
			Rsqn+=rsqn[i];
		}
		double[] numer=new double[Ng];
		double Numer=0;
		for(int i=0;i<numer.length;i++){
			numer[i]=n[i]*(p[i]-P)*wcor[i];
			Numer+=numer[i];
		}
		double[] denom=new double[Ng];
		double Denom=0;
		for(int i=0;i<denom.length;i++){
			denom[i]=n[i]*wcor[i]*wcor[i];
			Denom+=denom[i];
		}
		ret[0]=Numer*Numer/(Denom*P*(1-P));
		/*
		ret[0]=(N*Math.pow((N*(t[0][1]+2*t[0][2])-R*(n[1]+2*n[2])),2))/
		(R*S*(N*(n[1]+4*n[2])-Math.pow((n[1]+2*n[2]),2)));
	*/
		try{
			ret[1]=1-chidf1.cumulativeProbability(ret[0]);
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		//double s=Math.pow(3,2);
		//System.out.println("cnvtrend\t"+ret[0]+"\t"+ret[1]);
		return ret;
	}
	public double[] testMWW(boolean correct){
		/*
		 * Mann-Whiteney-Willcoxon
		 */
		double[] ret = new double[2];
		/*
		int maxCN=(Na-1)*2;
		
		int[][] t2= new int[2][maxCN+1];
		int a1=Na-1;
		int a2=Na-1;
		
		for(int i=0;i<t[0].length;i++){
			int numcp=a1+a2;
			t2[0][numcp]+=t[0][i];
			t2[1][numcp]+=t[1][i];
			a2--;
			if(a2<0){
				
				a1--;
				a2=a1;
			}
		}
		*/
		int count=0;
		int counted=1;
		double countedNum=S;
		if(R>S){
			count=1;
			counted=0;
			countedNum=R;
		}
		
		double U=0;
		//System.out.println("countedNum="+countedNum);
		for(int i=0;i<t2[count].length;i++){
			countedNum-=t2[counted][i];
			//System.out.println("i="+i+"\t"+countedNum+"\t"+t2[count][i]);
			U+=t2[count][i]*(countedNum+(double)t2[counted][i]/2);
		}
		double tmpU=R*S-U;
		U=Math.min(U,tmpU);
		double[] tie=new double[t2[0].length];
		double Tsum=0;
		for(int i=0;i<t2[0].length;i++){
			tie[i]=t2[0][i]+t2[1][i];
			Tsum+=tie[i]*tie[i]*tie[i]-tie[i];
		}
		double V=R*S*(N*N*N-N-Tsum)/12/(N*N-N);
		double E=R*S/2;
		double cor = 0;
		if(correct){
			cor=0.5;
		}else{
			
		}
		double Z=(Math.abs(U-E)-cor)/Math.sqrt(V);
		ret[0]=Z;
		System.out.print("U=\t"+U+"\t");
		System.out.print("V=\t"+V+"\t");
		System.out.print("E=\t"+E+"\t");
		org.apache.commons.math.distribution.NormalDistributionImpl norm = 
			new org.apache.commons.math.distribution.NormalDistributionImpl();
		try{
			ret[1]= (1-norm.cumulativeProbability(Z))*2;
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		//System.out.println("Z\t"+Z+"\t"+"p\t"+ret[1]);
		return ret;
	}
	public double[] testCNpool(){
		/*
		int maxCN=(Na-1)*2;
		
		int[][] t2= new int[2][maxCN+1];
		int a1=Na-1;
		int a2=Na-1;
		for(int i=0;i<t[0].length;i++){
			int numcp=a1+a2;
			t2[0][numcp]+=t[0][i];
			t2[1][numcp]+=t[1][i];
			a2--;
			if(a2<0){
				
				a1--;
				a2=a1;
			}
		}
		*/
		String st="";
		CNVcalculator cnv=new CNVcalculator(t2,st);
		double[] w= new double[t2[0].length];
		w[0]=2*(Na-1);
		for(int i=1;i<w.length;i++){
			w[i]=w[i-1]-1;
		}
		double[] ret = cnv.testcnvTrend(w);
		
		
		
		return ret;
	}
	public double[][] hweChicasecontall(){
		double[][] ret = new double[3][2];
		
		if(N==0){
			return ret;
		}
		if(R!=0){
			ret[0]=hweChi(t[0],t3[0]);	
		}
		if(S!=0){
			ret[1]=hweChi(t[1],t3[1]);
		}
		int[] n3 = new int[t3[0].length];
		for(int i=0;i<n3.length;i++){
			n3[i]=t3[0][i]+t3[1][i];
		}
		ret[2]=hweChi(n,n3);
		
		return ret;
	}

	public double[][] hweExactcasecontall(double[] sl){
		double[][] ret = new double[3][2];
		
		if(N==0){
			return ret;
		}
		if(R!=0){
			ret[0]=hweExact(t[0],t3[0],sl);	
		}
		if(S!=0){
			ret[1]=hweExact(t[1],t3[1],sl);
		}
		int[] n3 = new int[t3[0].length];
		for(int i=0;i<n3.length;i++){
			n3[i]=t3[0][i]+t3[1][i];
		}
		ret[2]=hweExact(n,n3,sl);
		
		return ret;
	}
	public static double[] hweChiPool(int[] t,double[] af){
		double[] ret = new double[2];
		double[] exp=new double[t.length];
		double sum =0;
		
		for(int i=0;i<t.length;i++){
			sum+=t[i];
		}
		for(int i=0;i<af.length;i++){
			for(int j=0;j<af.length;j++){
				int cp=i+j;
				if(cp<t.length){
					exp[cp]+=af[i]*af[j];
				}
				
			}
		}
		for(int i=0;i<exp.length;i++){
			exp[i]*=sum;
		}
		
		for(int i=0;i<exp.length;i++){
			if(exp[i]>0){
				ret[0]+=(t[i]-exp[i])*(t[i]-exp[i])/exp[i];
			}
			
		}
		//int df=(af.length)*(af.length+1)/2-(af.length);//自由度
		int df=(af.length)-1;//自由度
		System.out.print("df=\t"+df+"\t");
		org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidfY = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df);
		
		try{
			ret[1]=1-chidfY.cumulativeProbability(ret[0]);
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		return ret;
	}
	public double[] hweChiFromAf(double[] t,double[] af){
		
		double[] ret = new double[2];
		double[] exp=new double[t.length];
		double sum =0;
		for(int i=0;i<t.length;i++){
			sum+=t[i];
		}
		
		int counter=0;
		for(int i=0;i<af.length;i++){
			for(int j=i;j<af.length;j++){
				double tmp = af[i]*af[j]*sum;
				if(i==j){
					exp[counter]=tmp;
				}else{
					exp[counter]=tmp*2;
				}
				if(exp[counter]!=0){
					//System.out.println(t[counter]+"\t"+exp[counter]+"\t");
					double ttt =(t[counter]-exp[counter])*(t[counter]-exp[counter])/exp[counter];
					//System.out.println("ttt\t"+ttt); 
					ret[0]+=(t[counter]-exp[counter])*(t[counter]-exp[counter])/exp[counter];
					
				}
				counter++;
			}
		}
		/*
		for(int i=0;i<exp.length;i++){
			System.out.println(t[i]+"\t"+exp[i]);
		}
		*/
		try{
			ret[1]=1-chidfX.cumulativeProbability(ret[0]);
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		return ret;
	}

	public double[] hweChi(double[] t,int[] t3){
		
		double[] ret = new double[2];
		double[] exp=new double[t.length];
		double sum2 =0;
		double sum = 0;
		for(int i=0;i<t3.length;i++){
			sum2+=t3[i];
			
		}
		for(int i=0;i<t.length;i++){
			sum+=t[i];
		}
		int counter=0;
		for(int i=0;i<t3.length;i++){
			for(int j=i;j<t3.length;j++){
				double tmp = ((double)(t3[i])/sum2)*((double)(t3[j])/sum2)*sum;
				if(i==j){
					exp[counter]=tmp;
				}else{
					exp[counter]=tmp*2;
				}
				if(exp[counter]!=0){
					//System.out.println(t[counter]+"\t"+exp[counter]+"\t");
					ret[0]+=(t[counter]-exp[counter])*(t[counter]-exp[counter])/exp[counter];
					
				}
				counter++;
			}
		}
		/*
		for(int i=0;i<exp.length;i++){
			System.out.println(t[i]+"\t"+exp[i]);
		}
		*/
		try{
			ret[1]=1-chidfX.cumulativeProbability(ret[0]);
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		return ret;
	}

	public double[] hweExact(double[] t,int[] t3,double[] sl){
		int numhetero=0;
		int numchrom=0;
		int numsubj=0;
		int[] limit = new int[t.length];
		int a1=0;
		int a2=0;
		for(int i=0;i<limit.length;i++){
			numchrom+=2*t[i];
			numsubj+=t[i];
			if(a1==a2){
				limit[i]=t3[a1]/2;
			}else{
				limit[i]=Math.min(t3[a1], t3[a2]);
				numhetero+=t[i];
			}
			a2++;
			if(a2==t3.length){
				a1++;
				a2=a1;
			}
		}
		double prob=0;
		double commune=-sl[numchrom]+sl[numsubj];
		for(int i=0;i<t3.length;i++){
			commune+=sl[t3[i]];
		}
		double originalpart=numhetero*sl[2];
		for(int i=0;i<t.length;i++){
			originalpart-=sl[(int)(t[i])];
		}
		//System.out.print(originalpart+"\t");
		/*
		for(int i=0;i<limit.length;i++){
			System.out.println("limit\t"+limit[i]);
		}
		*/
		int[] data=new int[t.length];
		boolean yes=true;
		
		while(yes){
			int tmpnumhetero=0;
			double tablepart=0;
			boolean ok=true;
			int[] allelesum=new int[t3.length];
			int b1=0;
			int b2=0;
			for(int i=0;i<data.length;i++){
				allelesum[b1]+=data[i];
				allelesum[b2]+=data[i];
				if(allelesum[b1]>t3[b1]){
					ok=false;
					break;
				}
				if(allelesum[b2]>t3[b2]){
					ok=false;
					break;
				}
				if(b1!=b2){
					tmpnumhetero+=data[i];
				}
				b2++;
				if(b2>=t3.length){
					b1++;
					b2=b1;
				}
			}
			if(ok){
				for(int i=0;i<t3.length;i++){
					if(allelesum[i]!=t3[i]){
						ok=false;
					}
				}
			}
			//ok=true;
			if(ok){
				/*
				for(int i=0;i<data.length;i++){
					System.out.print(data[i]+"\t");
				}
				*/
				tablepart = tmpnumhetero*sl[2];
				for(int i=0;i<data.length;i++){
					tablepart-=sl[(int)(data[i])];
				}
				//System.out.print(tablepart+"\t");
				double tmpPr=Math.exp(commune+tablepart);
				//System.out.print(tmpPr+"\t");
				//System.out.println();
				if(tablepart<=originalpart){
					prob+=Math.exp(commune+tablepart);
				}
			}
			data[0]++;
			for(int i=0;i<limit.length-1;i++){
				if(data[i]>limit[i]){
					data[i]=0;
					data[i+1]++;
				}
			}
			if(data[data.length-1]==limit[data.length-1]+1){
				yes=false;
			}
		}
		double[] ret = new double[2];
		ret[1]=prob;
		ret[0]=Math.exp(commune+originalpart);
		//System.out.println("Prob\t="+ret[0]+"\t"+ret[1]);
		return ret;
	}
	public static double[] emCNV(int[] d,int max){
		/*
		 * dは、コピー数0からコピー数d.lengthまでのベクトル
		 */
		if(max==0){
			max=d.length;
		}
		/*
		else if(max<d.length){
			max=d.length;
		}
		*/
		double[] tmpret = new double[max];
		double[] ret = new double[max];
		double sum = 0;
		for(int i=0;i<d.length;i++){
			sum+=d[i];
		}
		sum*=2;
		//System.out.println("d.length="+d.length+" max="+max);
		for(int i=0;i<d.length;i++){
			double tmp = 1/(double)(i+1)/sum;
			for(int j=0;j<=i;j++){
				//System.out.println("i="+i+" j="+j);
				if(i-j<tmpret.length && j<tmpret.length){
					tmpret[j]+=tmp*d[i];
					tmpret[i-j]+=tmp*d[i];
				}
				
				
				
			}
			
		}
		/*
		for(int i=0;i<tmpret.length;i++){
			tmpret[i]=1/(double)tmpret.length;
		}
		*/
		
		ret = MiscUtil.DeepCopyDouble1(tmpret);
		//for(int i=0;i<ret.length;i++){
			//System.out.print(ret[i]+"\t");
		//}
		//System.out.println();
		tmpret=new double[max];
		int numiter=1000;
		for(int x=0;x<numiter;x++){
			for(int i=0;i<d.length;i++){
				double tmp = 1/(double)(i+1)/sum;
				double probsum=0;
				
				for(int j=0;j<=i;j++){
					if(i-j<tmpret.length && j<tmpret.length){
						//System.out.println("i="+i+" j="+j+" tmpretlen "+tmpret.length);
						double prob = ret[j]*ret[i-j];
						probsum+=prob;
					}
					
					
				}
				if(probsum==0){
					
				}else{
					for(int j=0;j<=i;j++){
						if(i-j<tmpret.length && j<tmpret.length){
							double prob = ret[j]*ret[i-j]/probsum/sum;
							//probsum+=prob;
							
							tmpret[j]+=prob*d[i];
							tmpret[i-j]+=prob*d[i];
						}
						
					}
				}
				
			}
			ret = MiscUtil.DeepCopyDouble1(tmpret);
			tmpret=new double[max];
			//for(int i=0;i<ret.length;i++){
				//System.out.print(ret[i]+"\t");
			//}
			//System.out.println();
		}
		return ret;
	}

}