2006-07-28

第４章（２）One-way ANOVA　多標本１変数パーミュテーションテストの例

C>=２群のすべてが等しいか、そうでないかの検定
全nデータが、 $n_1,n_2,...n_C$ 個ずつC群に分かれている。
パーミュテーションは $￥frac{n!}{n_1!n_2!...n_C!}$ 通り
統計量としては、 $S￥sum_{j=1}^C(￥bar{Y_j}-￥bar{Y_{all}})^2￥times n_j$ がわかりやすいが、それとpermutationally にequivalentな統計量は $T=￥sum_{j=1}^C(￥bar{Y_j})^2 ￥times n_j$ がある。また、one-way ANOVAで用いるFもT,Sと１対１の単調な関係にある。他方、ランクテストの統計量である、Kruskal-Wallisもあるが、これは、permutationally にequivalentではない。それを示したのが掲載図とこの図(リンク)である。
これまでと同様に、教科書中のデータサンプルについて、実行し、掲載図を作るにあたって用いたJavaソースはこちら



public class oneWayANOVAunivarMultiSamplePerm {
	public static void main(String[] args) {

		

//		異なるdelta値での統計量を算出

		int numseg = 1;

		double[] delta = new double[numseg];

		

		double initial=0;

		double dif = -initial*2/(double)numseg;

		for(int i=0;i<delta.length;i++){

			delta[i] = initial+dif*i;

		}

		//double[] delta = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};

		//double[] delta = {0};

		//for(int i=0;i<delta.length;i++){

			//delta[i] /= 0.1;

		//}

		//Data

		

		double[][] data = new double[3][];

		double[] tmp1 = {10.2,8.2,8.9,8.0,8.3,8.0};

		double[] tmp2 = {12.2,10.6,9.9,13.0,8.1,10.8,11.5};

		double[] tmp3 = {9.2,10.5,9.2,8.7,9.0};

		

		//double[] tmp1 = {3.42,3.84,3.96,3.76};

		//double[] tmp2 = {3.17,3.63,3.47,3.44,3.39};

		//double[] tmp3 = {3.64,3.72,3.91};

		/*

		double[][] data = new double[4][];

		double[] tmp1 = {9,9.5,9.75,10,13,9.5};

		double[] tmp2 = {11,10,10,11.75,10.5,15};

		double[] tmp3 = {11.5,12,9,11.5,13.25,13};

		double[] tmp4 = {13.25,11.5,12,13.5,11.5,12.35};

		*/

		data[0] = MiscUtil.DeepCopyDouble1(tmp1);

		data[1] = MiscUtil.DeepCopyDouble1(tmp2);

		data[2] = MiscUtil.DeepCopyDouble1(tmp3);

		//data[3] = MiscUtil.DeepCopyDouble1(tmp4);

		int[] categoryNum = new int[data.length];

		for(int i=0;i<categoryNum.length;i++){

			categoryNum[i] = data[i].length;

		}

		

		

		double[] concat = concatMult(data);

		

		//double[] x1 = {1,2,3,4,5};

		//double[] x2 = {11,12,13,14,15};

		//double[] x1 = {66,57,81,62,61,60,73,59,80,55,67,70};

		//double[] x2 = {64,58,45,43,37,56,44,42};

		

		//double[] concat = concat(x1,x2);

		int c = 100000000;

		

		int seed = (int)(Math.random()*c);

		

		//int seed = 10000;

		String mval ="";

		String mval2 ="";

		for(int w=0;w<delta.length;w++){

			MersenneTwisterFast mt = new MersenneTwisterFast(seed);

			//double[] concatMu = plusMu(concat,delta[w],x1.length);		

			double[] concatMu = MiscUtil.DeepCopyDouble1(concat);

			

			double kw = KruskalWallis(data);

			//System.out.println("mw " + mw);

			

			double DEV= 0.00000001;

			int perm =10000;

			

			//double to = sumUpto(concatMu,x1.length);

			double to = tAnova(data);

			double to2 = sAnova(data);

			double to3 = fRatio(data);

			System.out.println("fRatio " + to3);

			System.out.println("to to2 " + to + " " + to2);

			//System.out.println("mu " + delta[w]);

			//for(int i=0;i<concat.length;i++){

				//System.out.println("concat concatMu " + concat[i] + " " + concatMu[i]);

			//}

			//double[][] tmp0 = segment(concatMu,x1.length);

			//double to2 = difMean(tmp0[0],tmp0[1]);

			

			double[] t = new double[perm];

			double[] t2 = new double[perm];

			double[] fRatio = new double[perm];

			double[] kwt = new double[perm];

			for(int i=0;i<perm;i++){

				//shuffleDoubleMZ(x1,mt);

				//shuffleDoubleMZ(x2,mt);

				

				shuffleDoubleMZ(concatMu,mt);

				//concatMu = plusMu(concat,delta[w],x1.length);

				//concat = concat(x1,x2);

				//double[][] tmp = segment(concatMu,x1.length);

				double[][] tmp = segmentMult(concatMu,categoryNum);

				//t[i] = sumUpto(concatMu,x1.length);

				//t2[i] = difMean(tmp[0],tmp[1]);

				t[i]=tAnova(tmp);

				t2[i]=sAnova(tmp);

				fRatio[i]=fRatio(tmp);

				kwt[i]=KruskalWallis(tmp);

				//mwt[i] = MannWhitney(concat,x1.length);

				System.out.println("t\tmwt\t" + t[i] +"\t" + t2[i] + "\t" + kwt[i] + "\t" + fRatio[i]);

			}

			

//			度数分布・信頼区間など

			int[] h=signPerm.hist(t,20);

			double[][] hdouble = signPerm.histAccum(t,1000);

//			統計量の信頼区間を算出表示

			double confA = 0.000001;

			double[] confInterval = signPerm.confInterval(hdouble,confA);

			

			System.out.println("confident interval\t" + confInterval[0] + "\t" + confInterval[1]);

			

//			度数分布を標準出力

//			System.out.println(histSt);

			System.out.println("Histogram\n