Ewens sampling formulaを遺伝現象的に説明し直す

• Ewens sampling formulaは以下の式で表され、
  • [tex:Pr*1=\frac{n!\theta^k}{\theta^{\[n\]}}\prod_{i=1}^n \frac{1}{i^{m_i}m_i!}],
• このような式であらわされるような突然変異の係数について、サンプル数nのときの、異なるアレルの数の期待値は[tex:\sum length(m_1,m_2,...,m_n)\times Pr*2],であるが、これには別の簡便な式表現があって
  • 
• この遺伝現象に則した解説はこちらがわかりやすい
• Ewens sampling formulaによる、確率のべたな計算(すべての分割の確率の和はもちろん１になる）と、それに基づく分割の長さ()の期待値が、簡便な式と一致することを確かめるためのjavaソールはこちら。
• 使い方と出力例

		int n=5;
		double theta=0.0001;
		double alpha=0;

		Ewens.Ewens_All_Pr(n,theta,alpha);

		double ExpectedEwens=Ewens.ExpectedEwens(n,theta);
		System.out.println("exp(k)\t"+ExpectedEwens);

はじめに分割、次に和の構成
0	0	0	0	1		0	0	0	0	5	0.9997916954827282
1	0	0	1	0		0	0	0	1	4	1.2497396193534113E-4
0	1	1	0	0		0	0	0	2	3	8.331597462356081E-5
2	0	1	0	0		0	0	1	1	3	8.331597462356058E-9
1	2	0	0	0		0	0	1	2	2	6.2486980967670506E-9
3	1	0	0	0		0	1	1	1	2	4.1657987311780475E-13
5	0	0	0	0		1	1	1	1	1	4.1657987311780685E-18
expected	1.0002083190983988
exp(k)	1.0002083190983997

public static double ExpectedEwens(int n,double theta){
		double ret=0;
		for(int i=0;i<n;i++){
			ret+=theta/(theta+i);
		}

		return ret;
	}

	public static void Ewens_All_Pr(int n,double theta,double alpha){
		/*
		 * n-1本の区切りをn+1箇所に入れる。
		 * すべての区切りが同じ箇所に入っても可。
		 * ただし、前半の線分は後半の線分より長くてはいけない
		 * alpha=0はEwens
		 * さらにtheta=1はuniformly distributed random permutationからの整数分割
		 */
		int[] kugiri=new int[n-1];



		boolean loop=true;
		int count=0;
		System.out.println("はじめに分割、次に和の構成");
		double expected=0;
		while(loop){


			boolean tmpbool=true;
			for(int i=0;i<kugiri.length;i++){
				if(kugiri[i]<n){
					tmpbool=false;
					i=kugiri.length;
				}
			}
			if(tmpbool){
				loop=false;
			}
			int[] lengthVector=lengthFromKugiri(kugiri,n);

			boolean lengthOK=CheckLength(lengthVector);
			if(lengthOK){
				count++;

				int[] partition=fromLengthToPartition(lengthVector);
				int k=0;

				for(int i=0;i<partition.length;i++){
					System.out.print(partition[i]+"\t");
					k+=partition[i];
				}
				System.out.print("\t");
				for(int i=0;i<lengthVector.length;i++){
					System.out.print(lengthVector[i]+"\t");
				}
				double pr=Pitman(partition,theta,alpha);
				expected+=pr*k;
				System.out.print(pr);
				System.out.println();
			}

			kugiri=update(kugiri,n);
		}
		//System.out.println("count\t"+count);
		System.out.println("expected\t"+expected);
	}

*1:M_1,M_2,....,M_n)=(m_1,m_2,...,m_n

*2:M_1,M_2,....,M_n)=(m_1,m_2,...,m_n