[[HWE]Wright's F、３ジェノタイプ、アレル頻度 その他

Fについてはこちら。
ジェノタイプ頻度とFの関係については、こちら。
SNPがある。アレル頻度をとする。３ジェノタイプの頻度をとする。



なる関係にあり、

である。
あるアレル頻度において、はのときに最大で１となり、HWEを満足するときに0となり、のときに負の数をとって、最小となり、そのとき、,のいずれかが０となる。
今、同一集団において観測された、相互に独立な複数のSNPの３ジェノタイプ観測値について、共通のFが存在するものとして、その値を求めるとどうなるかについての雑なソースはこちら。

public class EstF {

	public static double THRES=0.0001;
	public static double GETA=0.0000000000001;
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO 自動生成されたメソッド・スタブ
		int N=10000;
		int Nsample=100;
		double fst =0.3;
		int seed = (int)(Math.random()*100000000);
		Utils.MersenneTwisterFast mz=new Utils.MersenneTwisterFast(seed);
		int[][] d2=new int[N][3];
		for(int i=0;i<d2.length;i++){
			double ftmp = fst *(1+mz.nextGaussian());
			if(ftmp<0){
				ftmp=0;
			}else if(ftmp>1){
				ftmp=1;
			}
			double af=Math.random();
			d2[i][0]=(int)(Nsample*(ftmp*af+(1-ftmp)*af*af));
			d2[i][1]=(int)(Nsample*(2*(1-ftmp)*af*(1-af)));
			d2[i][2]=(int)(Nsample*(ftmp*(1-af)+(1-ftmp)*(1-af)*(1-af)));
		}
		int[][] d ={{10,0,10},
				{20,0,10},
				{5,5,3}
		};
		double f = EstF(d2);
		System.out.println(f);
		double[] like ={1,5,1,0,4};
		int minID=MinFive(like);
		System.out.print(minID);

	}
	public static double EstF(int[][] d){
		double ret = 0;
		
		double[][] cal=new double[d.length][5];
		double minF =-1;
		/*
		 * [0] は有効か否か
		 * [1] は合計
		 * [2] はアレル頻度１
		 * [3] はアレル頻度２
		 * [4] はF
		 * 
		 */
		for(int i=0;i<d.length;i++){
			System.out.print(d[i][0]+"\t"+d[i][1]+"\t"+d[i][2]+"\t");
			int tmpsum=d[i][0]+d[i][1]+d[i][2];
			int tmpa1 =d[i][0]*2+d[i][1];
			int tmpa2 =d[i][2]*2+d[i][1];
			System.out.print(tmpsum+"\t"+tmpa1+"\t"+tmpa2+"\t");
			
			if(tmpsum==0 || tmpa1==0 || tmpa2==0){
				
			}else{
				cal[i][0]=1;
				cal[i][1]=tmpsum;
				cal[i][2]=(double)(tmpa1)/(double)(tmpsum*2);
				cal[i][3]=1-cal[i][2];
				double tmpMinF =-cal[i][3]/cal[i][2];
				if(cal[i][2]<cal[i][3]){
					tmpMinF=-cal[i][2]/cal[i][3];
				}
				cal[i][4]=1-d[i][1]/(2*cal[i][2]*cal[i][3]*cal[i][1]);
				System.out.print(cal[i][0]+"\t"+cal[i][1]+"\t"+cal[i][2]+"\t"+cal[i][3]+"\t"+cal[i][4]+"\t"+tmpMinF+"\n");
				
				if(tmpMinF>minF){
					minF=tmpMinF;
				}
				
			}
		}
		//ret =minF;
		minF+=GETA;
		double[] five={minF,0,(minF+1)/2,0,1};
		five[1]=(five[0]+five[2])/2;
		five[3]=(five[2]+five[4])/2;
		
		double[] like = new double[5];
		for(int i=0;i<like.length;i++){
			like[i]=Like(d,cal,five[i]);
		}
		
		while(five[4]-five[0]>THRES){
			Renew(d,cal,five,like);
			System.out.println(five[0]+"\t"+five[1]+"\t"+five[2]+"\t"+five[3]+"\t"+five[4]);
			System.out.println(like[0]+"\t"+like[1]+"\t"+like[2]+"\t"+like[3]+"\t"+like[4]);
		}
		
		return five[2];
	}
	
	public static double Like(int[][] d,double[][] cal,double f){
		double ret = 0;
		for(int i=0;i<cal.length;i++){
			if(cal[i][0]==1){
				double[] genfreq={
						Math.log(f*cal[i][2]+(1-f)*cal[i][2]*cal[i][2]),
						Math.log(2*(1-f)*cal[i][2]*cal[i][3]),
						Math.log(f*cal[i][3]+(1-f)*cal[i][3]*cal[i][3])
				};
				ret +=d[i][0]*genfreq[0]+d[i][1]*genfreq[1]+d[i][2]*genfreq[2];
			}
		}
		
		return ret;
	}
	
	public static void Renew(int[][] d, double[][] cal, double[] f,double[] like){
		int minID=MaxFive(like);
		if(minID==0){
			f[4]=f[2];
			f[2]=f[1];
			f[1]=(f[0]+f[2])/2;
			f[3]=(f[2]+f[4])/2;
			like[4]=like[2];
			like[2]=like[1];
			like[1]=Like(d,cal,f[1]);
			like[3]=Like(d,cal,f[3]);
		}else if(minID==1){
			f[4]=f[2];
			f[2]=f[1];
			f[1]=(f[0]+f[2])/2;
			f[3]=(f[2]+f[4])/2;
			like[4]=like[2];
			like[2]=like[1];
			like[1]=Like(d,cal,f[1]);
			like[3]=Like(d,cal,f[3]);
		}else if(minID==2){
			f[0]=f[1];
			f[4]=f[3];
			f[1]=(f[0]+f[2])/2;
			f[3]=(f[2]+f[4])/2;
			like[4]=like[3];
			like[0]=like[1];
			like[1]=Like(d,cal,f[1]);
			like[3]=Like(d,cal,f[3]);
		}else if(minID==3){
			f[0]=f[2];
			f[2]=f[3];
			f[1]=(f[0]+f[2])/2;
			f[3]=(f[2]+f[4])/2;
			like[2]=like[3];
			like[0]=like[2];
			like[1]=Like(d,cal,f[1]);
			like[3]=Like(d,cal,f[3]);
		}else if(minID==4){
			f[0]=f[2];
			f[2]=f[3];
			f[1]=(f[0]+f[2])/2;
			f[3]=(f[2]+f[4])/2;
			like[2]=like[3];
			like[0]=like[2];
			like[1]=Like(d,cal,f[1]);
			like[3]=Like(d,cal,f[3]);
		}
	}
	
	public static int MinFive(double[] like){
		
		int minID=5;
		if(like[1]<like[0]){
			if(like[2]<like[1]){
				if(like[3]<like[2]){
					if(like[4]<like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]<like[2]){
						minID=4;
					}else{
						minID=2;
					}
				}
			}else{
				if(like[3]<like[1]){
					if(like[4]<like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]<like[1]){
						minID=4;
					}else{
						minID=1;
					}
				}
			}
		}else{
			if(like[2]<like[0]){
				if(like[3]<like[2]){
					if(like[4]<like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]<like[2]){
						minID=4;
					}else{
						minID=2;
					}
				}
			}else{
				if(like[3]<like[0]){
					if(like[4]<like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]<like[0]){
						minID=4;
					}else{
						minID=0;
					}
				}
			}
		}
		return minID;
	}
	public static int MaxFive(double[] like){
		
		int minID=5;
		if(like[1]>like[0]){
			if(like[2]>like[1]){
				if(like[3]>like[2]){
					if(like[4]>like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]>like[2]){
						minID=4;
					}else{
						minID=2;
					}
				}
			}else{
				if(like[3]>like[1]){
					if(like[4]>like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]>like[1]){
						minID=4;
					}else{
						minID=1;
					}
				}
			}
		}else{
			if(like[2]>like[0]){
				if(like[3]>like[2]){
					if(like[4]>like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]>like[2]){
						minID=4;
					}else{
						minID=2;
					}
				}
			}else{
				if(like[3]>like[0]){
					if(like[4]>like[3]){
						minID=4;
					}else{
						minID=3;
					}
				}else{
					if(like[4]>like[0]){
						minID=4;
					}else{
						minID=0;
					}
				}
			}
		}
		return minID;
	}

}