静岡理工科大学 菅沼ホーム 目次 索引

分散分析

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,分散分析を行うためのものです.データの入力方法に関しては,プログラムの最後に提示してある入力例に対する説明を参考にして下さい.

  1. C++

    #include <stdio.h>
    
    double normal(double, double *);
    double p_normal(int *);
    double bisection(double(*)(double), double, double, double, double, int, int *);
    double normal_f(double);
    double f(double, double *, int, int);
    double p_f(int *);
    double gamma(double, int *);
    double newton(double(*)(double), double(*)(double), double, double,
                  double, int, int *);
    double f_f(double);
    double f_df(double);
    
    double p;         // α%値を計算するとき時α/100を設定
    int dof1, dof2;   // 自由度
    
    /**************************************/
    /* 分散分析(一元配置法と二元配置法) */
    /*      method : 1 or 2               */
    /*      Np[0] : 因子1の水準数         */
    /*        [1] : 因子2の水準数         */
    /*      N : データ数                  */
    /*      name[0] : 因子1の名前         */
    /*          [1] : 因子2の名前         */
    /*      a : a %値                    */
    /*      x : データ                    */
    /*           coded by Y.Suganuma      */
    /**************************************/
    void aov(int method, int *Np, int N, char name[][100], int a, double ***x)
    {
    	double *xi, *xj, **xij, xa, SP, SQ, SI, SE, VP, VQ, VI, VE, FP, FQ, FI;
    	int i1, i2, i3, sw;
    					// 一元配置法
    	if (method == 1) {
    
    		xi = new double [Np[0]];
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			xi[i1] = 0.0;
    			for (i2 = 0; i2 < N; i2++)
    				xi[i1] += x[i1][0][i2];
    			xi[i1] /= N;
    		}
    
    		xa = 0.0;
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			for (i2 = 0; i2 < N; i2++)
    				xa += x[i1][0][i2];
    		}
    		xa /= (Np[0] * N);
    
    		SP = 0.0;
    		for (i1 = 0; i1 < Np[0]; i1++)
    			SP += (xi[i1] - xa) * (xi[i1] - xa);
    		SP *= N;
    
    		SE = 0.0;
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			for (i2 = 0; i2 < N; i2++)
    				SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
    		}
    
    		VP = SP / (Np[0] - 1);
    		VE = SE / (Np[0] * (N - 1));
    
    		FP = VP / VE;
    
    		p    = 0.01 * a;
    		dof2 = Np[0] * (N - 1);
    
    		printf("全変動: 平方和 %.2f 自由度 %d\n", SP+SE, Np[0]*N-1);
    		dof1 = Np[0] - 1;
    		printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, p_f(&sw));
    		printf("水準内: 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*(N-1), VE);
    	}
    					// 二元配置法
    	else {
    
    		xi = new double [Np[0]];
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			xi[i1] = 0.0;
    			for (i2 = 0; i2 < Np[1]; i2++) {
    				for (i3 = 0; i3 < N; i3++)
    					xi[i1] += x[i1][i2][i3];
    			}
    			xi[i1] /= (Np[1] * N);
    		}
    
    		xj = new double [Np[1]];
    		for (i1 = 0; i1 < Np[1]; i1++) {
    			xj[i1] = 0.0;
    			for (i2 = 0; i2 < Np[0]; i2++) {
    				for (i3 = 0; i3 < N; i3++)
    					xj[i1] += x[i2][i1][i3];
    			}
    			xj[i1] /= (Np[0] * N);
    		}
    
    		xij = new double * [Np[0]];
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			xij[i1] = new double [Np[1]];
    			for (i2 = 0; i2 < Np[1]; i2++) {
    				xij[i1][i2] = 0.0;
    				for (i3 = 0; i3 < N; i3++)
    					xij[i1][i2] += x[i1][i2][i3];
    				xij[i1][i2] /= N;
    			}
    		}
    
    		xa = 0.0;
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			for (i2 = 0; i2 < Np[1]; i2++) {
    				for (i3 = 0; i3 < N; i3++)
    					xa += x[i1][i2][i3];
    			}
    		}
    		xa /= (Np[0] * Np[1] * N);
    
    		SP = 0.0;
    		for (i1 = 0; i1 < Np[0]; i1++)
    			SP += (xi[i1] - xa) * (xi[i1] - xa);
    		SP *= (Np[1] * N);
    
    		SQ = 0.0;
    		for (i1 = 0; i1 < Np[1]; i1++)
    			SQ += (xj[i1] - xa) * (xj[i1] - xa);
    		SQ *= (Np[0] * N);
    
    		SI = 0.0;
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			for (i2 = 0; i2 < Np[1]; i2++)
    				SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
    		}
    		SI *= N;
    
    		SE = 0.0;
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			for (i2 = 0; i2 < Np[1]; i2++) {
    				for (i3 = 0; i3 < N; i3++)
    					SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
    			}
    		}
    
    		VP = SP / (Np[0] - 1);
    		VQ = SQ / (Np[1] - 1);
    		VI = SI / ((Np[0] - 1) * (Np[1] - 1));
    		VE = SE / (Np[0] * Np[1] * (N - 1));
    
    		FP = VP / VE;
    		FQ = VQ / VE;
    		FI = VI / VE;
    
    		p    = 0.01 * a;
    		dof2 = Np[0] * Np[1] * (N - 1);
    
    		printf("全変動: 平方和 %.2f 自由度 %d\n", SP+SQ+SI+SE, Np[0]*Np[1]*N-1);
    		dof1 = Np[0] - 1;
    		printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, p_f(&sw));
    		dof1 = Np[1] - 1;
    		printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", name[1], SQ, Np[1]-1, VQ, FQ, a, p_f(&sw));
    		dof1 = (Np[0] - 1) * (Np[1] - 1);
    		printf("相互作用: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", SI, (Np[0]-1)*(Np[1]-1), VI, FI, a, p_f(&sw));
    		printf("水準内: 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*Np[1]*(N-1), VE);
    	}
    }
    
    int main()
    {
    	double ***x;
    	int i1, i2, i3, method, N, Np[2] = {1, 1}, a;
    	char name[2][100];
    
    	scanf("%d %d %d", &method, &N, &a);
    	if (method == 1 || method == 2) {
    		for (i1 = 0; i1 < method; i1++)
    			scanf("%s %d", name[i1], &Np[i1]);
    		x = new double ** [Np[0]];
    		for (i1 = 0; i1 < Np[0]; i1++) {
    			x[i1] = new double * [Np[1]];
    			for (i2 = 0; i2 < Np[1]; i2++)
    				x[i1][i2] = new double [N];
    		}
    		for (i1 = 0; i1 < Np[1]; i1++) {
    			for (i2 = 0; i2 < N; i2++) {
    				for (i3 = 0; i3 < Np[0]; i3++)
    					scanf("%lf", &x[i3][i1][i2]);
    			}
    		}
    		aov(method, Np, N, name, a, x);
    	}
    	else
    		printf("一元配置法,または,二元配置法だけです!\n");
    
    	return 0;
    }
    
    /****************************************/
    /* f分布の計算(P(X = ff), P(X < ff)) */
    /*      dd : P(X = ff)                  */
    /*      df1,df2 : 自由度                */
    /*      return : P(X < ff)              */
    /****************************************/
    #include <math.h>
    
    double f(double ff, double *dd, int df1, int df2)
    {
    	double pi = 4.0 * atan(1.0);
    	double pp, u, x;
    	int ia, ib, i1;
    
    	if (ff < 1.0e-10)
    		ff = 1.0e-10;
    
    	x  = ff * df1 / (ff * df1 + df2);
    
    	if(df1%2 == 0) {
    		if (df2%2 == 0) {
    			u  = x * (1.0 - x);
    			pp = x;
    			ia = 2;
    			ib = 2;
    		}
    		else {
    			u  = 0.5 * x * sqrt(1.0-x);
    			pp = 1.0 - sqrt(1.0-x);
    			ia = 2;
    			ib = 1;
    		}
    	}
    
    	else {
    		if (df2%2 == 0) {
    			u  = 0.5 * sqrt(x) * (1.0 - x);
    			pp = sqrt(x);
    			ia = 1;
    			ib = 2;
    		}
    		else {
    			u  = sqrt(x*(1.0-x)) / pi;
    			pp = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi;
    			ia = 1;
    			ib = 1;
    		}
    	}
    
    	if (ia != df1) {
    		for (i1 = ia; i1 <= df1-2; i1 += 2) {
    			pp -= 2.0 * u / i1;
    			u  *= x * (i1 + ib) / i1;
    		}
    	}
    
    	if (ib != df2) {
    		for (i1 = ib; i1 <= df2-2; i1 += 2) {
    			pp += 2.0 * u / i1;
    			u  *= (1.0 - x) * (i1 + df1) / i1;
    		}
    	}
    
    	*dd = u / ff;
    
    	return pp;
    }
    
    /****************************************/
    /* f分布のp%値(P(X > u) = 0.01p)   */
    /*      ind : >= 0 : normal(収束回数) */
    /*            = -1 : 収束しなかった     */
    /****************************************/
    #include <math.h>
    
    #define MAX 340
    
    double p_f(int *ind)
    {
    	double a, a1, b, b1, df1, df2, e, ff = 0.0, f0, x, y1, y2, yq;
    	int sw = 0;
    /*
         初期値計算の準備
              while文は,大きな自由度によるガンマ関数の
              オーバーフローを避けるため
    */
    	while (sw >= 0) {
    
    		df1 = 0.5 * (dof1 - sw);
    		df2 = 0.5 * dof2;
    		a   = 2.0 / (9.0 * (dof1 - sw));
    		a1  = 1.0 - a;
    		b   = 2.0 / (9.0 * dof2);
    		b1  = 1.0 - b;
    
    		yq  = p_normal(ind);
    
    		e   = b1 * b1 - b * yq * yq;
    
    		if (e > 0.8 || (dof1+dof2-sw) <= MAX)
    			sw = -1;
    		else {
    			sw += 1;
    			if ((dof1-sw) == 0)
    				sw = -2;
    		}
    	}
    
    	if (sw == -2)
    		*ind = -1;
    
    	else {
    /*
         f0 : 初期値
    */
    		if (e > 0.8) {
    			x  = (a1 * b1 + yq * sqrt(a1*a1*b+a*e)) / e;
    			f0 = pow(x, 3.0);
    		}
    		else {
    			y1 = pow((double)dof2, df2-1.0);
    			y2 = pow((double)dof1, df2);
    			x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
                     2.0 * y1 / y2 / p;
    			f0 = pow(x, 2.0/dof2);
    		}
    /*
         ニュートン法
    */
    		ff = newton(f_f, f_df, f0, 1.0e-6, 1.0e-10, 100, ind);
    	}
    
    	return ff;
    }
    
    /****************************************/
    /* Γ(x)の計算(ガンマ関数,近似式) */
    /*      ier : =0 : normal               */
    /*            =-1 : x=-n (n=0,1,2,・・・)  */
    /*      return : 結果                   */
    /****************************************/
    #include <math.h>
    
    double gamma(double x, int *ier)
    {
    	double err, g, s, t, v, w, y;
    	long k;
    
    	*ier = 0;
    
    	if (x > 5.0) {
    		v = 1.0 / x;
    		s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
                0.00078403922172) * v - 0.000229472093621) * v -
                0.00268132716049) * v + 0.00347222222222) * v +
                0.0833333333333) * v + 1.0;
    		g = 2.506628274631001 * exp(-x) * pow(x,x-0.5) * s;
    	}
    
    	else {
    
    		err = 1.0e-20;
    		w   = x;
    		t   = 1.0;
    
    		if (x < 1.5) {
    
    			if (x < err) {
    				k = (long)x;
    				y = (double)k - x;
    				if (fabs(y) < err || fabs(1.0-y) < err)
    					*ier = -1;
    			}
    
    			if (*ier == 0) {
    				while (w < 1.5) {
    					t /= w;
    					w += 1.0;
    				}
    			}
    		}
    
    		else {
    			if (w > 2.5) {
    				while (w > 2.5) {
    					w -= 1.0;
    					t *= w;
    				}
    			}
    		}
    
    		w -= 2.0;
    		g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
                 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
                 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
                 0.999999926;
    		g *= t;
    	}
    
    	return g;
    }
    
    /********************************/
    /* 1.0 - p - P(X > x)(関数値) */
    /********************************/
    double f_f(double x)
    {
    	double y;
    
    	return f(x, &y, dof1, dof2) - 1.0 + p;
    }
    
    /**************************/
    /* P(X = x)(関数の微分) */
    /**************************/
    double f_df(double x)
    {
    	double y, z;
    
    	z = f(x, &y, dof1, dof2);
    
    	return y;
    }
    
    /*****************************************************/
    /* Newton法による非線形方程式(f(x)=0)の解            */
    /*      fn : f(x)を計算する関数名                    */
    /*      dfn : f(x)の微分を計算する関数名             */
    /*      x0 : 初期値                                  */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    /*      max : 最大試行回数                           */
    /*      ind : 実際の試行回数                         */
    /*            (負の時は解を得ることができなかった) */
    /*      return : 解                                  */
    /*****************************************************/
    #include <math.h>
    
    double newton(double(*f)(double), double(*df)(double), double x0,
                  double eps1, double eps2, int max, int *ind)
    {
    	double g, dg, x, x1;
    	int sw;
    
    	x1   = x0;
    	x    = x1;
    	*ind = 0;
    	sw   = 0;
    
    	while (sw == 0 && *ind >= 0) {
    
    		sw    = 1;
    		*ind += 1;
    		g     = (*f)(x1);
    
    		if (fabs(g) > eps2) {
    			if (*ind <= max) {
    				dg = (*df)(x1);
    				if (fabs(dg) > eps2) {
    					x = x1 - g / dg;
    					if (fabs(x-x1) > eps1 && fabs(x-x1) > eps1*fabs(x)) {
    						x1 = x;
    						sw = 0;
    					}
    				}
    				else
    					*ind = -1;
    			}
    			else
    				*ind = -1;
    		}
    	}
    
    	return x;
    }
    
    /*************************************************/
    /* 標準正規分布N(0,1)の計算(P(X = x), P(X < x))*/
    /*      w : P(X = x)                             */
    /*      return : P(X < x)                        */
    /*************************************************/
    #include <math.h>
    
    double normal(double x, double *w)
    {
    	double pi = 4.0 * atan(1.0);
    	double y, z, P;
    /*
         確率密度関数(定義式)
    */
    	*w = exp(-0.5 * x * x) / sqrt(2.0*pi);
    /*
         確率分布関数(近似式を使用)
    */
    	y = 0.70710678118654 * fabs(x);
    	z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
            y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
            y * 0.0000430638)))));
    	P = 1.0 - pow(z,-16.0);
    
    	if (x < 0.0)
    		P = 0.5 - 0.5 * P;
    	else
    		P = 0.5 + 0.5 * P;
    
    	return P;
    }
    
    /******************************************************************/
    /* 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用) */
    /*      ind : >= 0 : normal(収束回数)                           */
    /*            = -1 : 収束しなかった                               */
    /******************************************************************/
    double p_normal(int *ind)
    {
    	double u;
    	int sw;
    
    	u    = bisection(normal_f, -7.0, 7.0, 1.0e-6, 1.0e-10, 100, &sw);
    	*ind = sw;
    
    	return u;
    }
    
    /******************************/
    /* 1.0 - p - P(X>x)(関数値) */
    /******************************/
    double normal_f(double x)
    {
    	double y;
    
    	return 1.0 - p - normal(x, &y);
    }
    
    /*********************************************************/
    /* 二分法による非線形方程式(f(x)=0)の解                  */
    /*      fn : f(x)を計算する関数名                        */
    /*      x1,x2 : 初期値                                   */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
    /*      max : 最大試行回数                               */
    /*      ind : 実際の試行回数                             */
    /*            (負の時は解を得ることができなかった)     */
    /*      return : 解                                      */
    /*********************************************************/
    #include <math.h>
    
    double bisection(double(*f)(double), double x1, double x2,
                     double eps1, double eps2, int max, int *ind)
    {
    	double f0, f1, f2, x0 = 0.0;
    	int sw;
    
    	f1 = (*f)(x1);
    	f2 = (*f)(x2);
    
    	if (f1*f2 > 0.0)
    		*ind = -1;
    
    	else {
    		*ind = 0;
    		if (f1*f2 == 0.0)
    			x0 = (f1 == 0.0) ? x1 : x2;
    		else {
    			sw = 0;
    			while (sw == 0 && *ind >= 0) {
    				sw    = 1;
    				*ind += 1;
    				x0    = 0.5 * (x1 + x2);
    				f0    = (*f)(x0);
    
    				if (fabs(f0) > eps2) {
    					if (*ind <= max) {
    						if (fabs(x1-x2) > eps1 && fabs(x1-x2) > eps1*fabs(x2)) {
    							sw = 0;
    							if (f0*f1 < 0.0) {
    								x2 = x0;
    								f2 = f0;
    							}
    							else {
    								x1 = x0;
    								f1 = f0;
    							}
    						}
    					}
    					else
    						*ind = -1;
    				}
    			}
    		}
    	}
    
    	return x0;
    }
    
    /*
    -------- 一元配置法に対するデータ例(コメント部分を除いて下さい)--------
    1 6 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    工場 3   // 因子の名前 水準の数
    3.1 4.7 5.1   // 各水準に対する1番目のデータ
    4.1 5.6 3.7   // 各水準に対する2番目のデータ
    3.3 4.3 4.5   // 各水準に対する3番目のデータ
    3.9 5.9 6.0   // 各水準に対する4番目のデータ
    3.7 6.1 3.9   // 各水準に対する5番目のデータ
    2.4 4.2 5.4   // 各水準に対する6番目のデータ
    
    -------- 二元配置法に対するデータ例(コメント部分を除いて下さい)--------
    2 3 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    薬剤 5   // 1番目の因子の名前 その水準の数
    品種 2   // 2番目の因子の名前 その水準の数
    3 4 12 -4 -4   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する1番目のデータ
    8 -8 31 12 19   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する2番目のデータ
    7 -5 8 0 23   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する3番目のデータ
    8 -10 9 10 15   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する1番目のデータ
    -5 11 26 -1 13   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する2番目のデータ
    10 -6 13 -7 -6   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する3番目のデータ
    */
    			

  2. Java

      アプレット版では,任意のデータに対して画面上で実行することができます.
    /****************************/
    /* 分散分析                 */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.io.*;
    import java.util.StringTokenizer;
    
    public class Test {
    	public static void main(String args[]) throws IOException
    	{
    		double x[][][];
    		int i1, i2, i3, method, N, a, Np[] = new int [2];
    		String line, name[] = new String [2];
    		StringTokenizer str;
    		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
    
    		Np[0]  = 1;
    		Np[1]  = 1;
    		line   = in.readLine();
    		str    = new StringTokenizer(line, " \r\n");
    		method = Integer.parseInt(str.nextToken());
    		N      = Integer.parseInt(str.nextToken());
    		a      = Integer.parseInt(str.nextToken());
    		if (method == 1 || method == 2) {
    			for (i1 = 0; i1 < method; i1++) {
    				line     = in.readLine();
    				str      = new StringTokenizer(line, " \r\n");
    				name[i1] = str.nextToken();
    				Np[i1]   = Integer.parseInt(str.nextToken());
    			}
    			x = new double [Np[0]][Np[1]][N];
    			for (i1 = 0; i1 < Np[1]; i1++) {
    				for (i2 = 0; i2 < N; i2++) {
    					line = in.readLine();
    					str  = new StringTokenizer(line, " \r\n");
    					for (i3 = 0; i3 < Np[0]; i3++) {
    						x[i3][i1][i2] = Double.parseDouble(str.nextToken());
    					}
    				}
    			}
    			aov(method, Np, N, name, a, x);
    		}
    		else
    			System.out.println("一元配置法,または,二元配置法だけです!");
    	}
    
    	/**************************************/
    	/* 分散分析(一元配置法と二元配置法) */
    	/*      method : 1 or 2               */
    	/*      Np[0] : 因子1の水準数         */
    	/*        [1] : 因子2の水準数         */
    	/*      N : データ数                  */
    	/*      name[0] : 因子1の名前         */
    	/*          [1] : 因子2の名前         */
    	/*      a : a %値                    */
    	/*      x : データ                    */
    	/*           coded by Y.Suganuma      */
    	/**************************************/
    	static void aov(int method, int Np[], int N, String name[], int a, double x[][][])
    	{
    		double xi[], xj[], xij[][], xa, SP, SQ, SI, SE, VP, VQ, VI, VE, FP, FQ, FI, p;
    		int i1, i2, i3, sw[] = new int [1], dof1, dof2;
    					// 一元配置法
    		if (method == 1) {
    
    			xi = new double [Np[0]];
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				xi[i1] = 0.0;
    				for (i2 = 0; i2 < N; i2++)
    					xi[i1] += x[i1][0][i2];
    				xi[i1] /= N;
    			}
    
    			xa = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < N; i2++)
    					xa += x[i1][0][i2];
    			}
    			xa /= (Np[0] * N);
    
    			SP = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++)
    				SP += (xi[i1] - xa) * (xi[i1] - xa);
    			SP *= N;
    
    			SE = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < N; i2++)
    					SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
    			}
    
    			VP = SP / (Np[0] - 1);
    			VE = SE / (Np[0] * (N - 1));
    
    			FP = VP / VE;
    
    			p    = 0.01 * a;
    			dof2 = Np[0] * (N - 1);
    
    			System.out.printf("全変動: 平方和 %.2f 自由度 %d\n", SP+SE, Np[0]*N-1);
    			dof1 = Np[0] - 1;
    			F kn = new F(dof1, dof2, p);
    			System.out.printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, kn.p_f(sw));
    			System.out.printf("水準内: 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*(N-1), VE);
    		}
    					// 二元配置法
    		else {
    
    			xi = new double [Np[0]];
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				xi[i1] = 0.0;
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						xi[i1] += x[i1][i2][i3];
    				}
    				xi[i1] /= (Np[1] * N);
    			}
    
    			xj = new double [Np[1]];
    			for (i1 = 0; i1 < Np[1]; i1++) {
    				xj[i1] = 0.0;
    				for (i2 = 0; i2 < Np[0]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						xj[i1] += x[i2][i1][i3];
    				}
    				xj[i1] /= (Np[0] * N);
    			}
    
    			xij = new double [Np[0]][Np[1]];
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					xij[i1][i2] = 0.0;
    					for (i3 = 0; i3 < N; i3++)
    						xij[i1][i2] += x[i1][i2][i3];
    					xij[i1][i2] /= N;
    				}
    			}
    
    			xa = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						xa += x[i1][i2][i3];
    				}
    			}
    			xa /= (Np[0] * Np[1] * N);
    
    			SP = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++)
    				SP += (xi[i1] - xa) * (xi[i1] - xa);
    			SP *= (Np[1] * N);
    
    			SQ = 0.0;
    			for (i1 = 0; i1 < Np[1]; i1++)
    				SQ += (xj[i1] - xa) * (xj[i1] - xa);
    			SQ *= (Np[0] * N);
    
    			SI = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++)
    					SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
    			}
    			SI *= N;
    
    			SE = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
    				}
    			}
    
    			VP = SP / (Np[0] - 1);
    			VQ = SQ / (Np[1] - 1);
    			VI = SI / ((Np[0] - 1) * (Np[1] - 1));
    			VE = SE / (Np[0] * Np[1] * (N - 1));
    
    			FP = VP / VE;
    			FQ = VQ / VE;
    			FI = VI / VE;
    
    			p    = 0.01 * a;
    			dof2 = Np[0] * Np[1] * (N - 1);
    
    			System.out.printf("全変動: 平方和 %.2f 自由度 %d\n", SP+SQ+SI+SE, Np[0]*Np[1]*N-1);
    			dof1 = Np[0] - 1;
    			F kn = new F(dof1, dof2, p);
    			System.out.printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, kn.p_f(sw));
    			dof1 = Np[1] - 1;
    			kn = new F(dof1, dof2, p);
    			System.out.printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", name[1], SQ, Np[1]-1, VQ, FQ, a, kn.p_f(sw));
    			dof1 = (Np[0] - 1) * (Np[1] - 1);
    			kn = new F(dof1, dof2, p);
    			System.out.printf("相互作用: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", SI, (Np[0]-1)*(Np[1]-1), VI, FI, a, kn.p_f(sw));
    			System.out.printf("水準内: 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*Np[1]*(N-1), VE);
    		}
    	}
    }
    
    /****************/
    /* 関数値の計算 */
    /****************/
    class F extends Newton_Bisection_Gamma
    {
    	int sw;
    	int dof1, dof2;    // 自由度
    	double p;   // α%値を計算するとき時α/100を設定
    					// コンストラクタ
    	F(int dof1, int dof2, double p)
    	{
    		this.dof1 = dof1;
    		this.dof2 = dof2;
    		this.p = p;
    	}
    					// 関数値の計算
    	double snx(double x)
    	{
    		double y = 0.0, w[] = new double [1];
    
    		switch (sw) {
    						// 関数値(f(x))の計算(正規分布)
    			case 0:
    				y = 1.0 - p - normal(x, w);
    				break;
    						// 関数値(f(x))の計算(F分布)
    			case 1:
    				y = f(x, w) - 1.0 + p;
    				break;
    		}
    
    		return y;
    	}
    					// 関数の微分の計算(F分布)
    	double dsnx(double x)
    	{
    		double y = 0.0, w[] = new double [1];
    		y = f(x, w);
    		y = w[0];
    		return y;
    	}
    
    	/****************************************/
    	/* f分布の計算(P(X = ff), P(X < ff)) */
    	/*      dd : P(X = ff)                  */
    	/*      return : P(X < ff)              */
    	/****************************************/
    	double f(double ff, double dd[])
    	{
    		double pi = Math.PI;
    		double pp, u, x;
    		int ia, ib, i1;
    
    		if (ff < 1.0e-10)
    			ff = 1.0e-10;
    
    		x  = ff * dof1 / (ff * dof1 + dof2);
    
    		if(dof1%2 == 0) {
    			if (dof2%2 == 0) {
    				u  = x * (1.0 - x);
    				pp = x;
    				ia = 2;
    				ib = 2;
    			}
    			else {
    				u  = 0.5 * x * Math.sqrt(1.0-x);
    				pp = 1.0 - Math.sqrt(1.0-x);
    				ia = 2;
    				ib = 1;
    			}
    		}
    
    		else {
    			if (dof2%2 == 0) {
    				u  = 0.5 * Math.sqrt(x) * (1.0 - x);
    				pp = Math.sqrt(x);
    				ia = 1;
    				ib = 2;
    			}
    			else {
    				u  = Math.sqrt(x*(1.0-x)) / pi;
    				pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
    				ia = 1;
    				ib = 1;
    			}
    		}
    
    		if (ia != dof1) {
    			for (i1 = ia; i1 <= dof1-2; i1 += 2) {
    				pp -= 2.0 * u / i1;
    				u  *= x * (i1 + ib) / i1;
    			}
    		}
    
    		if (ib != dof2) {
    			for (i1 = ib; i1 <= dof2-2; i1 += 2) {
    				pp += 2.0 * u / i1;
    				u  *= (1.0 - x) * (i1 + dof1) / i1;
    			}
    		}
    
    		dd[0] = u / ff;
    
    		return pp;
    	}
    
    	/****************************************/
    	/* f分布のp%値(P(X > u) = 0.01p)   */
    	/*      ind : >= 0 : normal(収束回数) */
    	/*            = -1 : 収束しなかった     */
    	/****************************************/
    	double p_f(int ind[])
    	{
    		double a = 0.0, a1 = 0.0, b = 0.0, b1 = 0.0, df1 = 0.0, df2 = 0.0,
                   e = 0.0, ff = 0.0, f0, x, y1, y2, yq = 0.0;
    		int sww = 0, MAX = 340;
    	/*
    	     初期値計算の準備
    	          while文は,大きな自由度によるガンマ関数の
    	          オーバーフローを避けるため
    	*/
    		while (sww >= 0) {
    
    			df1 = 0.5 * (dof1 - sww);
    			df2 = 0.5 * dof2;
    			a   = 2.0 / (9.0 * (dof1 - sww));
    			a1  = 1.0 - a;
    			b   = 2.0 / (9.0 * dof2);
    			b1  = 1.0 - b;
    
    			yq  = p_normal(ind);
    
    			e   = b1 * b1 - b * yq * yq;
    
    			if (e > 0.8 || (dof1+dof2-sww) <= MAX)
    				sww = -1;
    			else {
    				sww += 1;
    				if ((dof1-sww) == 0)
    					sww = -2;
    			}
    		}
    
    		if (sww == -2)
    			ind[0] = -1;
    	/*
    	     f0 : 初期値
    	*/
    		else {
    			if (e > 0.8) {
    				x  = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e;
    				f0 = Math.pow(x, 3.0);
    			}
    			else {
    				y1 = Math.pow((double)dof2, df2-1.0);
    				y2 = Math.pow((double)dof1, df2);
    				x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
                         2.0 * y1 / y2 / p;
    				f0 = Math.pow(x, 2.0/dof2);
    			}
    	/*
    	     ニュートン法
    	*/
    			sw = 1;
    			ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind);
    		}
    
    		return ff;
    	}
    
    	/*************************************************/
    	/* 標準正規分布N(0,1)の計算(P(X = x), P(X < x))*/
    	/*      w : P(X = x)                             */
    	/*      return : P(X < x)                        */
    	/*************************************************/
    	double normal(double x, double w[])
    	{
    		double y, z, P;
    	/*
    	     確率密度関数(定義式)
    	*/
    		w[0] = Math.exp(-0.5 * x * x) / Math.sqrt(2.0*Math.PI);
    	/*
    	     確率分布関数(近似式を使用)
    	*/
    		y = 0.70710678118654 * Math.abs(x);
    		z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
                y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
                y * 0.0000430638)))));
    		P = 1.0 - Math.pow(z, -16.0);
    
    		if (x < 0.0)
    			P = 0.5 - 0.5 * P;
    		else
    			P = 0.5 + 0.5 * P;
    
    		return P;
    	}
    
    	/******************************************************************/
    	/* 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用) */
    	/*      ind : >= 0 : normal(収束回数)                           */
    	/*            = -1 : 収束しなかった                               */
    	/******************************************************************/
    	double p_normal(int ind[])
    	{
    		double u;
    		int sww[] = new int [1];
    
    		sw     = 0;
    		u      = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, sww);
    		ind[0] = sww[0];
    
    		return u;
    	}
    }
    
    abstract class Newton_Bisection_Gamma
    {
    	abstract double snx(double x);
    	abstract double dsnx(double x);
    
    	/****************************************/
    	/* Γ(x)の計算(ガンマ関数,近似式) */
    	/*      ier : =0 : normal               */
    	/*            =-1 : x=-n (n=0,1,2,・・・)  */
    	/*      return : 結果                   */
    	/****************************************/
    	double gamma(double x, int ier[])
    	{
    		double err, g, s, t, v, w, y;
    		int k;
    
    		ier[0] = 0;
    
    		if (x > 5.0) {
    			v = 1.0 / x;
    			s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
                    0.00078403922172) * v - 0.000229472093621) * v -
                    0.00268132716049) * v + 0.00347222222222) * v +
                    0.0833333333333) * v + 1.0;
    			g = 2.506628274631001 * Math.exp(-x) * Math.pow(x,x-0.5) * s;
    		}
    
    		else {
    
    			err = 1.0e-20;
    			w   = x;
    			t   = 1.0;
    
    			if (x < 1.5) {
    
    				if (x < err) {
    					k = (int)x;
    					y = (double)k - x;
    					if (Math.abs(y) < err || Math.abs(1.0-y) < err)
    						ier[0] = -1;
    				}
    
    				if (ier[0] == 0) {
    					while (w < 1.5) {
    						t /= w;
    						w += 1.0;
    					}
    				}
    			}
    
    			else {
    				if (w > 2.5) {
    					while (w > 2.5) {
    						w -= 1.0;
    						t *= w;
    					}
    				}
    			}
    
    			w -= 2.0;
    			g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
                     0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
                     0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
                     0.999999926;
    			g *= t;
    		}
    
    		return g;
    	}
    
    	/*****************************************************/
    	/* Newton法による非線形方程式(f(x)=0)の解            */
    	/*      x1 : 初期値                                  */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    	/*      max : 最大試行回数                           */
    	/*      ind : 実際の試行回数                         */
    	/*            (負の時は解を得ることができなかった) */
    	/*      return : 解                                  */
    	/*****************************************************/
    	double newton(double x1, double eps1, double eps2, int max, int ind[])
    	{
    		double g, dg, x;
    		int sw;
    
    		x      = x1;
    		ind[0] = 0;
    		sw     = 0;
    
    		while (sw == 0 && ind[0] >= 0) {
    
    			ind[0]++;
    			sw = 1;
    			g  = snx(x1);
    
    			if (Math.abs(g) > eps2) {
    				if (ind[0] <= max) {
    					dg = dsnx(x1);
    					if (Math.abs(dg) > eps2) {
    						x = x1 - g / dg;
    						if (Math.abs(x-x1) > eps1 && Math.abs(x-x1) > eps1*Math.abs(x)) {
    							x1 = x;
    							sw = 0;
    						}
    					}
    					else
    						ind[0] = -1;
    				}
    				else
    					ind[0] = -1;
    			}
    		}
    
    		return x;
    	}
    
    	/*********************************************************/
    	/* 二分法による非線形方程式(f(x)=0)の解                  */
    	/*      x1,x2 : 初期値                                   */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
    	/*      max : 最大試行回数                               */
    	/*      ind : 実際の試行回数                             */
    	/*            (負の時は解を得ることができなかった)     */
    	/*      return : 解                                      */
    	/*********************************************************/
    	double bisection(double x1, double x2, double eps1, double eps2, int max, int ind[])
    	{
    		double f0, f1, f2, x0 = 0.0;
    		int sw;
    
    		f1 = snx(x1);
    		f2 = snx(x2);
    
    		if (f1*f2 > 0.0)
    			ind[0] = -1;
    
    		else {
    			ind[0] = 0;
    			if (f1*f2 == 0.0)
    				x0 = (f1 == 0.0) ? x1 : x2;
    			else {
    				sw = 0;
    				while (sw == 0 && ind[0] >= 0) {
    					sw      = 1;
    					ind[0] += 1;
    					x0      = 0.5 * (x1 + x2);
    					f0      = snx(x0);
    					if (Math.abs(f0) > eps2) {
    						if (ind[0] <= max) {
    							if (Math.abs(x1-x2) > eps1 && Math.abs(x1-x2) > eps1*Math.abs(x2)) {
    								sw = 0;
    								if (f0*f1 < 0.0) {
    									x2 = x0;
    									f2 = f0;
    								}
    								else {
    									x1 = x0;
    									f1 = f0;
    								}
    							}
    						}
    						else
    							ind[0] = -1;
    					}
    				}
    			}
    		}
    
    		return x0;
    	}
    }
    
    /*
    -------- 一元配置法に対するデータ例(コメント部分を除いて下さい)--------
    1 6 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    工場 3   // 因子の名前 水準の数
    3.1 4.7 5.1   // 各水準に対する1番目のデータ
    4.1 5.6 3.7   // 各水準に対する2番目のデータ
    3.3 4.3 4.5   // 各水準に対する3番目のデータ
    3.9 5.9 6.0   // 各水準に対する4番目のデータ
    3.7 6.1 3.9   // 各水準に対する5番目のデータ
    2.4 4.2 5.4   // 各水準に対する6番目のデータ
    
    -------- 二元配置法に対するデータ例(コメント部分を除いて下さい)--------
    2 3 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    薬剤 5   // 1番目の因子の名前 その水準の数
    品種 2   // 2番目の因子の名前 その水準の数
    3 4 12 -4 -4   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する1番目のデータ
    8 -8 31 12 19   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する2番目のデータ
    7 -5 8 0 23   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する3番目のデータ
    8 -10 9 10 15   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する1番目のデータ
    -5 11 26 -1 13   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する2番目のデータ
    10 -6 13 -7 -6   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する3番目のデータ
    */
    			
    アプレット版
    /****************************/
    /* 分散分析                 */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.awt.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.StringTokenizer;
    
    public class Var extends Applet implements ActionListener, ItemListener {
    	String str1 = "3.1 4.7 5.1\n4.1 5.6 3.7\n3.3 4.3 4.5\n3.9 5.9 6.0\n3.7 6.1 3.9\n2.4 4.2 5.4\n";
    	String str2 = "3 4 12 -4 -4\n8 -8 31 12 19\n7 -5 8 0 23\n8 -10 9 10 15\n-5 11 26 -1 13\n10 -6 13 -7 -6\n";
    	int meth = 1;
    	double p;         // α%値を計算するとき時α/100を設定
    	int dof1, dof2;   // 自由度
    	Choice ch;
    	TextField tx2, tx3, tx4, tx5, tx6, tx7;
    	Label lb;
    	TextArea ta1, ta2;
    	Button bt;
    	Font f = new Font("TimesRoman", Font.BOLD, 20);
    
    	/************/
    	/* 初期設定 */
    	/************/
    	public void init()
    	{
    					// レイアウト,背景色,フォント
    		setLayout(new BorderLayout(5, 5));
    		setBackground(new Color(225, 255, 225));
    		setFont(f);
    					// 上のパネル
    		Panel pn1 = new Panel();
    		add(pn1, BorderLayout.NORTH);
    
    		pn1.add(new Label("方法:"));
    		ch = new Choice ();
    		ch.addItem("一元配置法");
    		ch.addItem("二元配置法");
    		pn1.add(ch);
    		ch.addItemListener(this);
    
    		pn1.add(new Label("データ数:"));
    		tx2 = new TextField("6", 2);
    		pn1.add(tx2);
    
    		pn1.add(new Label("有意水準(%):"));
    		tx3 = new TextField("5", 3);
    		pn1.add(tx3);
    					// 中央のパネル(データ)
    		Panel pn2 = new Panel();
    		add(pn2, BorderLayout.CENTER);
    
    		pn2.add(new Label("因子名1と水準数:"));
    		tx4 = new TextField("因子1", 6);
    		pn2.add(tx4);
    		tx5 = new TextField("3", 2);
    		pn2.add(tx5);
    
    		lb = new Label("因子名2と水準数:");
    		lb.setVisible(false);
    		pn2.add(lb);
    		tx6 = new TextField("因子2", 6);
    		tx6.setVisible(false);
    		pn2.add(tx6);
    		tx7 = new TextField("2", 2);
    		tx7.setVisible(false);
    		pn2.add(tx7);
    
    		pn2.add(new Label("データ:"));
    		ta1 = new TextArea(str1, 10, 20);
    		pn2.add(ta1);
    
    		pn2.add(new Label(" "));
    		bt = new Button("OK");
    		bt.setBackground(Color.pink);
    		bt.addActionListener(this);
    		pn2.add(bt);
    					// 下のパネル(結果)
    		Panel pn3 = new Panel();
    		add(pn3, BorderLayout.SOUTH);
    
    		pn3.add(new Label("結果:"));
    		ta2 = new TextArea(7, 80);
    		pn3.add(ta2);
    	}
    
    	/************************/
    	/* 選択されたときの処理 */
    	/************************/
    	public void itemStateChanged(ItemEvent e)
    	{
    		if (e.getItemSelectable() == ch) {
    			meth = ch.getSelectedIndex() + 1;
    			if (meth == 1) {
    				tx2.setText("6");
    				tx5.setText("3");
    				lb.setVisible(false);
    				tx6.setVisible(false);
    				tx7.setVisible(false);
    				ta1.setText(str1);
    			}
    			else {
    				tx2.setText("3");
    				tx5.setText("5");
    				lb.setVisible(true);
    				tx6.setVisible(true);
    				tx7.setVisible(true);
    				ta1.setText(str2);
    			}
    			getParent().validate();
    		}
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    					// データ
    			int Np[] = new int [2];
    			String name[] = new String [2];
    			Np[0] = 1;
    			Np[1] = 1;
    			int N = Integer.parseInt(tx2.getText());
    			double a = Double.parseDouble(tx3.getText());
    			name[0] = tx4.getText();
    			Np[0] = Integer.parseInt(tx5.getText());
    			if (meth == 2) {
    				name[1] = tx6.getText();
    				Np[1] = Integer.parseInt(tx7.getText());
    			}
    			double x[][][] = new double [Np[0]][Np[1]][N];
    			String ss = ta1.getText();
    			StringTokenizer str = new StringTokenizer(ss, " \n");
    			for (int i1 = 0; i1 < Np[1]; i1++) {
    				for (int i2 = 0; i2 < N; i2++) {
    					for (int i3 = 0; i3 < Np[0]; i3++)
    						x[i3][i1][i2] = Double.parseDouble(str.nextToken());
    				}
    			}
    					// 計算
    			aov(meth, Np, N, name, a, x);
    		}
    	}
    
    	/****************/
    	/* 関数値の計算 */
    	/****************/
    	double snx(int sw, double x)
    	{
    
    		double y = 0.0, w[] = new double [1];
    
    		switch (sw) {
    						// 関数値(f(x))の計算(正規分布)
    			case 0:
    				y = 1.0 - p - normal(x, w);
    				break;
    						// 関数値(f(x))の計算(F分布)
    			case 1:
    				y = f(x, w, dof1, dof2) - 1.0 + p;
    				break;
    						// 関数の微分の計算(F分布)
    			case 2:
    				y = f(x, w, dof1, dof2);
    				y = w[0];
    				break;
    		}
    
    		return y;
    	}
    
    	/**************************************/
    	/* 分散分析(一元配置法と二元配置法) */
    	/*      method : 1 or 2               */
    	/*      Np[0] : 因子1の水準数         */
    	/*        [1] : 因子2の水準数         */
    	/*      N : データ数                  */
    	/*      name[0] : 因子1の名前         */
    	/*          [1] : 因子2の名前         */
    	/*      a : a %値                    */
    	/*      x : データ                    */
    	/*           coded by Y.Suganuma      */
    	/**************************************/
    	void aov(int method, int Np[], int N, String name[], double a, double x[][][])
    	{
    		double xi[], xj[], xij[][], xa, SP, SQ, SI, SE, VP, VQ, VI, VE, FP, FQ, FI;
    		int i1, i2, i3, sw[] = new int [1];
    					// 一元配置法
    		if (method == 1) {
    
    			xi = new double [Np[0]];
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				xi[i1] = 0.0;
    				for (i2 = 0; i2 < N; i2++)
    					xi[i1] += x[i1][0][i2];
    				xi[i1] /= N;
    			}
    
    			xa = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < N; i2++)
    					xa += x[i1][0][i2];
    			}
    			xa /= (Np[0] * N);
    
    			SP = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++)
    				SP += (xi[i1] - xa) * (xi[i1] - xa);
    			SP *= N;
    
    			SE = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < N; i2++)
    					SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
    			}
    
    			VP = SP / (Np[0] - 1);
    			VE = SE / (Np[0] * (N - 1));
    
    			FP = VP / VE;
    
    			p    = 0.01 * a;
    			dof2 = Np[0] * (N - 1);
    
    			ta2.setText("");
    			ta2.append("全変動: 平方和 " + String.format("%.2f",SP+SE) + " 自由度 " + String.format("%d",Np[0]*N-1) + "\n");
    			dof1 = Np[0] - 1;
    			ta2.append(String.format("%s",name[0]) + "の水準間: 平方和 " + String.format("%.2f",SP) + " 自由度 " + String.format("%d",Np[0]-1) + " 不偏分散 " + String.format("%.4f",VP) + " F " + String.format("%.2f",FP) + " " + String.format("%.1f",a) + "%値 " + String.format("%.2f",p_f(sw)) + "\n");
    			ta2.append("水準内: 平方和 " + String.format("%.2f",SE) + " 自由度 " + String.format("%d",Np[0]*(N-1)) + " 不偏分散 " + String.format("%.4f",VE) + "\n");
    		}
    					// 二元配置法
    		else {
    
    			xi = new double [Np[0]];
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				xi[i1] = 0.0;
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						xi[i1] += x[i1][i2][i3];
    				}
    				xi[i1] /= (Np[1] * N);
    			}
    
    			xj = new double [Np[1]];
    			for (i1 = 0; i1 < Np[1]; i1++) {
    				xj[i1] = 0.0;
    				for (i2 = 0; i2 < Np[0]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						xj[i1] += x[i2][i1][i3];
    				}
    				xj[i1] /= (Np[0] * N);
    			}
    
    			xij = new double [Np[0]][Np[1]];
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					xij[i1][i2] = 0.0;
    					for (i3 = 0; i3 < N; i3++)
    						xij[i1][i2] += x[i1][i2][i3];
    					xij[i1][i2] /= N;
    				}
    			}
    
    			xa = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						xa += x[i1][i2][i3];
    				}
    			}
    			xa /= (Np[0] * Np[1] * N);
    
    			SP = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++)
    				SP += (xi[i1] - xa) * (xi[i1] - xa);
    			SP *= (Np[1] * N);
    
    			SQ = 0.0;
    			for (i1 = 0; i1 < Np[1]; i1++)
    				SQ += (xj[i1] - xa) * (xj[i1] - xa);
    			SQ *= (Np[0] * N);
    
    			SI = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++)
    					SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
    			}
    			SI *= N;
    
    			SE = 0.0;
    			for (i1 = 0; i1 < Np[0]; i1++) {
    				for (i2 = 0; i2 < Np[1]; i2++) {
    					for (i3 = 0; i3 < N; i3++)
    						SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
    				}
    			}
    
    			VP = SP / (Np[0] - 1);
    			VQ = SQ / (Np[1] - 1);
    			VI = SI / ((Np[0] - 1) * (Np[1] - 1));
    			VE = SE / (Np[0] * Np[1] * (N - 1));
    
    			FP = VP / VE;
    			FQ = VQ / VE;
    			FI = VI / VE;
    
    			p    = 0.01 * a;
    			dof2 = Np[0] * Np[1] * (N - 1);
    
    			ta2.setText("");
    			ta2.append("全変動: 平方和 " + String.format("%.2f",SP+SQ+SI+SE) + " 自由度 " + String.format("%d",Np[0]*Np[1]*N-1) + "\n");
    			dof1 = Np[0] - 1;
    			ta2.append(String.format("%s",name[0]) + "の水準間: 平方和 " + String.format("%.2f",SP) + " 自由度 " + String.format("%d",Np[0]-1) + " 不偏分散 " + String.format("%.4f",VP) + " F " + String.format("%.2f",FP) + " " + String.format("%.1f",a) + "%値 " + String.format("%.2f",p_f(sw)) + "\n");
    			dof1 = Np[1] - 1;
    			ta2.append(String.format("%s",name[1]) + "の水準間: 平方和 " + String.format("%.2f",SQ) + " 自由度 " + String.format("%d",Np[1]-1) + " 不偏分散 " + String.format("%.4f",VQ) + " F " + String.format("%.2f",FQ) + " " + String.format("%.1f",a) + "%値 " + String.format("%.2f",p_f(sw)) + "\n");
    			dof1 = (Np[0] - 1) * (Np[1] - 1);
    			ta2.append("相互作用: 平方和 " + String.format("%.2f",SI) + " 自由度 " + String.format("%d",(Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + String.format("%.4f",VI) + " F " + String.format("%.2f",FI) + " " + String.format("%.1f",a) + "%値 " + String.format("%.2f",p_f(sw)) + "\n");
    			ta2.append("水準内: 平方和 " + String.format("%.2f",SE) + " 自由度 " + String.format("%d",Np[0]*Np[1]*(N-1)) + " 不偏分散 " + String.format("%.4f",VE) + "\n");
    		}
    	}
    
    	/****************************************/
    	/* f分布の計算(P(X = ff), P(X < ff)) */
    	/*      dd : P(X = ff)                  */
    	/*      df1,df2 : 自由度                */
    	/*      return : P(X < ff)              */
    	/****************************************/
    	double f(double ff, double dd[], int df1, int df2)
    	{
    		double pi = Math.PI;
    		double pp, u, x;
    		int ia, ib, i1;
    
    		if (ff < 1.0e-10)
    			ff = 1.0e-10;
    
    		x  = ff * df1 / (ff * df1 + df2);
    
    		if(df1%2 == 0) {
    			if (df2%2 == 0) {
    				u  = x * (1.0 - x);
    				pp = x;
    				ia = 2;
    				ib = 2;
    			}
    			else {
    				u  = 0.5 * x * Math.sqrt(1.0-x);
    				pp = 1.0 - Math.sqrt(1.0-x);
    				ia = 2;
    				ib = 1;
    			}
    		}
    
    		else {
    			if (df2%2 == 0) {
    				u  = 0.5 * Math.sqrt(x) * (1.0 - x);
    				pp = Math.sqrt(x);
    				ia = 1;
    				ib = 2;
    			}
    			else {
    				u  = Math.sqrt(x*(1.0-x)) / pi;
    				pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
    				ia = 1;
    				ib = 1;
    			}
    		}
    
    		if (ia != df1) {
    			for (i1 = ia; i1 <= df1-2; i1 += 2) {
    				pp -= 2.0 * u / i1;
    				u  *= x * (i1 + ib) / i1;
    			}
    		}
    
    		if (ib != df2) {
    			for (i1 = ib; i1 <= df2-2; i1 += 2) {
    				pp += 2.0 * u / i1;
    				u  *= (1.0 - x) * (i1 + df1) / i1;
    			}
    		}
    
    		dd[0] = u / ff;
    
    		return pp;
    	}
    
    	/****************************************/
    	/* f分布のp%値(P(X > u) = 0.01p)   */
    	/*      ind : >= 0 : normal(収束回数) */
    	/*            = -1 : 収束しなかった     */
    	/****************************************/
    	double p_f(int ind[])
    	{
    		double a = 0.0, a1 = 0.0, b = 0.0, b1 = 0.0, df1 = 0.0, df2 = 0.0,
                   e = 0.0, ff = 0.0, f0, x, y1, y2, yq = 0.0;
    		int sw = 0, MAX = 340;
    	/*
    	     初期値計算の準備
    	          while文は,大きな自由度によるガンマ関数の
    	          オーバーフローを避けるため
    	*/
    		while (sw >= 0) {
    
    			df1 = 0.5 * (dof1 - sw);
    			df2 = 0.5 * dof2;
    			a   = 2.0 / (9.0 * (dof1 - sw));
    			a1  = 1.0 - a;
    			b   = 2.0 / (9.0 * dof2);
    			b1  = 1.0 - b;
    
    			yq  = p_normal(ind);
    
    			e   = b1 * b1 - b * yq * yq;
    
    			if (e > 0.8 || (dof1+dof2-sw) <= MAX)
    				sw = -1;
    			else {
    				sw += 1;
    				if ((dof1-sw) == 0)
    					sw = -2;
    			}
    		}
    
    		if (sw == -2)
    			ind[0] = -1;
    	/*
    	     f0 : 初期値
    	*/
    		else {
    			if (e > 0.8) {
    				x  = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e;
    				f0 = Math.pow(x, 3.0);
    			}
    			else {
    				y1 = Math.pow((double)dof2, df2-1.0);
    				y2 = Math.pow((double)dof1, df2);
    				x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
                         2.0 * y1 / y2 / p;
    				f0 = Math.pow(x, 2.0/dof2);
    			}
    	/*
    	     ニュートン法
    	*/
    			ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind);
    		}
    
    		return ff;
    	}
    
    	/****************************************/
    	/* Γ(x)の計算(ガンマ関数,近似式) */
    	/*      ier : =0 : normal               */
    	/*            =-1 : x=-n (n=0,1,2,・・・)  */
    	/*      return : 結果                   */
    	/****************************************/
    	double gamma(double x, int ier[])
    	{
    		double err, g, s, t, v, w, y;
    		int k;
    
    		ier[0] = 0;
    
    		if (x > 5.0) {
    			v = 1.0 / x;
    			s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
                    0.00078403922172) * v - 0.000229472093621) * v -
                    0.00268132716049) * v + 0.00347222222222) * v +
                    0.0833333333333) * v + 1.0;
    			g = 2.506628274631001 * Math.exp(-x) * Math.pow(x,x-0.5) * s;
    		}
    
    		else {
    
    			err = 1.0e-20;
    			w   = x;
    			t   = 1.0;
    
    			if (x < 1.5) {
    
    				if (x < err) {
    					k = (int)x;
    					y = (double)k - x;
    					if (Math.abs(y) < err || Math.abs(1.0-y) < err)
    						ier[0] = -1;
    				}
    
    				if (ier[0] == 0) {
    					while (w < 1.5) {
    						t /= w;
    						w += 1.0;
    					}
    				}
    			}
    
    			else {
    				if (w > 2.5) {
    					while (w > 2.5) {
    						w -= 1.0;
    						t *= w;
    					}
    				}
    			}
    
    			w -= 2.0;
    			g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
                     0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
                     0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
                     0.999999926;
    			g *= t;
    		}
    
    		return g;
    	}
    
    	/*************************************************/
    	/* 標準正規分布N(0,1)の計算(P(X = x), P(X < x))*/
    	/*      w : P(X = x)                             */
    	/*      return : P(X < x)                        */
    	/*************************************************/
    	double normal(double x, double w[])
    	{
    		double y, z, P;
    	/*
    	     確率密度関数(定義式)
    	*/
    		w[0] = Math.exp(-0.5 * x * x) / Math.sqrt(2.0*Math.PI);
    	/*
    	     確率分布関数(近似式を使用)
    	*/
    		y = 0.70710678118654 * Math.abs(x);
    		z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
                y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
                y * 0.0000430638)))));
    		P = 1.0 - Math.pow(z, -16.0);
    
    		if (x < 0.0)
    			P = 0.5 - 0.5 * P;
    		else
    			P = 0.5 + 0.5 * P;
    
    		return P;
    	}
    
    	/******************************************************************/
    	/* 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用) */
    	/*      ind : >= 0 : normal(収束回数)                           */
    	/*            = -1 : 収束しなかった                               */
    	/******************************************************************/
    	double p_normal(int ind[])
    	{
    		double u;
    		int sw[] = new int [1];
    
    		u        = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, sw);
    		ind[0]   = sw[0];
    
    		return u;
    	}
    
    	/*****************************************************/
    	/* Newton法による非線形方程式(f(x)=0)の解            */
    	/*      x1 : 初期値                                  */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    	/*      max : 最大試行回数                           */
    	/*      ind : 実際の試行回数                         */
    	/*            (負の時は解を得ることができなかった) */
    	/*      return : 解                                  */
    	/*****************************************************/
    	double newton(double x1, double eps1, double eps2, int max, int ind[])
    	{
    		double g, dg, x;
    		int sw;
    
    		x      = x1;
    		ind[0] = 0;
    		sw     = 0;
    
    		while (sw == 0 && ind[0] >= 0) {
    
    			ind[0]++;
    			sw = 1;
    			g  = snx(1, x1);
    
    			if (Math.abs(g) > eps2) {
    				if (ind[0] <= max) {
    					dg = snx(2, x1);
    					if (Math.abs(dg) > eps2) {
    						x = x1 - g / dg;
    						if (Math.abs(x-x1) > eps1 && Math.abs(x-x1) > eps1*Math.abs(x)) {
    							x1 = x;
    							sw = 0;
    						}
    					}
    					else
    						ind[0] = -1;
    				}
    				else
    					ind[0] = -1;
    			}
    		}
    
    		return x;
    	}
    
    	/*********************************************************/
    	/* 二分法による非線形方程式(f(x)=0)の解                  */
    	/*      x1,x2 : 初期値                                   */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
    	/*      max : 最大試行回数                               */
    	/*      ind : 実際の試行回数                             */
    	/*            (負の時は解を得ることができなかった)     */
    	/*      return : 解                                      */
    	/*********************************************************/
    	double bisection(double x1, double x2, double eps1, double eps2, int max, int ind[])
    	{
    		double f0, f1, f2, x0 = 0.0;
    		int sw;
    
    		f1 = snx(0, x1);
    		f2 = snx(0, x2);
    
    		if (f1*f2 > 0.0)
    			ind[0] = -1;
    
    		else {
    			ind[0] = 0;
    			if (f1*f2 == 0.0)
    				x0 = (f1 == 0.0) ? x1 : x2;
    			else {
    				sw = 0;
    				while (sw == 0 && ind[0] >= 0) {
    					sw      = 1;
    					ind[0] += 1;
    					x0      = 0.5 * (x1 + x2);
    					f0      = snx(0, x0);
    					if (Math.abs(f0) > eps2) {
    						if (ind[0] <= max) {
    							if (Math.abs(x1-x2) > eps1 && Math.abs(x1-x2) > eps1*Math.abs(x2)) {
    								sw = 0;
    								if (f0*f1 < 0.0) {
    									x2 = x0;
    									f2 = f0;
    								}
    								else {
    									x1 = x0;
    									f1 = f0;
    								}
    							}
    						}
    						else
    							ind[0] = -1;
    					}
    				}
    			}
    		}
    
    		return x0;
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,任意のデータに対して画面上で実行することができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>分散分析</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		method = 0;
    		dof1 = 0;
    		dof2 = 0;
    		p = 0.0;
    		function main()
    		{
    						// データ
    			let Np = new Array();
    			Np[0] = 1;
    			Np[1] = 1;
    			let name = new Array();
    			let N = parseInt(document.getElementById("N").value);
    			let a = parseFloat(document.getElementById("a").value);
    			name[0] = document.getElementById("name1").value;
    			Np[0] = parseInt(document.getElementById("np_1").value);
    			let meth = method + 1;
    			if (meth == 2) {
    				name[1] = document.getElementById("name2").value;
    				Np[1] = parseInt(document.getElementById("np_2").value);
    			}
    			let x = new Array();
    			for (let i1 = 0; i1 < Np[0]; i1++) {
    				x[i1] = new Array();
    				for (let i2 = 0; i2 < Np[1]; i2++)
    					x[i1][i2] = new Array();
    			}
    			let ss = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
    			let k = 0;
    			for (let i1 = 0; i1 < Np[1]; i1++) {
    				for (let i2 = 0; i2 < N; i2++) {
    					for (let i3 = 0; i3 < Np[0]; i3++) {
    						x[i3][i1][i2] = parseFloat(ss[k]);
    						k++;
    					}
    				}
    			}
    						// 計算
    			aov(meth, Np, N, name, a, x);
    		}
    
    		/****************/
    		/* 関数値の計算 */
    		/****************/
    		function snx(sw, x)
    		{
    
    			let y = 0.0;
    			let w = new Array();
    
    			switch (sw) {
    						// 関数値(f(x))の計算(正規分布)
    				case 0:
    					y = 1.0 - p - normal(x, w);
    					break;
    						// 関数値(f(x))の計算(F分布)
    				case 1:
    					y = f(x, w, dof1, dof2) - 1.0 + p;
    					break;
    						// 関数の微分の計算(F分布)
    				case 2:
    					y = f(x, w, dof1, dof2);
    					y = w[0];
    					break;
    			}
    
    			return y;
    		}
    
    		/**************************************/
    		/* 分散分析(一元配置法と二元配置法) */
    		/*      method : 1 or 2               */
    		/*      Np[0] : 因子1の水準数         */
    		/*        [1] : 因子2の水準数         */
    		/*      N : データ数                  */
    		/*      name[0] : 因子1の名前         */
    		/*          [1] : 因子2の名前         */
    		/*      a : a %値                    */
    		/*      x : データ                    */
    		/*           coded by Y.Suganuma      */
    		/**************************************/
    		function aov(method, Np, N, name, a, x)
    		{
    			let xa;
    			let SP;
    			let SQ;
    			let SI;
    			let SE;
    			let VP;
    			let VQ;
    			let VI;
    			let VE;
    			let FP;
    			let FQ;
    			let FI;
    			let i1;
    			let i2;
    			let i3;
    			let sw = new Array();
    						// 一元配置法
    			if (method == 1) {
    
    				let xi = new Array();
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					xi[i1] = 0.0;
    					for (i2 = 0; i2 < N; i2++)
    						xi[i1] += x[i1][0][i2];
    					xi[i1] /= N;
    				}
    
    				xa = 0.0;
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					for (i2 = 0; i2 < N; i2++)
    						xa += x[i1][0][i2];
    				}
    				xa /= (Np[0] * N);
    
    				SP = 0.0;
    				for (i1 = 0; i1 < Np[0]; i1++)
    					SP += (xi[i1] - xa) * (xi[i1] - xa);
    				SP *= N;
    
    				SE = 0.0;
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					for (i2 = 0; i2 < N; i2++)
    						SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
    				}
    
    				VP = SP / (Np[0] - 1);
    				VE = SE / (Np[0] * (N - 1));
    
    				FP = VP / VE;
    
    				p    = 0.01 * a;
    				dof2 = Np[0] * (N - 1);
    
    				str = "";
    				str = str + "全変動: 平方和 " + (SP+SE) + " 自由度 " + (Np[0]*N-1) + "\n";
    				dof1 = Np[0] - 1;
    				str = str + name[0] + "の水準間: 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "%値 " + p_f(sw) + "\n";
    				str = str + "水準内: 平方和 " + SE + " 自由度 " + (Np[0]*(N-1)) + " 不偏分散 " + VE + "\n";
    				document.getElementById("ans").value = str;
    			}
    						// 二元配置法
    			else {
    
    				let xi = new Array();
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					xi[i1] = 0.0;
    					for (i2 = 0; i2 < Np[1]; i2++) {
    						for (i3 = 0; i3 < N; i3++)
    							xi[i1] += x[i1][i2][i3];
    					}
    					xi[i1] /= (Np[1] * N);
    				}
    
    				let xj = new Array();
    				for (i1 = 0; i1 < Np[1]; i1++) {
    					xj[i1] = 0.0;
    					for (i2 = 0; i2 < Np[0]; i2++) {
    						for (i3 = 0; i3 < N; i3++)
    							xj[i1] += x[i2][i1][i3];
    					}
    					xj[i1] /= (Np[0] * N);
    				}
    
    				let xij = new Array();
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					xij[i1] = new Array();
    					for (i2 = 0; i2 < Np[1]; i2++) {
    						xij[i1][i2] = 0.0;
    						for (i3 = 0; i3 < N; i3++)
    							xij[i1][i2] += x[i1][i2][i3];
    						xij[i1][i2] /= N;
    					}
    				}
    
    				xa = 0.0;
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					for (i2 = 0; i2 < Np[1]; i2++) {
    						for (i3 = 0; i3 < N; i3++)
    							xa += x[i1][i2][i3];
    					}
    				}
    				xa /= (Np[0] * Np[1] * N);
    
    				SP = 0.0;
    				for (i1 = 0; i1 < Np[0]; i1++)
    					SP += (xi[i1] - xa) * (xi[i1] - xa);
    				SP *= (Np[1] * N);
    
    				SQ = 0.0;
    				for (i1 = 0; i1 < Np[1]; i1++)
    					SQ += (xj[i1] - xa) * (xj[i1] - xa);
    				SQ *= (Np[0] * N);
    
    				SI = 0.0;
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					for (i2 = 0; i2 < Np[1]; i2++)
    						SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
    				}
    				SI *= N;
    
    				SE = 0.0;
    				for (i1 = 0; i1 < Np[0]; i1++) {
    					for (i2 = 0; i2 < Np[1]; i2++) {
    						for (i3 = 0; i3 < N; i3++)
    							SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
    					}
    				}
    
    				VP = SP / (Np[0] - 1);
    				VQ = SQ / (Np[1] - 1);
    				VI = SI / ((Np[0] - 1) * (Np[1] - 1));
    				VE = SE / (Np[0] * Np[1] * (N - 1));
    
    				FP = VP / VE;
    				FQ = VQ / VE;
    				FI = VI / VE;
    
    				p    = 0.01 * a;
    				dof2 = Np[0] * Np[1] * (N - 1);
    
    				str = "";
    				str = str + "全変動: 平方和 " + (SP+SQ+SI+SE) + " 自由度 " + (Np[0]*Np[1]*N-1) + "\n";
    				dof1 = Np[0] - 1;
    				str = str + name[0] + "の水準間: 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "%値 " + p_f(sw) + "\n";
    				dof1 = Np[1] - 1;
    				str = str + name[1] + "の水準間: 平方和 " + SQ + " 自由度 " + (Np[1]-1) + " 不偏分散 " + VQ + " F " + FQ + " " + a + "%値 " + p_f(sw) + "\n";
    				dof1 = (Np[0] - 1) * (Np[1] - 1);
    				str = str + "相互作用: 平方和 " + SI + " 自由度 " + ((Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + VI + " F " + FI + " " + a + "%値 " + p_f(sw) + "\n";
    				str = str + "水準内: 平方和 " + SE + " 自由度 " + (Np[0]*Np[1]*(N-1)) + " 不偏分散 " + VE + "\n";
    				document.getElementById("ans").value = str;
    			}
    		}
    
    		/****************************************/
    		/* f分布の計算(P(X = ff), P(X < ff)) */
    		/*      dd : P(X = ff)                  */
    		/*      df1,df2 : 自由度                */
    		/*      return : P(X < ff)              */
    		/****************************************/
    		function f(ff, dd, df1, df2)
    		{
    			let pi = Math.PI;
    			let pp;
    			let u;
    			let x;
    			let ia;
    			let ib;
    			let i1;
    
    			if (ff < 1.0e-10)
    				ff = 1.0e-10;
    
    			x  = ff * df1 / (ff * df1 + df2);
    
    			if(df1%2 == 0) {
    				if (df2%2 == 0) {
    					u  = x * (1.0 - x);
    					pp = x;
    					ia = 2;
    					ib = 2;
    				}
    				else {
    					u  = 0.5 * x * Math.sqrt(1.0-x);
    					pp = 1.0 - Math.sqrt(1.0-x);
    					ia = 2;
    					ib = 1;
    				}
    			}
    
    			else {
    				if (df2%2 == 0) {
    					u  = 0.5 * Math.sqrt(x) * (1.0 - x);
    					pp = Math.sqrt(x);
    					ia = 1;
    					ib = 2;
    				}
    				else {
    					u  = Math.sqrt(x*(1.0-x)) / pi;
    					pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
    					ia = 1;
    					ib = 1;
    				}
    			}
    
    			if (ia != df1) {
    				for (i1 = ia; i1 <= df1-2; i1 += 2) {
    					pp -= 2.0 * u / i1;
    					u  *= x * (i1 + ib) / i1;
    				}
    			}
    
    			if (ib != df2) {
    				for (i1 = ib; i1 <= df2-2; i1 += 2) {
    					pp += 2.0 * u / i1;
    					u  *= (1.0 - x) * (i1 + df1) / i1;
    				}
    			}
    
    			dd[0] = u / ff;
    
    			return pp;
    		}
    
    		/****************************************/
    		/* f分布のp%値(P(X > u) = 0.01p)   */
    		/*      ind : >= 0 : normal(収束回数) */
    		/*            = -1 : 収束しなかった     */
    		/****************************************/
    		function p_f(ind)
    		{
    			let a = 0.0;
    			let a1 = 0.0;
    			let b = 0.0;
    			let b1 = 0.0;
    			let df1 = 0.0;
    			let df2 = 0.0;
    			let e = 0.0;
    			let ff = 0.0;
    			let f0;
    			let x;
    			let y1;
    			let y2;
    			let yq = 0.0;
    			let sw = 0;
    			let MAX = 340;
    		/*
    		     初期値計算の準備
    	    	      while文は,大きな自由度によるガンマ関数の
    	        	  オーバーフローを避けるため
    		*/
    			while (sw >= 0) {
    
    				df1 = 0.5 * (dof1 - sw);
    				df2 = 0.5 * dof2;
    				a   = 2.0 / (9.0 * (dof1 - sw));
    				a1  = 1.0 - a;
    				b   = 2.0 / (9.0 * dof2);
    				b1  = 1.0 - b;
    
    				yq  = p_normal(ind);
    
    				e   = b1 * b1 - b * yq * yq;
    
    				if (e > 0.8 || (dof1+dof2-sw) <= MAX)
    					sw = -1;
    				else {
    					sw += 1;
    					if ((dof1-sw) == 0)
    						sw = -2;
    				}
    			}
    
    			if (sw == -2)
    				ind[0] = -1;
    		/*
    		     f0 : 初期値
    		*/
    			else {
    				if (e > 0.8) {
    					x  = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e;
    					f0 = Math.pow(x, 3.0);
    				}
    				else {
    					y1 = Math.pow(dof2, df2-1.0);
    					y2 = Math.pow(dof1, df2);
    					x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) * 2.0 * y1 / y2 / p;
    					f0 = Math.pow(x, 2.0/dof2);
    				}
    		/*
    		     ニュートン法
    		*/
    				ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind);
    			}
    
    			return ff;
    		}
    
    		/****************************************/
    		/* Γ(x)の計算(ガンマ関数,近似式) */
    		/*      ier : =0 : normal               */
    		/*            =-1 : x=-n (n=0,1,2,・・・)  */
    		/*      return : 結果                   */
    		/****************************************/
    		function gamma(x, ier)
    		{
    			let err;
    			let g;
    			let s;
    			let t;
    			let v;
    			let w;
    			let y;
    			let k;
    
    			ier[0] = 0;
    
    			if (x > 5.0) {
    				v = 1.0 / x;
    				s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
                    	0.00078403922172) * v - 0.000229472093621) * v -
                    	0.00268132716049) * v + 0.00347222222222) * v +
                    	0.0833333333333) * v + 1.0;
    				g = 2.506628274631001 * Math.exp(-x) * Math.pow(x,x-0.5) * s;
    			}
    
    			else {
    
    				err = 1.0e-20;
    				w   = x;
    				t   = 1.0;
    
    				if (x < 1.5) {
    
    					if (x < err) {
    						k = Math.floor(x);
    						y = k - x;
    						if (Math.abs(y) < err || Math.abs(1.0-y) < err)
    							ier[0] = -1;
    					}
    
    					if (ier[0] == 0) {
    						while (w < 1.5) {
    							t /= w;
    							w += 1.0;
    						}
    					}
    				}
    
    				else {
    					if (w > 2.5) {
    						while (w > 2.5) {
    							w -= 1.0;
    							t *= w;
    						}
    					}
    				}
    
    				w -= 2.0;
    				g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
                     	0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
                     	0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
                     	0.999999926;
    				g *= t;
    			}
    
    			return g;
    		}
    
    		/*************************************************/
    		/* 標準正規分布N(0,1)の計算(P(X = x), P(X < x))*/
    		/*      w : P(X = x)                             */
    		/*      return : P(X < x)                        */
    		/*************************************************/
    		function normal(x, w)
    		{
    			let y;
    			let z;
    			let P;
    		/*
    		     確率密度関数(定義式)
    		*/
    			w[0] = Math.exp(-0.5 * x * x) / Math.sqrt(2.0*Math.PI);
    		/*
    		     確率分布関数(近似式を使用)
    		*/
    			y = 0.70710678118654 * Math.abs(x);
    			z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
                	y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
                	y * 0.0000430638)))));
    			P = 1.0 - Math.pow(z, -16.0);
    
    			if (x < 0.0)
    				P = 0.5 - 0.5 * P;
    			else
    				P = 0.5 + 0.5 * P;
    
    			return P;
    		}
    
    		/******************************************************************/
    		/* 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用) */
    		/*      ind : >= 0 : normal(収束回数)                           */
    		/*            = -1 : 収束しなかった                               */
    		/******************************************************************/
    		function p_normal(ind)
    		{
    			let u;
    			let sw = new Array();
    
    			u        = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, sw);
    			ind[0]   = sw[0];
    
    			return u;
    		}
    
    		/*****************************************************/
    		/* Newton法による非線形方程式(f(x)=0)の解            */
    		/*      x1 : 初期値                                  */
    		/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    		/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    		/*      max : 最大試行回数                           */
    		/*      ind : 実際の試行回数                         */
    		/*            (負の時は解を得ることができなかった) */
    		/*      return : 解                                  */
    		/*****************************************************/
    		function newton(x1, eps1, eps2, max, ind)
    		{
    			let g;
    			let dg;
    			let x;
    			let sw;
    
    			x      = x1;
    			ind[0] = 0;
    			sw     = 0;
    
    			while (sw == 0 && ind[0] >= 0) {
    
    				ind[0]++;
    				sw = 1;
    				g  = snx(1, x1);
    
    				if (Math.abs(g) > eps2) {
    					if (ind[0] <= max) {
    						dg = snx(2, x1);
    						if (Math.abs(dg) > eps2) {
    							x = x1 - g / dg;
    							if (Math.abs(x-x1) > eps1 && Math.abs(x-x1) > eps1*Math.abs(x)) {
    								x1 = x;
    								sw = 0;
    							}
    						}
    						else
    							ind[0] = -1;
    					}
    					else
    						ind[0] = -1;
    				}
    			}
    
    			return x;
    		}
    
    		/*********************************************************/
    		/* 二分法による非線形方程式(f(x)=0)の解                  */
    		/*      x1,x2 : 初期値                                   */
    		/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
    		/*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
    		/*      max : 最大試行回数                               */
    		/*      ind : 実際の試行回数                             */
    		/*            (負の時は解を得ることができなかった)     */
    		/*      return : 解                                      */
    		/*********************************************************/
    		function bisection(x1, x2, eps1, eps2, max, ind)
    		{
    			let f0;
    			let f1;
    			let f2;
    			let x0 = 0.0;
    			let sw;
    
    			f1 = snx(0, x1);
    			f2 = snx(0, x2);
    
    			if (f1*f2 > 0.0)
    				ind[0] = -1;
    
    			else {
    				ind[0] = 0;
    				if (f1*f2 == 0.0)
    					x0 = (f1 == 0.0) ? x1 : x2;
    				else {
    					sw = 0;
    					while (sw == 0 && ind[0] >= 0) {
    						sw      = 1;
    						ind[0] += 1;
    						x0      = 0.5 * (x1 + x2);
    						f0      = snx(0, x0);
    						if (Math.abs(f0) > eps2) {
    							if (ind[0] <= max) {
    								if (Math.abs(x1-x2) > eps1 && Math.abs(x1-x2) > eps1*Math.abs(x2)) {
    									sw = 0;
    									if (f0*f1 < 0.0) {
    										x2 = x0;
    										f2 = f0;
    									}
    									else {
    										x1 = x0;
    										f1 = f0;
    									}
    								}
    							}
    							else
    								ind[0] = -1;
    						}
    					}
    				}
    			}
    
    			return x0;
    		}
    
    		/********/
    		/* 方法 */
    		/********/
    		function set_m()
    		{
    			let sel = document.getElementById("method");
    			for (let i1 = 0; i1 < 2; i1++) {
    				if (sel.options[i1].selected) {
    					method = i1;
    					break;
    				}
    			}
    							// 一元配置法
    			if (method == 0) {
    				document.getElementById("name2_t").style.display = "none";
    				document.getElementById("N").value = "6";
    				document.getElementById("np_1").value = "3";
    				document.getElementById("data").value = "3.1 4.7 5.1\n4.1 5.6 3.7\n3.3 4.3 4.5\n3.9 5.9 6.0\n3.7 6.1 3.9\n2.4 4.2 5.4";
    			}
    							// 二元配置法
    			else {
    				document.getElementById("name2_t").style.display = "";
    				document.getElementById("N").value = "3";
    				document.getElementById("np_1").value = "5";
    				document.getElementById("data").value = "3 4 12 -4 -4\n8 -8 31 12 19\n7 -5 8 0 23\n8 -10 9 10 15\n-5 11 26 -1 13\n10 -6 13 -7 -6";
    			}
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>分散分析</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,一元配置法及び一元配置法によって分散分析を行う場合に対する値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		方法:<SELECT ID="method" onChange="set_m()" STYLE="font-size:100%">
    			<OPTION SELECTED>一元配置法
    			<OPTION>二元配置法
    		</SELECT> 
    		データ数:<INPUT ID="N" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="6"> 
    		有意水準(%):<INPUT ID="a" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="5"><BR><BR>
    		<SPAN ID="name1_t">因子1と水準数:<INPUT ID="name1" STYLE="font-size: 100%;" TYPE="text" SIZE="15" VALUE="因子1"> <INPUT ID="np_1" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="3"> </SPAN>
    		<SPAN ID="name2_t" STYLE="display: none">因子2と水準数:<INPUT ID="name2" STYLE="font-size: 100%;" TYPE="text" SIZE="15" VALUE="因子2"> <INPUT ID="np_2" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="2"></SPAN><BR><BR>
    		データ:<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">
    3.1 4.7 5.1
    4.1 5.6 3.7
    3.3 4.3 4.5
    3.9 5.9 6.0
    3.7 6.1 3.9
    2.4 4.2 5.4
    		</TEXTAREA>
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
    		<TEXTAREA ID="ans" COLS="80" ROWS="10" STYLE="font-size: 100%;"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /**************************************/
    /* 分散分析(一元配置法と二元配置法) */
    /*           coded by Y.Suganuma      */
    /**************************************/
    
    	$Np   = array(1, 1);
    	$name = array(2);
    	$p    = 0.0;
    	$dof1 = 1;
    	$dof2 = 1;
    
    	fscanf(STDIN, "%d %d %d", $method, $N, $a);
    	if ($method == 1 || $method == 2) {
    		for ($i1 = 0; $i1 < $method; $i1++)
    			fscanf(STDIN, "%s %d", $name[$i1], $Np[$i1]);
    		$x = array($Np[0]);
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			$x[$i1] = array($Np[1]);
    			for ($i2 = 0; $i2 < $Np[1]; $i2++)
    				$x[$i1][$i2] = array($N);
    		}
    		for ($i1 = 0; $i1 < $Np[1]; $i1++) {
    			for ($i2 = 0; $i2 < $N; $i2++) {
    				$str            = trim(fgets(STDIN));
    				$x[0][$i1][$i2] = floatval(strtok($str, " "));
    				for ($i3 = 1; $i3 < $Np[0]; $i3++)
    					$x[$i3][$i1][$i2] = floatval(strtok(" "));
    			}
    		}
    		aov($method, $Np, $N, $name, $a, $x);
    	}
    	else
    		printf("一元配置法,または,二元配置法だけです!\n");
    
    /**************************************/
    /* 分散分析(一元配置法と二元配置法) */
    /*      method : 1 or 2               */
    /*      Np[0] : 因子1の水準数         */
    /*        [1] : 因子2の水準数         */
    /*      N : データ数                  */
    /*      name[0] : 因子1の名前         */
    /*          [1] : 因子2の名前         */
    /*      a : a %値                    */
    /*      x : データ                    */
    /*           coded by Y.Suganuma      */
    /**************************************/
    function aov($method, $Np, $N, $name, $a, $x)
    {
    	global $p, $dof1, $dof2;
    					// 一元配置法
    	if ($method == 1) {
    
    		$xi = array($Np[0]);
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			$xi[$i1] = 0.0;
    			for ($i2 = 0; $i2 < $N; $i2++)
    				$xi[$i1] += $x[$i1][0][$i2];
    			$xi[$i1] /= $N;
    		}
    
    		$xa = 0.0;
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			for ($i2 = 0; $i2 < $N; $i2++)
    				$xa += $x[$i1][0][$i2];
    		}
    		$xa /= ($Np[0] * $N);
    
    		$SP = 0.0;
    		for ($i1 = 0; $i1 < $Np[0]; $i1++)
    			$SP += ($xi[$i1] - $xa) * ($xi[$i1] - $xa);
    		$SP *= $N;
    
    		$SE = 0.0;
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			for ($i2 = 0; $i2 < $N; $i2++)
    				$SE += ($x[$i1][0][$i2] - $xi[$i1]) * ($x[$i1][0][$i2] - $xi[$i1]);
    		}
    
    		$VP = $SP / ($Np[0] - 1);
    		$VE = $SE / ($Np[0] * ($N - 1));
    
    		$FP = $VP / $VE;
    
    		$p    = 0.01 * $a;
    		$dof2 = $Np[0] * ($N - 1);
    
    		printf("全変動: 平方和 %.2f 自由度 %d\n", $SP+$SE, $Np[0]*$N-1);
    		$dof1 = $Np[0] - 1;
    		printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", $name[0], $SP, $Np[0]-1, $VP, $FP, $a, p_f($sw));
    		printf("水準内: 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", $SE, $Np[0]*($N-1), $VE);
    	}
    					// 二元配置法
    	else {
    
    		$xi = array($Np[0]);
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			$xi[$i1] = 0.0;
    			for ($i2 = 0; $i2 < $Np[1]; $i2++) {
    				for ($i3 = 0; $i3 < $N; $i3++)
    					$xi[$i1] += $x[$i1][$i2][$i3];
    			}
    			$xi[$i1] /= ($Np[1] * $N);
    		}
    
    		$xj = array($Np[1]);
    		for ($i1 = 0; $i1 < $Np[1]; $i1++) {
    			$xj[$i1] = 0.0;
    			for ($i2 = 0; $i2 < $Np[0]; $i2++) {
    				for ($i3 = 0; $i3 < $N; $i3++)
    					$xj[$i1] += $x[$i2][$i1][$i3];
    			}
    			$xj[$i1] /= ($Np[0] * $N);
    		}
    
    		$xij = array($Np[0]);
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			$xij[$i1] = array($Np[1]);
    			for ($i2 = 0; $i2 < $Np[1]; $i2++) {
    				$xij[$i1][$i2] = 0.0;
    				for ($i3 = 0; $i3 < $N; $i3++)
    					$xij[$i1][$i2] += $x[$i1][$i2][$i3];
    				$xij[$i1][$i2] /= $N;
    			}
    		}
    
    		$xa = 0.0;
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			for ($i2 = 0; $i2 < $Np[1]; $i2++) {
    				for ($i3 = 0; $i3 < $N; $i3++)
    					$xa += $x[$i1][$i2][$i3];
    			}
    		}
    		$xa /= ($Np[0] * $Np[1] * $N);
    
    		$SP = 0.0;
    		for ($i1 = 0; $i1 < $Np[0]; $i1++)
    			$SP += ($xi[$i1] - $xa) * ($xi[$i1] - $xa);
    		$SP *= ($Np[1] * $N);
    
    		$SQ = 0.0;
    		for ($i1 = 0; $i1 < $Np[1]; $i1++)
    			$SQ += ($xj[$i1] - $xa) * ($xj[$i1] - $xa);
    		$SQ *= ($Np[0] * $N);
    
    		$SI = 0.0;
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			for ($i2 = 0; $i2 < $Np[1]; $i2++)
    				$SI += ($xij[$i1][$i2] - $xi[$i1] - $xj[$i2] + $xa) * ($xij[$i1][$i2] - $xi[$i1] - $xj[$i2] + $xa);
    		}
    		$SI *= $N;
    
    		$SE = 0.0;
    		for ($i1 = 0; $i1 < $Np[0]; $i1++) {
    			for ($i2 = 0; $i2 < $Np[1]; $i2++) {
    				for ($i3 = 0; $i3 < $N; $i3++)
    					$SE += ($x[$i1][$i2][$i3] - $xij[$i1][$i2]) * ($x[$i1][$i2][$i3] - $xij[$i1][$i2]);
    			}
    		}
    
    		$VP = $SP / ($Np[0] - 1);
    		$VQ = $SQ / ($Np[1] - 1);
    		$VI = $SI / (($Np[0] - 1) * ($Np[1] - 1));
    		$VE = $SE / ($Np[0] * $Np[1] * ($N - 1));
    
    		$FP = $VP / $VE;
    		$FQ = $VQ / $VE;
    		$FI = $VI / $VE;
    
    		$p    = 0.01 * $a;
    		$dof2 = $Np[0] * $Np[1] * ($N - 1);
    
    		printf("全変動: 平方和 %.2f 自由度 %d\n", $SP+$SQ+$SI+$SE, $Np[0]*$Np[1]*$N-1);
    		$dof1 = $Np[0] - 1;
    		printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", $name[0], $SP, $Np[0]-1, $VP, $FP, $a, p_f($sw));
    		$dof1 = $Np[1] - 1;
    		printf("%sの水準間: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", $name[1], $SQ, $Np[1]-1, $VQ, $FQ, $a, p_f($sw));
    		$dof1 = ($Np[0] - 1) * ($Np[1] - 1);
    		printf("相互作用: 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d%値 %.2f\n", $SI, ($Np[0]-1)*($Np[1]-1), $VI, $FI, $a, p_f($sw));
    		printf("水準内: 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", $SE, $Np[0]*$Np[1]*($N-1), $VE);
    	}
    }
    
    /*********************************************************/
    /* 二分法による非線形方程式(f(x)=0)の解                  */
    /*      f : f(x)を計算する関数名                         */
    /*      x1,x2 : 初期値                                   */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
    /*      max : 最大試行回数                               */
    /*      ind : 実際の試行回数                             */
    /*            (負の時は解を得ることができなかった)     */
    /*      return : 解                                      */
    /*********************************************************/
    function bisection($f, $x1, $x2, $eps1, $eps2, $max, &$ind)
    {
    	$x0 = 0.0;
    	$f1 = $f($x1);
    	$f2 = $f($x2);
    
    	if ($f1*$f2 > 0.0)
    		$ind = -1;
    
    	else {
    		$ind = 0;
    		if ($f1*$f2 == 0.0)
    			$x0 = ($f1 == 0.0) ? $x1 : $x2;
    		else {
    			$sw = 0;
    			while ($sw == 0 && $ind >= 0) {
    				$sw   = 1;
    				$ind += 1;
    				$x0   = 0.5 * ($x1 + $x2);
    				$f0   = $f($x0);
    
    				if (abs($f0) > $eps2) {
    					if ($ind <= $max) {
    						if (abs($x1-$x2) > $eps1 && abs($x1-$x2) > $eps1*abs($x2)) {
    							$sw = 0;
    							if ($f0*$f1 < 0.0) {
    								$x2 = $x0;
    								$f2 = $f0;
    							}
    							else {
    								$x1 = $x0;
    								$f1 = $f0;
    							}
    						}
    					}
    					else
    						$ind = -1;
    				}
    			}
    		}
    	}
    
    	return $x0;
    }
    
    /*************************************************/
    /* 標準正規分布N(0,1)の計算(P(X = x), P(X < x))*/
    /*      w : P(X = x)                             */
    /*      return : P(X < x)                        */
    /*************************************************/
    function normal($x, &$w)
    {
    	$pi = 4.0 * atan(1.0);
    /*
         確率密度関数(定義式)
    */
    	$w = exp(-0.5 * $x * $x) / sqrt(2.0 * $pi);
    /*
         確率分布関数(近似式を使用)
    */
    	$y = 0.70710678118654 * abs($x);
    	$z = 1.0 + $y * (0.0705230784 + $y * (0.0422820123 +
            $y * (0.0092705272 + $y * (0.0001520143 + $y * (0.0002765672 +
            $y * 0.0000430638)))));
    	$P = 1.0 - pow($z, -16.0);
    
    	if ($x < 0.0)
    		$P = 0.5 - 0.5 * $P;
    	else
    		$P = 0.5 + 0.5 * $P;
    
    	return $P;
    }
    
    /******************************************************************/
    /* 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用) */
    /*      ind : >= 0 : normal(収束回数)                           */
    /*            = -1 : 収束しなかった                               */
    /******************************************************************/
    function p_normal(&$ind)
    {
    	$u   = bisection("normal_f", -7.0, 7.0, 1.0e-6, 1.0e-10, 100, $sw);
    	$ind = $sw;
    
    	return $u;
    }
    
    /******************************/
    /* 1.0 - p - P(X>x)(関数値) */
    /******************************/
    function normal_f($x)
    {
    	global $p;
    	return 1.0 - $p - normal($x, $y);
    }
    
    /*****************************************************/
    /* Newton法による非線形方程式(f(x)=0)の解            */
    /*      f : f(x)を計算する関数名                     */
    /*      df : f(x)の微分を計算する関数名              */
    /*      x0 : 初期値                                  */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    /*      max : 最大試行回数                           */
    /*      ind : 実際の試行回数                         */
    /*            (負の時は解を得ることができなかった) */
    /*      return : 解                                  */
    /*****************************************************/
    function newton($f, $df, $x0, $eps1, $eps2, $max, &$ind)
    {
    	$x1  = $x0;
    	$x   = $x1;
    	$ind = 0;
    	$sw  = 0;
    
    	while ($sw == 0 && $ind >= 0) {
    
    		$sw   = 1;
    		$ind += 1;
    		$g    = $f($x1);
    
    		if (abs($g) > $eps2) {
    			if ($ind <= $max) {
    				$dg = $df($x1);
    				if (abs($dg) > $eps2) {
    					$x = $x1 - $g / $dg;
    					if (abs($x-$x1) > $eps1 && abs($x-$x1) > $eps1*abs($x)) {
    						$x1 = $x;
    						$sw = 0;
    					}
    				}
    				else
    					$ind = -1;
    			}
    			else
    				$ind = -1;
    		}
    	}
    
    	return $x;
    }
    
    /****************************************/
    /* Γ(x)の計算(ガンマ関数,近似式) */
    /*      ier : =0 : normal               */
    /*            =-1 : x=-n (n=0,1,2,・・・)  */
    /*      return : 結果                   */
    /****************************************/
    function gamma($x, &$ier)
    {
    	$ier = 0;
    
    	if ($x > 5.0) {
    		$v = 1.0 / $x;
    		$s = ((((((-0.000592166437354 * $v + 0.0000697281375837) * $v +
                0.00078403922172) * $v - 0.000229472093621) * $v -
                0.00268132716049) * $v + 0.00347222222222) * $v +
                0.0833333333333) * $v + 1.0;
    		$g = 2.506628274631001 * exp(-$x) * pow($x,$x-0.5) * $s;
    	}
    
    	else {
    
    		$err = 1.0e-20;
    		$w   = $x;
    		$t   = 1.0;
    
    		if ($x < 1.5) {
    
    			if ($x < $err) {
    				$k = intval($x);
    				$y = floatval($k) - $x;
    				if (abs($y) < $err || abs(1.0-$y) < $err)
    					$ier = -1;
    			}
    
    			if ($ier == 0) {
    				while ($w < 1.5) {
    					$t /= $w;
    					$w += 1.0;
    				}
    			}
    		}
    
    		else {
    			if ($w > 2.5) {
    				while ($w > 2.5) {
    					$w -= 1.0;
    					$t *= $w;
    				}
    			}
    		}
    
    		$w -= 2.0;
    		$g  = (((((((0.0021385778 * $w - 0.0034961289) * $w + 
                 0.0122995771) * $w - 0.00012513767) * $w + 0.0740648982) * $w +
                 0.0815652323) * $w + 0.411849671) * $w + 0.422784604) * $w +
                 0.999999926;
    		$g *= $t;
    	}
    
    	return $g;
    }
    
    /****************************************/
    /* f分布の計算(P(X = ff), P(X < ff)) */
    /*      dd : P(X = ff)                  */
    /*      df1,df2 : 自由度                */
    /*      return : P(X < ff)              */
    /****************************************/
    
    function f($ff, &$dd, $df1, $df2)
    {
    	$pi = 4.0 * atan(1.0);
    
    	if ($ff < 1.0e-10)
    		$ff = 1.0e-10;
    
    	$x  = $ff * $df1 / ($ff * $df1 + $df2);
    
    	if($df1%2 == 0) {
    		if ($df2%2 == 0) {
    			$u  = $x * (1.0 - $x);
    			$pp = $x;
    			$ia = 2;
    			$ib = 2;
    		}
    		else {
    			$u  = 0.5 * $x * sqrt(1.0-$x);
    			$pp = 1.0 - sqrt(1.0-$x);
    			$ia = 2;
    			$ib = 1;
    		}
    	}
    
    	else {
    		if ($df2%2 == 0) {
    			$u  = 0.5 * sqrt($x) * (1.0 - $x);
    			$pp = sqrt($x);
    			$ia = 1;
    			$ib = 2;
    		}
    		else {
    			$u  = sqrt($x*(1.0-$x)) / $pi;
    			$pp = 1.0 - 2.0 * atan2(sqrt(1.0-$x), sqrt($x)) / $pi;
    			$ia = 1;
    			$ib = 1;
    		}
    	}
    
    	if ($ia != $df1) {
    		for ($i1 = $ia; $i1 <= $df1-2; $i1 += 2) {
    			$pp -= 2.0 * $u / $i1;
    			$u  *= $x * ($i1 + $ib) / $i1;
    		}
    	}
    
    	if ($ib != $df2) {
    		for ($i1 = $ib; $i1 <= $df2-2; $i1 += 2) {
    			$pp += 2.0 * $u / $i1;
    			$u  *= (1.0 - $x) * ($i1 + $df1) / $i1;
    		}
    	}
    
    	$dd = $u / $ff;
    
    	return $pp;
    }
    
    /****************************************/
    /* f分布のp%値(P(X > u) = 0.01p)   */
    /*      ind : >= 0 : normal(収束回数) */
    /*            = -1 : 収束しなかった     */
    /****************************************/
    
    function p_f(&$ind)
    {
    	global $p, $dof1, $dof2;
    
    	$ff  = 0.0;
    	$sw  = 0;
    	$MAX = 340;
    /*
         初期値計算の準備
              while文は,大きな自由度によるガンマ関数の
              オーバーフローを避けるため
    */
    	while ($sw >= 0) {
    
    		$df1 = 0.5 * ($dof1 - $sw);
    		$df2 = 0.5 * $dof2;
    		$a   = 2.0 / (9.0 * ($dof1 - $sw));
    		$a1  = 1.0 - $a;
    		$b   = 2.0 / (9.0 * $dof2);
    		$b1  = 1.0 - $b;
    
    		$yq  = p_normal($ind);
    
    		$e   = $b1 * $b1 - $b * $yq * $yq;
    
    		if ($e > 0.8 || ($dof1+$dof2-$sw) <= $MAX)
    			$sw = -1;
    		else {
    			$sw += 1;
    			if (($dof1-$sw) == 0)
    				$sw = -2;
    		}
    	}
    
    	if ($sw == -2)
    		$ind = -1;
    
    	else {
    /*
         f0 : 初期値
    */
    		if ($e > 0.8) {
    			$x  = ($a1 * $b1 + $yq * sqrt($a1*$a1*$b+$a*$e)) / $e;
    			$f0 = pow($x, 3.0);
    		}
    		else {
    			$y1 = pow(floatval($dof2), $df2-1.0);
    			$y2 = pow(floatval($dof1), $df2);
    			$x  = gamma($df1+$df2, $ind) / gamma($df1, $ind) / gamma($df2, $ind) *
                     2.0 * $y1 / $y2 / $p;
    			$f0 = pow($x, 2.0/$dof2);
    		}
    /*
         ニュートン法
    */
    		$ff = newton("f_f", "f_df", $f0, 1.0e-6, 1.0e-10, 100, $ind);
    	}
    
    	return $ff;
    }
    
    /********************************/
    /* 1.0 - p - P(X > x)(関数値) */
    /********************************/
    function f_f($x)
    {
    	global $p, $dof1, $dof2;
    	return f($x, $y, $dof1, $dof2) - 1.0 + $p;
    }
    
    /**************************/
    /* P(X = x)(関数の微分) */
    /**************************/
    function f_df($x)
    {
    	global $dof1, $dof2;
    
    	$z = f($x, $y, $dof1, $dof2);
    
    	return $y;
    }
    
    /*
    -------- 一元配置法に対するデータ例(コメント部分を除いて下さい)--------
    1 6 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    工場 3   // 因子の名前 水準の数
    3.1 4.7 5.1   // 各水準に対する1番目のデータ
    4.1 5.6 3.7   // 各水準に対する2番目のデータ
    3.3 4.3 4.5   // 各水準に対する3番目のデータ
    3.9 5.9 6.0   // 各水準に対する4番目のデータ
    3.7 6.1 3.9   // 各水準に対する5番目のデータ
    2.4 4.2 5.4   // 各水準に対する6番目のデータ
    
    -------- 二元配置法に対するデータ例(コメント部分を除いて下さい)--------
    2 3 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    薬剤 5   // 1番目の因子の名前 その水準の数
    品種 2   // 2番目の因子の名前 その水準の数
    3 4 12 -4 -4   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する1番目のデータ
    8 -8 31 12 19   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する2番目のデータ
    7 -5 8 0 23   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する3番目のデータ
    8 -10 9 10 15   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する1番目のデータ
    -5 11 26 -1 13   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する2番目のデータ
    10 -6 13 -7 -6   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する3番目のデータ
    */
    
    ?>
    			

  5. Ruby

    ############################
    # 分散分析
    #      coded by Y.Suganuma
    ############################
    
    ############################################
    # Γ(x)の計算(ガンマ関数,近似式)
    #      ier : =0 : normal
    #            =-1 : x=-n (n=0,1,2,・・・)
    #      return : 結果
    #      coded by Y.Suganuma
    ############################################
    
    def gamma(x, ier)
    
    	ier[0] = 0
    
    	if x > 5.0
    		v = 1.0 / x
    		s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v + 0.00078403922172) * v - 0.000229472093621) * v - 0.00268132716049) * v + 0.00347222222222) * v + 0.0833333333333) * v + 1.0
    		g = 2.506628274631001 * Math.exp(-x) * (x ** (x-0.5)) * s
    
    	else
    
    		err = 1.0e-20
    		w   = x
    		t   = 1.0
    
    		if x < 1.5
    
    			if x < err
    				k = Integer(x)
    				y = Float(k) - x
    				if y.abs() < err or (1.0-y).abs() < err
    					ier[0] = -1
    				end
    			end
    
    			if ier[0] == 0
    				while w < 1.5
    					t /= w
    					w += 1.0
    				end
    			end
    
    		else
    			if w > 2.5
    				while w > 2.5
    					w -= 1.0
    					t *= w
    				end
    			end
    		end
    
    		w -= 2.0
    		g  = (((((((0.0021385778 * w - 0.0034961289) * w + 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w + 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w + 0.999999926
    		g *= t
    	end
    
    	return g
    end
    
    ################################################
    # 標準正規分布N(0,1)の計算(P(X = x), P(X < x))
    #      w : P(X = x)
    #      return : P(X < x)
    ################################################
    
    def normal(x, w)
    			# 確率密度関数(定義式)
    	w[0] = Math.exp(-0.5 * x * x) / Math.sqrt(2.0*Math::PI)
    			# 確率分布関数(近似式を使用)
    	y = 0.70710678118654 * x.abs()
    	z = 1.0 + y * (0.0705230784 + y * (0.0422820123 + y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 + y * 0.0000430638)))))
    	pp = 1.0 - z ** (-16.0)
    
    	if x < 0.0
    		pp = 0.5 - 0.5 * pp
    	else
    		pp = 0.5 + 0.5 * pp
    	end
    
    	return pp
    end
    
    #######################################
    # f分布の計算(P(X = ff), P(X < ff))
    #      dd : P(X = ff)
    #      df1,df2 : 自由度
    #      return : P(X < ff)
    #######################################
    
    def f(ff, dd, df1, df2)
    
    	if ff < 1.0e-10
    		ff = 1.0e-10
    	end
    
    	x = ff * df1 / (ff * df1 + df2)
    
    	if df1%2 == 0
    		if df2%2 == 0
    			u  = x * (1.0 - x)
    			pp = x
    			ia = 2
    			ib = 2
    		else
    			u  = 0.5 * x * Math.sqrt(1.0-x)
    			pp = 1.0 - Math.sqrt(1.0-x)
    			ia = 2
    			ib = 1
    		end
    
    	else
    		if df2%2 == 0
    			u  = 0.5 * Math.sqrt(x) * (1.0 - x)
    			pp = Math.sqrt(x)
    			ia = 1
    			ib = 2
    		else
    			u  = Math.sqrt(x*(1.0-x)) / Math::PI
    			pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / Math::PI
    			ia = 1
    			ib = 1
    		end
    	end
    
    	if ia != df1
    		i1 = ia
    		while i1 < df1-1
    			pp -= 2.0 * u / i1
    			u  *= x * (i1 + ib) / i1
    			i1 += 2
    		end
    	end
    
    	if ib != df2
    		i1 = ib
    		while i1 < df2-1
    			pp += 2.0 * u / i1
    			u  *= (1.0 - x) * (i1 + df1) / i1
    			i1 += 2
    		end
    	end
    
    	dd[0] = u / ff
    
    	return pp
    end
    
    ############################################
    # 二分法による非線形方程式(f(x)=0)の解
    #      x1,x2 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      ind : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      fn : f(x)を計算する関数名
    #      return : 解
    #      coded by Y.Suganuma
    ############################################
    
    def bisection(x1, x2, eps1, eps2, max, ind, &fn)
    
    	x0 = 0.0
    	f1 = fn.call(x1)
    	f2 = fn.call(x2)
    
    	if f1*f2 > 0.0
    		ind[0] = -1
    
    	else
    		ind[0] = 0
    		if f1*f2 == 0.0
    			if f1 == 0.0
    				x0 = x1
    			else
    				x0 = x2
    			end
    		else
    			sw = 0
    			while sw == 0 && ind[0] >= 0
    				sw      = 1
    				ind[0] += 1
    				x0     = 0.5 * (x1 + x2)
    				f0     = fn.call(x0)
    
    				if f0.abs() > eps2
    					if ind[0] <= max
    						if (x1-x2).abs() > eps1 && (x1-x2).abs() > eps1*x2.abs()
    							sw = 0
    							if f0*f1 < 0.0
    								x2 = x0
    								f2 = f0
    							else
    								x1 = x0
    								f1 = f0
    							end
    						end
    					else
    						ind[0] = -1
    					end
    				end
    			end
    		end
    	end
    
    	return x0
    end
    
    ############################################
    # Newton法による非線形方程式(f(x)=0)の解
    #      x0 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      ind : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      fn : f(x)とその微分を計算する関数名
    #      return : 解
    #      coded by Y.Suganuma
    ############################################
    def newton(x0, eps1, eps2, max, ind, &fn) 
    
    	x1     = x0
    	x      = x1
    	ind[0] = 0
    	sw     = 0
    
    	while sw == 0 and ind[0] >= 0 
    
    		sw      = 1
    		ind[0] += 1
    		g       = fn.call(0, x1)
    
    		if g.abs() > eps2 
    			if ind[0] <= max 
    				dg = fn.call(1, x1)
    				if dg.abs() > eps2 
    					x = x1 - g / dg
    					if (x-x1).abs() > eps1 && (x-x1).abs() > eps1*x.abs() 
    						x1 = x
    						sw = 0
    					end
    				else 
    					ind[0] = -1
    				end
    			else 
    				ind[0] = -1
    			end
    		end
    	end
    
    	return x
    end
    
    ################################################################
    # 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用)
    #      ind : >= 0 : normal(収束回数)
    #            = -1 : 収束しなかった
    ################################################################
    def p_normal(ind)
    
    	normal_snx = Proc.new { |x|
    		y = Array.new(1)
    		1.0 - $p - normal(x, y)
    	}
    
    	u = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, ind, &normal_snx)
    	return u
    end
    
    ######################################
    # f分布のp%値(P(X > u) = 0.01p)
    #      ind : >= 0 : normal(収束回数)
    #            = -1 : 収束しなかった
    ######################################
    
    def p_f(ind)
    
    	f_snx = Proc.new { |sw, x|
    		y = Array.new(1)
    		z = f(x, y, $dof1, $dof2)
    		if sw == 0
    			z - 1.0 + $p
    		else
    			z = y[0]
    		end
    	}
    
    	max = 340
    	ff  = 0.0
    	sw  = 0
    			# 初期値計算の準備
    			# while文は,大きな自由度によるガンマ関数の
    			# オーバーフローを避けるため
    	while sw >= 0
    
    		df1 = 0.5 * ($dof1 - sw)
    		df2 = 0.5 * $dof2
    		a   = 2.0 / (9.0 * ($dof1 - sw))
    		a1  = 1.0 - a
    		b   = 2.0 / (9.0 * $dof2)
    		b1  = 1.0 - b
    
    		yq  = p_normal(ind)
    
    		e   = b1 * b1 - b * yq * yq
    
    		if e > 0.8 or $dof1+$dof2-sw <= max
    			sw = -1
    		else
    			sw += 1
    			if ($dof1-sw) == 0
    				sw = -2
    			end
    		end
    	end
    
    	if sw == -2
    		ind[0] = -1
    
    	else
    			# f0 : 初期値
    		if e > 0.8
    			x  = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e
    			f0 = x ** 3.0
    		else
    			y1 = Float($dof2) ** (df2-1.0)
    			y2 = Float($dof1) ** df2
    			x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) * 2.0 * y1 / y2 / $p
    			f0 = x ** (2.0/$dof2)
    		end
    			# ニュートン法
    		ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind, &f_snx)
    	end
    
    	return ff
    end
    
    #######################################
    # 分散分析(一元配置法と二元配置法)
    #      method : 1 or 2
    #      np[0] : 因子1の水準数
    #        [1] : 因子2の水準数
    #      nn : データ数
    #      name[0] : 因子1の名前
    #          [1] : 因子2の名前
    #      a : a %値
    #      x : データ
    #           coded by Y.Suganuma
    #######################################
    def aov(method, np, nn, name, a, x)
    
    			# 一元配置法
    
    	if method == 1
    
    		xi = Array.new(np[0])
    		for i1 in 0 ... np[0]
    			xi[i1] = 0.0
    			for i2 in 0 ... nn
    				xi[i1] += x[i1][0][i2]
    			end
    			xi[i1] /= nn
    		end
    
    		xa = 0.0
    		for i1 in 0 ... np[0]
    			for i2 in 0 ... nn
    				xa += x[i1][0][i2]
    			end
    		end
    		xa /= (np[0] * nn)
    
    		sp = 0.0
    		for i1 in 0 ... np[0]
    			sp += (xi[i1] - xa) * (xi[i1] - xa)
    		end
    		sp *= nn
    
    		se = 0.0
    		for i1 in 0 ... np[0]
    			for i2 in 0 ... nn
    				se += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1])
    			end
    		end
    
    		vp = sp / (np[0] - 1)
    		ve = se / (np[0] * (nn - 1))
    
    		fp = vp / ve
    
    		$p    = 0.01 * a
    		$dof2 = np[0] * (nn - 1)
    
    		sw = Array.new(1)
    		print("全変動: 平方和 " + String(sp+se) + " 自由度 " + String(np[0]*nn-1) + "\n")
    		$dof1 = np[0] - 1
    		print(name[0])
    		print("の水準間: 平方和 " + String(sp) + " 自由度 " + String(np[0]-1) + " 不偏分散 " + String(vp) + " F " + String(fp) + " " + String(a) + "%値 " + String(p_f(sw)) + "\n")
    		print("水準内: 平方和 " + String(se) + " 自由度 " + String(np[0]*(nn-1)) + " 不偏分散 " + String(ve) + "\n")
    
    			# 二元配置法
    
    	else
    
    		xi = Array.new(np[0])
    		for i1 in 0 ... np[0]
    			xi[i1] = 0.0
    			for i2 in 0 ... np[1]
    				for i3 in 0 ... nn
    					xi[i1] += x[i1][i2][i3]
    				end
    			end
    			xi[i1] /= (np[1] * nn)
    		end
    
    		xj = Array.new(np[1])
    		for i1 in 0 ... np[1]
    			xj[i1] = 0.0
    			for i2 in 0 ... np[0]
    				for i3 in 0 ... nn
    					xj[i1] += x[i2][i1][i3]
    				end
    			end
    			xj[i1] /= (np[0] * nn)
    		end
    
    		xij = Array.new(np[0])
    		for i1 in 0 ... np[0]
    			xij[i1] = Array.new(np[1])
    			for i2 in 0 ... np[1]
    				xij[i1][i2] = 0.0
    				for i3 in 0 ... nn
    					xij[i1][i2] += x[i1][i2][i3]
    				end
    				xij[i1][i2] /= nn
    			end
    		end
    
    		xa = 0.0
    		for i1 in 0 ... np[0]
    			for i2 in 0 ... np[1]
    				for i3 in 0 ... nn
    					xa += x[i1][i2][i3]
    				end
    			end
    		end
    		xa /= (np[0] * np[1] * nn)
    
    		sp = 0.0
    		for i1 in 0 ... np[0]
    			sp += (xi[i1] - xa) * (xi[i1] - xa)
    		end
    		sp *= (np[1] * nn)
    
    		sq = 0.0
    		for i1 in 0 ... np[1]
    			sq += (xj[i1] - xa) * (xj[i1] - xa)
    		end
    		sq *= (np[0] * nn)
    
    		si = 0.0
    		for i1 in 0 ... np[0]
    			for i2 in 0 ... np[1]
    				si += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa)
    			end
    		end
    		si *= nn
    
    		se = 0.0
    		for i1 in 0 ... np[0]
    			for i2 in 0 ... np[1]
    				for i3 in 0 ... nn
    					se += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2])
    				end
    			end
    		end
    
    		vp = sp / (np[0] - 1)
    		vq = sq / (np[1] - 1)
    		vi = si / ((np[0] - 1) * (np[1] - 1))
    		ve = se / (np[0] * np[1] * (nn - 1))
    
    		fp = vp / ve
    		fq = vq / ve
    		fi = vi / ve
    
    		$p    = 0.01 * a
    		$dof2 = np[0] * np[1] * (nn - 1)
    
    		sw = Array.new(1)
    		print("全変動: 平方和 " + String(sp+sq+si+se) + " 自由度 " + String(np[0]*np[1]*nn-1) + "\n")
    		$dof1 = np[0] - 1
    		print(name[0])
    		print("の水準間: 平方和 " + String(sp) + " 自由度 " + String(np[0]-1) + " 不偏分散 " + String(vp) + " F " + String(fp) + " " + String(a) + "%値 " + String(p_f(sw)) + "\n")
    		$dof1 = np[1] - 1
    		print(name[1])
    		print("の水準間: 平方和 " + String(sq) + " 自由度 " + String(np[1]-1) + " 不偏分散 " + String(vq) + " F " + String(fq) + " " + String(a) + "%値 " + String(p_f(sw)) + "\n")
    		$dof1 = (np[0] - 1) * (np[1] - 1)
    		print("相互作用: 平方和 " + String(si) + " 自由度 " + String((np[0]-1)*(np[1]-1)) + " 不偏分散 " + String(vi) + " F " + String(fi) + " " + String(a) + "%値 " + String(p_f(sw)) + "\n")
    		print("水準内: 平方和 " + String(se) + " 自由度 " + String(np[0]*np[1]*(nn-1)) + " 不偏分散 " + String(ve) + "\n")
    	end
    end
    
    $p    = 0.0
    $dof1 = 1
    $dof2 = 2
    
    name   = ["", ""]
    np     = [1, 1]
    ss     = gets().split(" ")
    method = Integer(ss[0])   # 因子の数
    nn     = Integer(ss[1])   # 各水準におけるデータ数
    a      = Float(ss[2])   # 有意水準(%)
    
    if method == 1 or method == 2
    	for i1 in 0 ... method
    		ss       = gets().split(" ")
    		name[i1] = ss[0]
    		np[i1]   = Integer(ss[1])
    	end
    	x = Array.new(np[0])
    	for i1 in 0 ... np[0]
    		x[i1] = Array.new(np[1])
    		for i2 in 0 ... np[1]
    			x[i1][i2] = Array.new(nn)
    		end
    	end
    	for i1 in 0 ... np[1]
    		for i2 in 0 ... nn
    			ss = gets().split(" ")
    			for i3 in 0 ... np[0]
    				x[i3][i1][i2] = Float(ss[i3])
    			end
    		end
    	end
    	aov(method, np, nn, name, a, x)
    else
    	print("一元配置法,または,二元配置法だけです!\n")
    end
    
    =begin
    -------- 一元配置法に対するデータ例(コメント部分を除いて下さい)--------
    1 6 5   # 因子の数 各水準におけるデータ数 有意水準(%)
    工場 3   # 因子の名前 水準の数
    3.1 4.7 5.1   # 各水準に対する1番目のデータ
    4.1 5.6 3.7   # 各水準に対する2番目のデータ
    3.3 4.3 4.5   # 各水準に対する3番目のデータ
    3.9 5.9 6.0   # 各水準に対する4番目のデータ
    3.7 6.1 3.9   # 各水準に対する5番目のデータ
    2.4 4.2 5.4   # 各水準に対する6番目のデータ
    
    -------- 二元配置法に対するデータ例(コメント部分を除いて下さい)--------
    2 3 5   # 因子の数 各水準におけるデータ数 有意水準(%)
    薬剤 5   # 1番目の因子の名前 その水準の数
    品種 2   # 2番目の因子の名前 その水準の数
    3 4 12 -4 -4   # 1番目の因子の各水準に対して,2番目の因子の水準1に対する1番目のデータ
    8 -8 31 12 19   # 1番目の因子の各水準に対して,2番目の因子の水準1に対する2番目のデータ
    7 -5 8 0 23   # 1番目の因子の各水準に対して,2番目の因子の水準1に対する3番目のデータ
    8 -10 9 10 15   # 1番目の因子の各水準に対して,2番目の因子の水準2に対する1番目のデータ
    -5 11 26 -1 13   # 1番目の因子の各水準に対して,2番目の因子の水準2に対する2番目のデータ
    10 -6 13 -7 -6   # 1番目の因子の各水準に対して,2番目の因子の水準2に対する3番目のデータ
    =end
    			

  6. Python

    # -*- coding: UTF-8 -*-
    import numpy as np
    import sys
    from math import *
    
    ############################################
    # Γ(x)の計算(ガンマ関数,近似式)
    #      ier : =0 : normal
    #            =-1 : x=-n (n=0,1,2,・・・)
    #      return : 結果
    #      coded by Y.Suganuma
    ############################################
    
    def gamma(x, ier) :
    
    	ier[0] = 0
    
    	if x > 5.0 :
    		v = 1.0 / x
    		s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v + 0.00078403922172) * v - 0.000229472093621) * v - 0.00268132716049) * v + 0.00347222222222) * v + 0.0833333333333) * v + 1.0
    		g = 2.506628274631001 * exp(-x) * pow(x,x-0.5) * s
    
    	else :
    
    		err = 1.0e-20
    		w   = x
    		t   = 1.0
    
    		if x < 1.5 :
    
    			if x < err :
    				k = int(x)
    				y = float(k) - x
    				if abs(y) < err or abs(1.0-y) < err :
    					ier[0] = -1
    
    			if ier[0] == 0 :
    				while w < 1.5 :
    					t /= w
    					w += 1.0
    
    		else :
    			if w > 2.5 :
    				while w > 2.5 :
    					w -= 1.0
    					t *= w
    
    		w -= 2.0
    		g  = (((((((0.0021385778 * w - 0.0034961289) * w + 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w + 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w + 0.999999926
    		g *= t
    
    	return g
    
    ################################################
    # 標準正規分布N(0,1)の計算(P(X = x), P(X < x))
    #      w : P(X = x)
    #      return : P(X < x)
    ################################################
    
    def normal(x, w) :
    			# 確率密度関数(定義式)
    	w[0] = exp(-0.5 * x * x) / sqrt(2.0*pi)
    			# 確率分布関数(近似式を使用)
    	y = 0.70710678118654 * abs(x)
    	z = 1.0 + y * (0.0705230784 + y * (0.0422820123 + y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 + y * 0.0000430638)))))
    	P = 1.0 - z ** (-16.0)
    
    	if x < 0.0 :
    		P = 0.5 - 0.5 * P
    	else :
    		P = 0.5 + 0.5 * P
    
    	return P
    
    #######################################
    # f分布の計算(P(X = ff), P(X < ff))
    #      dd : P(X = ff)
    #      df1,df2 : 自由度
    #      return : P(X < ff)
    #######################################
    
    def f(ff, dd, df1, df2) :
    
    	if ff < 1.0e-10 :
    		ff = 1.0e-10
    
    	x = ff * df1 / (ff * df1 + df2)
    
    	if df1%2 == 0 :
    		if df2%2 == 0 :
    			u  = x * (1.0 - x)
    			pp = x
    			ia = 2
    			ib = 2
    		else :
    			u  = 0.5 * x * sqrt(1.0-x)
    			pp = 1.0 - sqrt(1.0-x)
    			ia = 2
    			ib = 1
    
    	else :
    		if df2%2 == 0 :
    			u  = 0.5 * sqrt(x) * (1.0 - x)
    			pp = sqrt(x)
    			ia = 1
    			ib = 2
    		else :
    			u  = sqrt(x*(1.0-x)) / pi
    			pp = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi
    			ia = 1
    			ib = 1
    
    	if ia != df1 :
    		for i1 in range(ia, df1-1, 2) :
    			pp -= 2.0 * u / i1
    			u  *= x * (i1 + ib) / i1
    
    	if ib != df2 :
    		for i1 in range(ib, df2-1, 2) :
    			pp += 2.0 * u / i1
    			u  *= (1.0 - x) * (i1 + df1) / i1
    
    	dd[0] = u / ff
    
    	return pp
    
    ############################################
    # 二分法による非線形方程式(f(x)=0)の解
    #      fn : f(x)を計算する関数名
    #      x1,x2 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      ind : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      return : 解
    #      coded by Y.Suganuma
    ############################################
    
    def bisection(fn, x1, x2, eps1, eps2, max, ind) :
    
    	x0 = 0.0
    	f1 = fn(x1)
    	f2 = fn(x2)
    
    	if f1*f2 > 0.0 :
    		ind[0] = -1
    
    	else :
    		ind[0] = 0
    		if f1*f2 == 0.0 :
    			if f1 == 0.0 :
    				x0 = x1
    			else :
    				x0 = x2
    		else :
    			sw = 0
    			while sw == 0 and ind[0] >= 0 :
    				sw      = 1
    				ind[0] += 1
    				x0     = 0.5 * (x1 + x2)
    				f0     = fn(x0)
    
    				if abs(f0) > eps2 :
    					if ind[0] <= max :
    						if abs(x1-x2) > eps1 and abs(x1-x2) > eps1*abs(x2) :
    							sw = 0
    							if f0*f1 < 0.0 :
    								x2 = x0
    								f2 = f0
    							else :
    								x1 = x0
    								f1 = f0
    					else :
    						ind[0] = -1
    
    	return x0
    
    ############################################
    # Newton法による非線形方程式(f(x)=0)の解
    #      fn : f(x)を計算する関数名
    #      dfn : f(x)の微分を計算する関数名
    #      x0 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      ind : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      return : 解
    #      coded by Y.Suganuma
    ############################################
    
    def newton(fn, dfn, x0, eps1, eps2, max, ind) :
    
    	x1     = x0
    	x      = x1
    	ind[0] = 0
    	sw     = 0
    
    	while sw == 0 and ind[0] >= 0 :
    
    		sw      = 1
    		ind[0] += 1
    		g       = fn(x1)
    
    		if abs(g) > eps2 :
    			if ind[0] <= max :
    				dg = dfn(x1)
    				if abs(dg) > eps2 :
    					x = x1 - g / dg
    					if abs(x-x1) > eps1 and abs(x-x1) > eps1*abs(x) :
    						x1 = x
    						sw = 0
    				else :
    					ind[0] = -1
    			else :
    				ind[0] = -1
    
    	return x
    
    ##########################################
    # 1.0 - p - P(X>x)(関数値, 標準正規分布)
    ##########################################
    def normal_f(x) :
    	y = np.empty(1, np.float)
    	return 1.0 - p - normal(x, y)
    
    ################################################################
    # 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用)
    #      ind : >= 0 : normal(収束回数)
    #            = -1 : 収束しなかった
    ################################################################
    def p_normal(ind) :
    	u = bisection(normal_f, -7.0, 7.0, 1.0e-6, 1.0e-10, 100, ind)
    	return u
    
    #######################################
    # 1.0 - p - P(X > x)(関数値, f 分布)
    #######################################
    def f_f(x) :
    	y = np.empty(1, np.float)
    	return f(x, y, dof1, dof2) - 1.0 + p
    
    #################################
    # P(X = x)(関数の微分, f 分布)
    #################################
    def f_df(x) :
    	y = np.empty(1, np.float)
    	z = f(x, y, dof1, dof2)
    	return y[0]
    
    ######################################
    # f分布のp%値(P(X > u) = 0.01p)
    #      ind : >= 0 : normal(収束回数)
    #            = -1 : 収束しなかった
    ######################################
    
    def p_f(ind) :
    
    	MAX = 340
    	ff  = 0.0
    	sw  = 0
    			# 初期値計算の準備
    			# while文は,大きな自由度によるガンマ関数の
    			# オーバーフローを避けるため
    	while sw >= 0 :
    
    		df1 = 0.5 * (dof1 - sw)
    		df2 = 0.5 * dof2
    		a   = 2.0 / (9.0 * (dof1 - sw))
    		a1  = 1.0 - a
    		b   = 2.0 / (9.0 * dof2)
    		b1  = 1.0 - b
    
    		yq  = p_normal(ind)
    
    		e   = b1 * b1 - b * yq * yq
    
    		if e > 0.8 or dof1+dof2-sw <= MAX :
    			sw = -1
    		else :
    			sw += 1
    			if (dof1-sw) == 0 :
    				sw = -2
    
    	if sw == -2 :
    		ind[0] = -1
    
    	else :
    			# f0 : 初期値
    		if e > 0.8 :
    			x  = (a1 * b1 + yq * sqrt(a1*a1*b+a*e)) / e
    			f0 = x ** 3.0
    		else :
    			y1 = float(dof2) ** (df2-1.0)
    			y2 = float(dof1) ** df2
    			x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) * 2.0 * y1 / y2 / p
    			f0 = x ** (2.0/dof2)
    			# ニュートン法
    		ff = newton(f_f, f_df, f0, 1.0e-6, 1.0e-10, 100, ind)
    
    	return ff
    
    #######################################
    # 分散分析(一元配置法と二元配置法)
    #      method : 1 or 2
    #      Np[0] : 因子1の水準数
    #        [1] : 因子2の水準数
    #      N : データ数
    #      name[0] : 因子1の名前
    #          [1] : 因子2の名前
    #      a : a %値
    #      x : データ
    #           coded by Y.Suganuma
    #######################################
    def aov(method, Np, N, name, a, x) :
    
    	global p, dof1, dof2
    
    			# 一元配置法
    
    	if method == 1 :
    
    		xi = np.empty(Np[0], np.float)
    		for i1 in range(0, Np[0]) :
    			xi[i1] = 0.0
    			for i2 in range(0, N) :
    				xi[i1] += x[i1][0][i2]
    			xi[i1] /= N
    
    		xa = 0.0
    		for i1 in range(0, Np[0]) :
    			for i2 in range(0, N) :
    				xa += x[i1][0][i2]
    		xa /= (Np[0] * N)
    
    		SP = 0.0
    		for i1 in range(0, Np[0]) :
    			SP += (xi[i1] - xa) * (xi[i1] - xa)
    		SP *= N
    
    		SE = 0.0
    		for i1 in range(0, Np[0]) :
    			for i2 in range(0, N) :
    				SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1])
    
    		VP = SP / (Np[0] - 1)
    		VE = SE / (Np[0] * (N - 1))
    
    		FP = VP / VE
    
    		p    = 0.01 * a
    		dof2 = Np[0] * (N - 1)
    
    		sw = np.empty(1, np.int)
    		print("全変動: 平方和 " + str(SP+SE) + " 自由度 " + str(Np[0]*N-1))
    		dof1 = Np[0] - 1
    		print(name[0] + "の水準間: 平方和 " + str(SP) + " 自由度 " + str(Np[0]-1) + " 不偏分散 " + str(VP) + " F " + str(FP) + " " + str(a) + "%値 " + str(p_f(sw)))
    		print("水準内: 平方和 " + str(SE) + " 自由度 " + str(Np[0]*(N-1)) + " 不偏分散 " + str(VE))
    
    			# 二元配置法
    
    	else :
    
    		xi = np.empty(Np[0], np.float)
    		for i1 in range(0, Np[0]) :
    			xi[i1] = 0.0
    			for i2 in range(0, Np[1]) :
    				for i3 in range(0, N) :
    					xi[i1] += x[i1][i2][i3]
    			xi[i1] /= (Np[1] * N)
    
    		xj = np.empty(Np[1], np.float)
    		for i1 in range(0, Np[1]) :
    			xj[i1] = 0.0
    			for i2 in range(0, Np[0]) :
    				for i3 in range(0, N) :
    					xj[i1] += x[i2][i1][i3]
    			xj[i1] /= (Np[0] * N)
    
    		xij = np.empty((Np[0], Np[1]), np.float)
    		for i1 in range(0, Np[0]) :
    			for i2 in range(0, Np[1]) :
    				xij[i1][i2] = 0.0
    				for i3 in range(0, N) :
    					xij[i1][i2] += x[i1][i2][i3]
    				xij[i1][i2] /= N
    
    		xa = 0.0
    		for i1 in range(0, Np[0]) :
    			for i2 in range(0, Np[1]) :
    				for i3 in range(0, N) :
    					xa += x[i1][i2][i3]
    		xa /= (Np[0] * Np[1] * N)
    
    		SP = 0.0
    		for i1 in range(0, Np[0]) :
    			SP += (xi[i1] - xa) * (xi[i1] - xa)
    		SP *= (Np[1] * N)
    
    		SQ = 0.0
    		for i1 in range(0, Np[1]) :
    			SQ += (xj[i1] - xa) * (xj[i1] - xa)
    		SQ *= (Np[0] * N)
    
    		SI = 0.0
    		for i1 in range(0, Np[0]) :
    			for i2 in range(0, Np[1]) :
    				SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa)
    		SI *= N
    
    		SE = 0.0
    		for i1 in range(0, Np[0]) :
    			for i2 in range(0, Np[1]) :
    				for i3 in range(0, N) :
    					SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2])
    
    		VP = SP / (Np[0] - 1)
    		VQ = SQ / (Np[1] - 1)
    		VI = SI / ((Np[0] - 1) * (Np[1] - 1))
    		VE = SE / (Np[0] * Np[1] * (N - 1))
    
    		FP = VP / VE
    		FQ = VQ / VE
    		FI = VI / VE
    
    		p    = 0.01 * a
    		dof2 = Np[0] * Np[1] * (N - 1)
    
    		sw = np.empty(1, np.int)
    		print("全変動: 平方和 " + str(SP+SQ+SI+SE) + " 自由度 " + str(Np[0]*Np[1]*N-1))
    		dof1 = Np[0] - 1
    		print(name[0] + "の水準間: 平方和 " + str(SP) + " 自由度 " + str(Np[0]-1) + " 不偏分散 " + str(VP) + " F " + str(FP) + " " + str(a) + "%値 " + str(p_f(sw)))
    		dof1 = Np[1] - 1
    		print(name[1] + "の水準間: 平方和 " + str(SQ) + " 自由度 " + str(Np[1]-1) + " 不偏分散 " + str(VQ) + " F " + str(FQ) + " " + str(a) + "%値 " + str(p_f(sw)))
    		dof1 = (Np[0] - 1) * (Np[1] - 1)
    		print("相互作用: 平方和 " + str(SI) + " 自由度 " + str((Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + str(VI) + " F " + str(FI) + " " + str(a) + "%値 " + str(p_f(sw)))
    		print("水準内: 平方和 " + str(SE) + " 自由度 " + str(Np[0]*Np[1]*(N-1)) + " 不偏分散 " + str(VE))
    
    p    = 0.0
    dof1 = 1
    dof2 = 2
    
    ############################
    # 分散分析
    #      coded by Y.Suganuma
    ############################
    
    name   = ["", ""]
    Np     = np.array([1, 1], np.int)
    line   = sys.stdin.readline()
    ss     = line.split()
    method = int(ss[0])   # 因子の数
    N      = int(ss[1])   # 各水準におけるデータ数
    a      = float(ss[2])   # 有意水準(%)
    
    if method == 1 or method == 2 :
    	for i1 in range(0, method) :
    		line   = sys.stdin.readline()
    		ss     = line.split()
    		name[i1] = ss[0]
    		Np[i1]   = int(ss[1])
    	x = np.empty((Np[0], Np[1], N), np.float)
    	for i1 in range(0, Np[1]) :
    		for i2 in range(0, N) :
    			line = sys.stdin.readline()
    			ss   = line.split()
    			for i3 in range(0, Np[0]) :
    				x[i3][i1][i2] = float(ss[i3])
    	aov(method, Np, N, name, a, x)
    else :
    	print("一元配置法,または,二元配置法だけです!")
    
    """
    -------- 一元配置法に対するデータ例(コメント部分を除いて下さい)--------
    1 6 5   # 因子の数 各水準におけるデータ数 有意水準(%)
    工場 3   # 因子の名前 水準の数
    3.1 4.7 5.1   # 各水準に対する1番目のデータ
    4.1 5.6 3.7   # 各水準に対する2番目のデータ
    3.3 4.3 4.5   # 各水準に対する3番目のデータ
    3.9 5.9 6.0   # 各水準に対する4番目のデータ
    3.7 6.1 3.9   # 各水準に対する5番目のデータ
    2.4 4.2 5.4   # 各水準に対する6番目のデータ
    
    -------- 二元配置法に対するデータ例(コメント部分を除いて下さい)--------
    2 3 5   # 因子の数 各水準におけるデータ数 有意水準(%)
    薬剤 5   # 1番目の因子の名前 その水準の数
    品種 2   # 2番目の因子の名前 その水準の数
    3 4 12 -4 -4   # 1番目の因子の各水準に対して,2番目の因子の水準1に対する1番目のデータ
    8 -8 31 12 19   # 1番目の因子の各水準に対して,2番目の因子の水準1に対する2番目のデータ
    7 -5 8 0 23   # 1番目の因子の各水準に対して,2番目の因子の水準1に対する3番目のデータ
    8 -10 9 10 15   # 1番目の因子の各水準に対して,2番目の因子の水準2に対する1番目のデータ
    -5 11 26 -1 13   # 1番目の因子の各水準に対して,2番目の因子の水準2に対する2番目のデータ
    10 -6 13 -7 -6   # 1番目の因子の各水準に対して,2番目の因子の水準2に対する3番目のデータ
    """
    			

  7. C#

    /****************************/
    /* 分散分析                 */
    /*      coded by Y.Suganuma */
    /****************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	double p;     // α%値を計算するとき時α/100を設定
    	int dof1, dof2;      // 自由度
    
    	public Test1()
    	{
    		int[] Np       = {1, 1};
    		char[] charSep = new char[] {' '};
    		String[] str   = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
    		int method     = int.Parse(str[0]);
    		int N          = int.Parse(str[1]);
    		int a          = int.Parse(str[2]);
    		String[] name  = new String [2];
    		if (method == 1 || method == 2) {
    			for (int i1 = 0; i1 < method; i1++) {
    				str      = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
    				name[i1] = str[0].Trim();
    				Np[i1]   = int.Parse(str[1]);
    			}
    			double[][][]x = new double [Np[0]][][];
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				x[i1] = new double [Np[1]][];
    				for (int i2 = 0; i2 < Np[1]; i2++)
    					x[i1][i2] = new double [N];
    			}
    			for (int i1 = 0; i1 < Np[1]; i1++) {
    				for (int i2 = 0; i2 < N; i2++) {
    					str  = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
    					for (int i3 = 0; i3 < Np[0]; i3++) {
    						x[i3][i1][i2] = double.Parse(str[i3]);
    					}
    				}
    			}
    			aov(method, Np, N, name, a, x);
    		}
    		else
    			Console.WriteLine("一元配置法,または,二元配置法だけです!");
    	}
    
    	/**************************************/
    	/* 分散分析(一元配置法と二元配置法) */
    	/*      method : 1 or 2               */
    	/*      Np[0] : 因子1の水準数         */
    	/*        [1] : 因子2の水準数         */
    	/*      N : データ数                  */
    	/*      name[0] : 因子1の名前         */
    	/*          [1] : 因子2の名前         */
    	/*      a : a %値                    */
    	/*      x : データ                    */
    	/*           coded by Y.Suganuma      */
    	/**************************************/
    	void aov(int method, int[] Np, int N, String[] name, int a, double[][][] x)
    	{
    		int sw = 0;
    					// 一元配置法
    		if (method == 1) {
    
    			double[] xi = new double [Np[0]];
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				xi[i1] = 0.0;
    				for (int i2 = 0; i2 < N; i2++)
    					xi[i1] += x[i1][0][i2];
    				xi[i1] /= N;
    			}
    
    			double xa = 0.0;
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				for (int i2 = 0; i2 < N; i2++)
    					xa += x[i1][0][i2];
    			}
    			xa /= (Np[0] * N);
    
    			double SP = 0.0;
    			for (int i1 = 0; i1 < Np[0]; i1++)
    				SP += (xi[i1] - xa) * (xi[i1] - xa);
    			SP *= N;
    
    			double SE = 0.0;
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				for (int i2 = 0; i2 < N; i2++)
    					SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
    			}
    
    			double VP = SP / (Np[0] - 1);
    			double VE = SE / (Np[0] * (N - 1));
    
    			double FP = VP / VE;
    
    			p         = 0.01 * a;
    			dof2      = Np[0] * (N - 1);
    
    			Console.WriteLine("全変動: 平方和 " + (SP+SE) + " 自由度 " + (Np[0]*N-1));
    			dof1 = Np[0] - 1;
    			Console.WriteLine(name[0] + "の水準間: 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "%値 " + p_f(ref sw));
    			Console.WriteLine("水準内: 平方和 " + SE + " 自由度 " + (Np[0]*(N-1))+ " 不偏分散 " + VE);
    		}
    					// 二元配置法
    		else {
    
    			double[] xi = new double [Np[0]];
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				xi[i1] = 0.0;
    				for (int i2 = 0; i2 < Np[1]; i2++) {
    					for (int i3 = 0; i3 < N; i3++)
    						xi[i1] += x[i1][i2][i3];
    				}
    				xi[i1] /= (Np[1] * N);
    			}
    
    			double[] xj = new double [Np[1]];
    			for (int i1 = 0; i1 < Np[1]; i1++) {
    				xj[i1] = 0.0;
    				for (int i2 = 0; i2 < Np[0]; i2++) {
    					for (int i3 = 0; i3 < N; i3++)
    						xj[i1] += x[i2][i1][i3];
    				}
    				xj[i1] /= (Np[0] * N);
    			}
    
    			double[][] xij = new double [Np[0]][];
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				xij[i1] = new double [Np[1]];
    				for (int i2 = 0; i2 < Np[1]; i2++) {
    					xij[i1][i2] = 0.0;
    					for (int i3 = 0; i3 < N; i3++)
    						xij[i1][i2] += x[i1][i2][i3];
    					xij[i1][i2] /= N;
    				}
    			}
    
    			double xa = 0.0;
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				for (int i2 = 0; i2 < Np[1]; i2++) {
    					for (int i3 = 0; i3 < N; i3++)
    						xa += x[i1][i2][i3];
    				}
    			}
    			xa /= (Np[0] * Np[1] * N);
    
    			double SP = 0.0;
    			for (int i1 = 0; i1 < Np[0]; i1++)
    				SP += (xi[i1] - xa) * (xi[i1] - xa);
    			SP *= (Np[1] * N);
    
    			double SQ = 0.0;
    			for (int i1 = 0; i1 < Np[1]; i1++)
    				SQ += (xj[i1] - xa) * (xj[i1] - xa);
    			SQ *= (Np[0] * N);
    
    			double SI = 0.0;
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				for (int i2 = 0; i2 < Np[1]; i2++)
    					SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
    			}
    			SI *= N;
    
    			double SE = 0.0;
    			for (int i1 = 0; i1 < Np[0]; i1++) {
    				for (int i2 = 0; i2 < Np[1]; i2++) {
    					for (int i3 = 0; i3 < N; i3++)
    						SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
    				}
    			}
    
    			double VP = SP / (Np[0] - 1);
    			double VQ = SQ / (Np[1] - 1);
    			double VI = SI / ((Np[0] - 1) * (Np[1] - 1));
    			double VE = SE / (Np[0] * Np[1] * (N - 1));
    
    			double FP = VP / VE;
    			double FQ = VQ / VE;
    			double FI = VI / VE;
    
    			p         = 0.01 * a;
    			dof2      = Np[0] * Np[1] * (N - 1);
    
    			Console.WriteLine("全変動: 平方和 " + (SP+SQ+SI+SE) + " 自由度 " + (Np[0]*Np[1]*N-1));
    			dof1 = Np[0] - 1;
    			Console.WriteLine(name[0] + "の水準間: 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "%値 " + p_f(ref sw));
    			dof1 = Np[1] - 1;
    			Console.WriteLine(name[1] + "の水準間: 平方和 " + SQ + " 自由度 " + (Np[1]-1) + " 不偏分散 " + VQ + " F " + FQ + " " + a + "%値 " + p_f(ref sw));
    			dof1 = (Np[0] - 1) * (Np[1] - 1);
    			Console.WriteLine("相互作用: 平方和 " + SI + " 自由度 " + ((Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + VI + " F " + FI + " " + a + "%値 " + p_f(ref sw));
    			Console.WriteLine("水準内: 平方和 " + SE + " 自由度 " + (Np[0]*Np[1]*(N-1)) + " 不偏分散 " + VE);
    		}
    	}
    	/****************************************/
    	/* f分布の計算(P(X = ff), P(X < ff)) */
    	/*      dd : P(X = ff)                  */
    	/*      df1,df2 : 自由度                */
    	/*      return : P(X < ff)              */
    	/****************************************/
    	double f(double ff, ref double dd, int df1, int df2)
    	{
    		double pp, u;
    		int ia, ib;
    
    		if (ff < 1.0e-10)
    			ff = 1.0e-10;
    
    		double x = ff * df1 / (ff * df1 + df2);
    
    		if (df1%2 == 0) {
    			if (df2%2 == 0) {
    				u  = x * (1.0 - x);
    				pp = x;
    				ia = 2;
    				ib = 2;
    			}
    			else {
    				u  = 0.5 * x * Math.Sqrt(1.0-x);
    				pp = 1.0 - Math.Sqrt(1.0-x);
    				ia = 2;
    				ib = 1;
    			}
    		}
    
    		else {
    			if (df2%2 == 0) {
    				u  = 0.5 * Math.Sqrt(x) * (1.0 - x);
    				pp = Math.Sqrt(x);
    				ia = 1;
    				ib = 2;
    			}
    			else {
    				u  = Math.Sqrt(x*(1.0-x)) / Math.PI;
    				pp = 1.0 - 2.0 * Math.Atan2(Math.Sqrt(1.0-x), Math.Sqrt(x)) / Math.PI;
    				ia = 1;
    				ib = 1;
    			}
    		}
    
    		if (ia != df1) {
    			for (int i1 = ia; i1 <= df1-2; i1 += 2) {
    				pp -= 2.0 * u / i1;
    				u  *= x * (i1 + ib) / i1;
    			}
    		}
    
    		if (ib != df2) {
    			for (int i1 = ib; i1 <= df2-2; i1 += 2) {
    				pp += 2.0 * u / i1;
    				u  *= (1.0 - x) * (i1 + df1) / i1;
    			}
    		}
    
    		dd = u / ff;
    
    		return pp;
    	}
    
    	/****************************************/
    	/* f分布のp%値(P(X > u) = 0.01p)   */
    	/*      ind : >= 0 : normal(収束回数) */
    	/*            = -1 : 収束しなかった     */
    	/****************************************/
    	double p_f(ref int ind)
    	{
    		int MAX  = 340;
    		double a = 0.0, a1 = 0.0, b = 0.0, b1 = 0.0, df1 = 0.0, df2 = 0.0, e = 0.0, ff = 0.0, yq = 0.0;
    		int sw   = 0;
    	/*
    	     初期値計算の準備
    	          while文は,大きな自由度によるガンマ関数の
    	          オーバーフローを避けるため
    	*/
    		while (sw >= 0) {
    
    			df1 = 0.5 * (dof1 - sw);
    			df2 = 0.5 * dof2;
    			a   = 2.0 / (9.0 * (dof1 - sw));
    			a1  = 1.0 - a;
    			b   = 2.0 / (9.0 * dof2);
    			b1  = 1.0 - b;
    
    			yq  = p_normal(ref ind);
    
    			e   = b1 * b1 - b * yq * yq;
    
    			if (e > 0.8 || (dof1+dof2-sw) <= MAX)
    				sw = -1;
    			else {
    				sw += 1;
    				if ((dof1-sw) == 0)
    					sw = -2;
    			}
    		}
    
    		if (sw == -2)
    			ind = -1;
    
    		else {
    	/*
    	     f0 : 初期値
    	*/
    			double f0 = 0.0;
    			if (e > 0.8) {
    				double x  = (a1 * b1 + yq * Math.Sqrt(a1*a1*b+a*e)) / e;
    				f0        = Math.Pow(x, 3.0);
    			}
    			else {
    				double y1 = Math.Pow((double)dof2, df2-1.0);
    				double y2 = Math.Pow((double)dof1, df2);
    				double x  = gamma(df1+df2, ref ind) / gamma(df1, ref ind) / gamma(df2, ref ind) *
    	                 2.0 * y1 / y2 / p;
    				f0 = Math.Pow(x, 2.0/dof2);
    			}
    	/*
    	     ニュートン法
    	*/
    			ff = newton(f0, 1.0e-6, 1.0e-10, 100, ref ind, f_f, f_df);
    		}
    
    		return ff;
    	}
    
    	/****************************************/
    	/* Γ(x)の計算(ガンマ関数,近似式) */
    	/*      ier : =0 : normal               */
    	/*            =-1 : x=-n (n=0,1,2,・・・)  */
    	/*      return : 結果                   */
    	/****************************************/
    	static double gamma(double x, ref int ier)
    	{
    		double g;
    		ier = 0;
    
    		if (x > 5.0) {
    			double v = 1.0 / x;
    			double s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
                            0.00078403922172) * v - 0.000229472093621) * v -
                            0.00268132716049) * v + 0.00347222222222) * v +
                            0.0833333333333) * v + 1.0;
    			g        = 2.506628274631001 * Math.Exp(-x) * Math.Pow(x,x-0.5) * s;
    		}
    
    		else {
    
    			double err = 1.0e-20;
    			double w   = x;
    			double t   = 1.0;
    
    			if (x < 1.5) {
    
    				if (x < err) {
    					int k    = (int)x;
    					double y = (double)k - x;
    					if (Math.Abs(y) < err || Math.Abs(1.0-y) < err)
    						ier = -1;
    				}
    
    				if (ier == 0) {
    					while (w < 1.5) {
    						t /= w;
    						w += 1.0;
    					}
    				}
    			}
    
    			else {
    				if (w > 2.5) {
    					while (w > 2.5) {
    						w -= 1.0;
    						t *= w;
    					}
    				}
    			}
    
    			w -= 2.0;
    			g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
                     0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
                     0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
                     0.999999926;
    			g *= t;
    		}
    
    		return g;
    	}
    
    	/********************************/
    	/* 1.0 - p - P(X > x)(関数値) */
    	/********************************/
    	double f_f(double x)
    	{
    		double y = 0.0;
    		return f(x, ref y, dof1, dof2) - 1.0 + p;
    	}
    
    	/**************************/
    	/* P(X = x)(関数の微分) */
    	/**************************/
    	double f_df(double x)
    	{
    		double y = 0.0;
    		double z = f(x, ref y, dof1, dof2);
    
    		return y;
    	}
    
    	/*****************************************************/
    	/* Newton法による非線形方程式(f(x)=0)の解            */
    	/*      x1 : 初期値                                  */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    	/*      max : 最大試行回数                           */
    	/*      ind : 実際の試行回数                         */
    	/*            (負の時は解を得ることができなかった) */
    	/*      fn : 関数値を計算する関数                    */
    	/*      dfn : 関数の微分値を計算する関数             */
    	/*      return : 解                                  */
    	/*****************************************************/
    	double newton(double x1, double eps1, double eps2, int max, ref int ind,
    	              Func fn, Func dfn)
    	{
    		double x = x1;
    		int sw   = 0;
    		ind      = 0;
    
    		while (sw == 0 && ind >= 0) {
    
    			ind++;
    			sw        = 1;
    			double g  = fn(x1);
    
    			if (Math.Abs(g) > eps2) {
    				if (ind <= max) {
    					double dg = dfn(x1);
    					if (Math.Abs(dg) > eps2) {
    						x = x1 - g / dg;
    						if (Math.Abs(x-x1) > eps1 && Math.Abs(x-x1) > eps1*Math.Abs(x)) {
    							x1 = x;
    							sw = 0;
    						}
    					}
    					else
    						ind = -1;
    				}
    				else
    					ind = -1;
    			}
    		}
    
    		return x;
    	}
    
    	/*************************************************/
    	/* 標準正規分布N(0,1)の計算(P(X = x), P(X < x))*/
    	/*      w : P(X = x)                             */
    	/*      return : P(X < x)                        */
    	/*************************************************/
    	double normal(double x, ref double w)
    	{
    			// 確率密度関数(定義式)
    		w = Math.Exp(-0.5 * x * x) / Math.Sqrt(2.0*Math.PI);
    			// 確率分布関数(近似式を使用)
    		double y = 0.70710678118654 * Math.Abs(x);
    		double z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
    		           y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
    		           y * 0.0000430638)))));
    		double P = 1.0 - Math.Pow(z,-16.0);
    
    		if (x < 0.0)
    			P = 0.5 - 0.5 * P;
    		else
    			P = 0.5 + 0.5 * P;
    
    		return P;
    	}
    
    	/******************************************************************/
    	/* 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用) */
    	/*      ind : >= 0 : normal(収束回数)                           */
    	/*            = -1 : 収束しなかった                               */
    	/******************************************************************/
    	double p_normal(ref int ind)
    	{
    		int sw   = 0;
    		double u = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, ref sw, normal_f);
    		ind      = sw;
    
    		return u;
    	}
    
    	/******************************/
    	/* 1.0 - p - P(X>x)(関数値) */
    	/******************************/
    	double normal_f(double x)
    	{
    		double y = 0.0;
    		return 1.0 - p - normal(x, ref y);
    	}
    
    	/*********************************************************/
    	/* 二分法による非線形方程式(f(x)=0)の解                  */
    	/*      x1,x2 : 初期値                                   */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)           */
    	/*      max : 最大試行回数                               */
    	/*      ind : 実際の試行回数                             */
    	/*            (負の時は解を得ることができなかった)     */
    	/*      fn : 関数値を計算する関数                        */
    	/*      return : 解                                      */
    	/*********************************************************/
    	double bisection(double x1, double x2, double eps1, double eps2, int max,
    	                 ref int ind, Func fn)
    	{
    		double f1 = fn(x1);
    		double f2 = fn(x2);
    		double x0 = 0.0;
    
    		if (f1*f2 > 0.0)
    			ind = -1;
    
    		else {
    			ind = 0;
    			if (f1*f2 == 0.0)
    				x0 = (f1 == 0.0) ? x1 : x2;
    			else {
    				int sw = 0;
    				while (sw == 0 && ind >= 0) {
    					sw        = 1;
    					ind      += 1;
    					x0        = 0.5 * (x1 + x2);
    					double f0 = fn(x0);
    					if (Math.Abs(f0) > eps2) {
    						if (ind <= max) {
    							if (Math.Abs(x1-x2) > eps1 && Math.Abs(x1-x2) > eps1*Math.Abs(x2)) {
    								sw = 0;
    								if (f0*f1 < 0.0) {
    									x2 = x0;
    									f2 = f0;
    								}
    								else {
    									x1 = x0;
    									f1 = f0;
    								}
    							}
    						}
    						else
    							ind = -1;
    					}
    				}
    			}
    		}
    
    		return x0;
    	}
    }
    
    /*
    -------- 一元配置法に対するデータ例(コメント部分を除いて下さい)--------
    1 6 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    工場 3   // 因子の名前 水準の数
    3.1 4.7 5.1   // 各水準に対する1番目のデータ
    4.1 5.6 3.7   // 各水準に対する2番目のデータ
    3.3 4.3 4.5   // 各水準に対する3番目のデータ
    3.9 5.9 6.0   // 各水準に対する4番目のデータ
    3.7 6.1 3.9   // 各水準に対する5番目のデータ
    2.4 4.2 5.4   // 各水準に対する6番目のデータ
    
    -------- 二元配置法に対するデータ例(コメント部分を除いて下さい)--------
    2 3 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    薬剤 5   // 1番目の因子の名前 その水準の数
    品種 2   // 2番目の因子の名前 その水準の数
    3 4 12 -4 -4   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する1番目のデータ
    8 -8 31 12 19   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する2番目のデータ
    7 -5 8 0 23   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する3番目のデータ
    8 -10 9 10 15   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する1番目のデータ
    -5 11 26 -1 13   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する2番目のデータ
    10 -6 13 -7 -6   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する3番目のデータ
    */
    			

  8. VB

    '**************************'
    ' 分散分析                 '
    '      coded by Y.Suganuma '
    '**************************'
    Imports System.Text.RegularExpressions
    
    Module Test
    
    	Dim p As Double     ' α%値を計算するとき時α/100を設定
    	Dim dof1 As Integer     ' 自由度
    	Dim dof2 As Integer     ' 自由度
    
    	Sub Main()
    
    		Dim Np() As Integer = {1, 1}
    		Dim MS As Regex     = New Regex("\s+") 
    
    		Dim str() As String   = MS.Split(Console.ReadLine().Trim())
    		Dim method As Integer = Integer.Parse(str(0))
    		Dim N As Integer      = Integer.Parse(str(1))
    		Dim a As Integer      = Integer.Parse(str(2))
    		Dim name(2) As String
    		If method = 1 or method = 2
    			For i1 As Integer = 0 To method-1
    				str      = MS.Split(Console.ReadLine().Trim())
    				name(i1) = str(0)
    				Np(i1)   = Integer.Parse(str(1))
    			Next
    			Dim x(Np(0),Np(1),N) As Double
    			For i1 As Integer = 0 To Np(1)-1
    				For i2 As Integer = 0 To N-1
    					str = MS.Split(Console.ReadLine().Trim())
    					For i3 As Integer = 0 To Np(0)-1
    						x(i3,i1,i2) = Double.Parse(str(i3))
    					Next
    				Next
    			Next
    			aov(method, Np, N, name, a, x)
    		Else
    			Console.WriteLine("一元配置法,または,二元配置法だけです!")
    		End If
    
    	End Sub
    
    	'************************************'
    	' 分散分析(一元配置法と二元配置法) '
    	'      method : 1 or 2               '
    	'      Np(0) : 因子1の水準数         '
    	'        (1) : 因子2の水準数         '
    	'      N : データ数                  '
    	'      name(0) : 因子1の名前         '
    	'          (1) : 因子2の名前         '
    	'      a : a %値                    '
    	'      x : データ                    '
    	'           coded by Y.Suganuma      '
    	'************************************'
    	Sub aov(method As Integer, Np() As Integer, N As Integer, name() As String, a As Integer, x(,,) As Double)
    
    					' P(X > x) - 1 + p(関数値)(ラムダ式)
    		Dim f_f = Function(v1) As Double
    			Dim v2 As Double = 0.0
    			Return f(v1, v2, dof1, dof2) - 1.0 + p
    		End Function
    					' P(X = x)(関数の微分)(ラムダ式)
    		Dim f_df = Function(v1) As Double
    			Dim v2 As Double = 0.0
    			Dim v3 As Double = f(v1, v2, dof1, dof2)
    			Return v2
    		End Function
    					' 1.0 - p - P(X>x)(関数値)(ラムダ式)
    		Dim normal_f = Function(v1) As Double
    			Dim v2 As Double = 0.0
    			Return 1.0 - p - normal(v1, v2)
    		End Function
    
    		Dim sw As Integer = 0
    					' 一元配置法
    		If method = 1
    
    			Dim xi(Np(0)) As Double
    			For i1 As Integer = 0 To Np(0)-1
    				xi(i1) = 0.0
    				For i2 As Integer = 0 To N-1
    					xi(i1) += x(i1,0,i2)
    				Next
    				xi(i1) /= N
    			Next
    
    			Dim xa As Double = 0.0
    			For i1 As Integer = 0 To Np(0)-1
    				For i2 As Integer = 0 To N-1
    					xa += x(i1,0,i2)
    				Next
    			Next
    			xa /= (Np(0) * N)
    
    			Dim SP As Double = 0.0
    			For i1 As Integer = 0 To Np(0)-1
    				SP += (xi(i1) - xa) * (xi(i1) - xa)
    			Next
    			SP *= N
    
    			Dim SE As Double = 0.0
    			For i1 As Integer = 0 To Np(0)-1
    				For i2 As Integer = 0 To N-1
    					SE += (x(i1,0,i2) - xi(i1)) * (x(i1,0,i2) - xi(i1))
    				Next
    			Next
    
    			Dim VP As Double = SP / (Np(0) - 1)
    			Dim VE As Double = SE / (Np(0) * (N - 1))
    
    			Dim FP As Double = VP / VE
    
    			p         = 0.01 * a
    			dof2      = Np(0) * (N - 1)
    
    			Console.WriteLine("全変動: 平方和 " & (SP+SE) & " 自由度 " & (Np(0)*N-1))
    			dof1 = Np(0) - 1
    			Console.WriteLine(name(0) & "の水準間: 平方和 " & SP & " 自由度 " & (Np(0)-1) & " 不偏分散 " & VP & " F " & FP & " " & a & "%値 " & p_f(sw, f_f, f_df, normal_f))
    			Console.WriteLine("水準内: 平方和 " & SE & " 自由度 " & (Np(0)*(N-1)) & " 不偏分散 " & VE)
    					' 二元配置法
    		Else
    
    			Dim xi(Np(0)) As Double
    			For i1 As Integer = 0 To Np(0)-1
    				xi(i1) = 0.0
    				For i2 As Integer = 0 To Np(1)-1
    					For i3 As Integer = 0 To N-1
    						xi(i1) += x(i1,i2,i3)
    					Next
    				Next
    				xi(i1) /= (Np(1) * N)
    			Next
    
    			Dim xj(Np(1)) As Double
    			For i1 As Integer = 0 To Np(1)-1
    				xj(i1) = 0.0
    				For i2 As Integer = 0 To Np(0)-1
    					For i3 As Integer = 0 To N-1
    						xj(i1) += x(i2,i1,i3)
    					Next
    				Next
    				xj(i1) /= (Np(0) * N)
    			Next
    
    			Dim xij(Np(0),Np(1)) As Double
    			For i1 As Integer = 0 To Np(0)-1
    				For i2 As Integer = 0 To Np(1)-1
    					xij(i1,i2) = 0.0
    					For i3 As Integer = 0 To N-1
    						xij(i1,i2) += x(i1,i2,i3)
    					Next
    					xij(i1,i2) /= N
    				Next
    			Next
    
    			Dim xa As Double = 0.0
    			For i1 As Integer = 0 To Np(0)-1
    				For i2 As Integer = 0 To Np(1)-1
    					For i3 As Integer = 0 To N-1
    						xa += x(i1,i2,i3)
    					Next
    				Next
    			Next
    			xa /= (Np(0) * Np(1) * N)
    
    			Dim SP As Double = 0.0
    			For i1 As Integer = 0 To Np(0)-1
    				SP += (xi(i1) - xa) * (xi(i1) - xa)
    			Next
    			SP *= (Np(1) * N)
    
    			Dim SQ As Double = 0.0
    			For i1 As Integer = 0 To Np(1)-1
    				SQ += (xj(i1) - xa) * (xj(i1) - xa)
    			Next
    			SQ *= (Np(0) * N)
    
    			Dim SI As Double = 0.0
    			For i1 As Integer = 0 To Np(0)-1
    				For i2 As Integer = 0 To Np(1)-1
    					SI += (xij(i1,i2) - xi(i1) - xj(i2) + xa) * (xij(i1,i2) - xi(i1) - xj(i2) + xa)
    				Next
    			Next
    			SI *= N
    
    			Dim SE As Double = 0.0
    			For i1 As Integer = 0 To Np(0)-1
    				For i2 As Integer = 0 To Np(1)-1
    					For i3 As Integer = 0 To N-1
    						SE += (x(i1,i2,i3) - xij(i1,i2)) * (x(i1,i2,i3) - xij(i1,i2))
    					Next
    				Next
    			Next
    
    			Dim VP As Double = SP / (Np(0) - 1)
    			Dim VQ As Double = SQ / (Np(1) - 1)
    			Dim VI As Double = SI / ((Np(0) - 1) * (Np(1) - 1))
    			Dim VE As Double = SE / (Np(0) * Np(1) * (N - 1))
    
    			Dim FP As Double = VP / VE
    			Dim FQ As Double = VQ / VE
    			Dim FI As Double = VI / VE
    
    			p    = 0.01 * a
    			dof2 = Np(0) * Np(1) * (N - 1)
    
    			Console.WriteLine("全変動: 平方和 " & (SP+SQ+SI+SE) & " 自由度 " & (Np(0)*Np(1)*N-1))
    			dof1 = Np(0) - 1
    			Console.WriteLine(name(0) & "の水準間: 平方和 " & SP & " 自由度 " & (Np(0)-1) & " 不偏分散 " & VP & " F " & FP & " " & a & "%値 " & p_f(sw, f_f, f_df, normal_f))
    			dof1 = Np(1) - 1
    			Console.WriteLine(name(1) & "の水準間: 平方和 " & SQ & " 自由度 " & (Np(1)-1) & " 不偏分散 " & VQ & " F " & FQ & " " & a & "%値 " & p_f(sw, f_f, f_df, normal_f))
    			dof1 = (Np(0) - 1) * (Np(1) - 1)
    			Console.WriteLine("相互作用: 平方和 " & SI & " 自由度 " & ((Np(0)-1)*(Np(1)-1)) & " 不偏分散 " & VI & " F " & FI & " " & a & "%値 " & p_f(sw, f_f, f_df, normal_f))
    			Console.WriteLine("水準内: 平方和 " & SE & " 自由度 " & (Np(0)*Np(1)*(N-1)) & " 不偏分散 " & VE)
    		End If
    
    	End Sub
    
    	'**************************************'
    	' f分布の計算(P(X = ff), P(X < ff)) '
    	'      dd : P(X = ff)                  '
    	'      df1,df2 : 自由度                '
    	'      return : P(X < ff)              '
    	'**************************************'
    	Function f(ff As Double, ByRef dd As Double, df1 As Integer, df2 As Integer)
    
    		Dim pp As Double
    		Dim u As Double
    		Dim ia As Integer
    		Dim ib As Integer
    
    		If ff < 1.0e-10
    			ff = 1.0e-10
    		End If
    
    		Dim x As Double = ff * df1 / (ff * df1 + df2)
    
    		If (df1 Mod 2) = 0
    			If (df2 Mod 2) = 0
    				u  = x * (1.0 - x)
    				pp = x
    				ia = 2
    				ib = 2
    			Else
    				u  = 0.5 * x * Math.Sqrt(1.0-x)
    				pp = 1.0 - Math.Sqrt(1.0-x)
    				ia = 2
    				ib = 1
    			End If
    
    		Else
    			If (df2 Mod 2) = 0
    				u  = 0.5 * Math.Sqrt(x) * (1.0 - x)
    				pp = Math.Sqrt(x)
    				ia = 1
    				ib = 2
    			Else
    				u  = Math.Sqrt(x*(1.0-x)) / Math.PI
    				pp = 1.0 - 2.0 * Math.Atan2(Math.Sqrt(1.0-x), Math.Sqrt(x)) / Math.PI
    				ia = 1
    				ib = 1
    			End If
    		End If
    
    		If ia <> df1
    			For i1 As Integer = ia To df1-2 Step 2
    				pp -= 2.0 * u / i1
    				u  *= x * (i1 + ib) / i1
    			Next
    		End If
    
    		If ib <> df2
    			For i1 As Integer = ib To df2-2 Step 2
    				pp += 2.0 * u / i1
    				u  *= (1.0 - x) * (i1 + df1) / i1
    			Next
    		End If
    
    		dd = u / ff
    
    		Return pp
    
    	End Function
    
    	'**************************************'
    	' f分布のp%値(P(X > u) = 0.01p)   '
    	'      ind : >= 0 : normal(収束回数) '
    	'            = -1 : 収束しなかった     '
    	'**************************************'
    	Function p_f(ByRef ind As Integer, f_f As Func(Of Double, Double),
    	             f_df As Func(Of Double, Double), normal_f As Func(Of Double, Double))
    
    		Dim MAX As Integer = 340
    		Dim a As Double    = 0.0
    		Dim a1 As Double   = 0.0
    		Dim b As Double    = 0.0
    		Dim b1 As Double   = 0.0
    		Dim df1 As Double  = 0.0
    		Dim df2 As Double  = 0.0
    		Dim e As Double    = 0.0
    		Dim ff As Double   = 0.0
    		Dim yq As Double   = 0.0
    		Dim sw As Integer  = 0
    				'
    				' 初期値計算の準備
    				'      while文は,大きな自由度によるガンマ関数の
    				'      オーバーフローを避けるため
    				'
    		Do While sw >= 0
    
    			df1 = 0.5 * (dof1 - sw)
    			df2 = 0.5 * dof2
    			a   = 2.0 / (9.0 * (dof1 - sw))
    			a1  = 1.0 - a
    			b   = 2.0 / (9.0 * dof2)
    			b1  = 1.0 - b
    
    			yq  = p_normal(ind, normal_f)
    
    			e   = b1 * b1 - b * yq * yq
    
    			If e > 0.8 or (dof1+dof2-sw) <= MAX
    				sw = -1
    			Else
    				sw += 1
    				If (dof1-sw) = 0
    					sw = -2
    				End If
    			End If
    		Loop
    
    		If sw = -2
    			ind = -1
    
    		Else
    				'
    				' f0 : 初期値
    				'
    			Dim f0 As Double = 0.0
    			If e > 0.8
    				Dim x As Double  = (a1 * b1 + yq * Math.Sqrt(a1*a1*b+a*e)) / e
    				f0               = Math.Pow(x, 3.0)
    			Else
    				Dim y1 As Double = Math.Pow(dof2, df2-1.0)
    				Dim y2 As Double = Math.Pow(dof1, df2)
    				Dim x As Double  = gamma(df1+df2, ind) / gamma(df1, ind) / 
    				                   gamma(df2, ind) * 2.0 * y1 / y2 / p
    				f0 = Math.Pow(x, 2.0/dof2)
    			End If
    				'
    				' ニュートン法
    				'
    			ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind, f_f, f_df)
    		End If
    
    		Return ff
    
    	End Function
    
    	'**************************************'
    	' Γ(x)の計算(ガンマ関数,近似式) '
    	'      ier : =0 : normal               '
    	'            =-1 : x=-n (n=0,1,2,・・・)  '
    	'      return : 結果                   '
    	'**************************************'
    	Function gamma(x As Double, ByRef ier As Integer)
    
    		Dim g As Double
    		ier = 0
    
    		If x > 5.0
    			Dim v As Double = 1.0 / x
    			Dim s As Double = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
    			                  0.00078403922172) * v - 0.000229472093621) * v -
    			                  0.00268132716049) * v + 0.00347222222222) * v +
    			                  0.0833333333333) * v + 1.0
    			g = 2.506628274631001 * Math.Exp(-x) * Math.Pow(x,x-0.5) * s
    
    		Else
    			Dim err As Double = 1.0e-20
    			Dim w As Double   = x
    			Dim t As Double   = 1.0
    
    			If x < 1.5
    
    				If x < err
    					Dim k As Integer = Math.Floor(x)
    					Dim y As Double  = k - x
    					If Math.Abs(y) < err or Math.Abs(1.0-y) < err
    						ier = -1
    					End If
    				End If
    
    				If ier = 0
    					Do While w < 1.5
    						t /= w
    						w += 1.0
    					Loop
    				End If
    
    			Else
    				If w > 2.5
    					Do While w > 2.5
    						w -= 1.0
    						t *= w
    					Loop
    				End If
    			End If
    
    			w -= 2.0
    			g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
    			     0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
    			     0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
    			     0.999999926
    			g *= t
    		End If
    
    		Return g
    
    	End Function
    
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' Newton法による非線形方程式(f(x)=0)の解            '
    	'      x1 : 初期値                                  '
    	'      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   '
    	'      eps2 : 終了条件2(|f(x(k))|<eps2)       '
    	'      max : 最大試行回数                           '
    	'      ind : 実際の試行回数                         '
    	'            (負の時は解を得ることができなかった) '
    	'      fn : 関数値を計算する関数                    '
    	'      dfn : 関数の微分値を計算する関数             '
    	'      return : 解                                  '
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function newton(x1 As Double, eps1 As Double, eps2 As Double, max As Integer,
    	              ByRef ind As Integer,
    	              fn As Func(Of Double, Double), dfn As Func(Of Double, Double))
    
    		Dim x As Double   = x1
    		Dim sw As Integer = 0
    		ind = 0
    
    		Do While sw = 0 and ind >= 0
    
    			ind += 1
    			sw   = 1
    			Dim g As Double = fn(x1)
    
    			If Math.Abs(g) > eps2
    				If ind <= max
    					Dim dg As Double = dfn(x1)
    					If Math.Abs(dg) > eps2
    						x = x1 - g / dg
    						If Math.Abs(x-x1) > eps1 and Math.Abs(x-x1) > eps1*Math.Abs(x)
    							x1 = x
    							sw = 0
    						End If
    					Else
    						ind = -1
    					End If
    				Else
    					ind = -1
    				End If
    			End If
    		Loop
    
    		Return x
    
    	End Function
    
    	'***********************************************'
    	' 標準正規分布N(0,1)の計算(P(X = x), P(X < x))'
    	'      w : P(X = x)                             '
    	'      return : P(X < x)                        '
    	'***********************************************'
    	Function normal(x As Double, ByRef w As Double)
    
    			' 確率密度関数(定義式)
    		w = Math.Exp(-0.5 * x * x) / Math.Sqrt(2.0*Math.PI)
    			' 確率分布関数(近似式を使用)
    		Dim y As Double = 0.70710678118654 * Math.Abs(x)
    		Dim z As Double = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
    		                  y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
    		                  y * 0.0000430638)))))
    		Dim P As Double = 1.0 - Math.Pow(z,-16.0)
    
    		If x < 0.0
    			P = 0.5 - 0.5 * P
    		Else
    			P = 0.5 + 0.5 * P
    		End If
    
    		Return P
    
    	End Function
    
    	'****************************************************************'
    	' 標準正規分布N(0,1)のp%値(P(X > u) = 0.01p)(二分法を使用) '
    	'      ind : >= 0 : normal(収束回数)                           '
    	'            = -1 : 収束しなかった                               '
    	'****************************************************************'
    	Function p_normal(ByRef ind As Integer, fn As Func(Of Double, Double))
    
    		Dim sw As Integer = 0
    		Dim u As Double   = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, sw, fn)
    		ind               = sw
    
    		Return u
    
    	End Function
    
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' 二分法による非線形方程式(f(x)=0)の解              '
    	'      x1,x2 : 初期値                               '
    	'      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   '
    	'      eps2 : 終了条件2(|f(x(k))|<eps2)       '
    	'      max : 最大試行回数                           '
    	'      ind : 実際の試行回数                         '
    	'            (負の時は解を得ることができなかった) '
    	'      fn : 関数値を計算する関数                    '
    	'      return : 解                                  '
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function bisection(x1 As Double, x2 As Double, eps1 As Double, eps2 As Double,
    	                   max As Integer, ByRef ind As Integer,
    	                   fn As Func(Of Double, Double))
    
    		Dim f1 As Double = fn(x1)
    		Dim f2 As Double = fn(x2)
    		Dim x0 As Double = 0.0
    
    		If f1*f2 > 0.0
    			ind = -1
    
    		Else
    			ind = 0
    			If f1*f2 = 0.0
    				If f1 = 0.0
    					x0 = x1
    				Else
    					x0 = x2
    				End If
    			Else
    				Dim sw As Integer = 0
    				Do While sw = 0 and ind >= 0
    					sw   = 1
    					ind += 1
    					x0   = 0.5 * (x1 + x2)
    					Dim f0 As Double = fn(x0)
    					If Math.Abs(f0) > eps2
    						If ind <= max
    							If Math.Abs(x1-x2) > eps1 and Math.Abs(x1-x2) > eps1*Math.Abs(x2)
    								sw = 0
    								If f0*f1 < 0.0
    									x2 = x0
    									f2 = f0
    								Else
    									x1 = x0
    									f1 = f0
    								End If
    							End If
    						Else
    							ind = -1
    						End If
    					End If
    				Loop
    			End If
    		End If
    
    		Return x0
    
    	End Function
    
    End Module
    
    /*
    -------- 一元配置法に対するデータ例(コメント部分を除いて下さい)--------
    1 6 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    工場 3   // 因子の名前 水準の数
    3.1 4.7 5.1   // 各水準に対する1番目のデータ
    4.1 5.6 3.7   // 各水準に対する2番目のデータ
    3.3 4.3 4.5   // 各水準に対する3番目のデータ
    3.9 5.9 6.0   // 各水準に対する4番目のデータ
    3.7 6.1 3.9   // 各水準に対する5番目のデータ
    2.4 4.2 5.4   // 各水準に対する6番目のデータ
    
    -------- 二元配置法に対するデータ例(コメント部分を除いて下さい)--------
    2 3 5   // 因子の数 各水準におけるデータ数 有意水準(%)
    薬剤 5   // 1番目の因子の名前 その水準の数
    品種 2   // 2番目の因子の名前 その水準の数
    3 4 12 -4 -4   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する1番目のデータ
    8 -8 31 12 19   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する2番目のデータ
    7 -5 8 0 23   // 1番目の因子の各水準に対して,2番目の因子の水準1に対する3番目のデータ
    8 -10 9 10 15   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する1番目のデータ
    -5 11 26 -1 13   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する2番目のデータ
    10 -6 13 -7 -6   // 1番目の因子の各水準に対して,2番目の因子の水準2に対する3番目のデータ
    */
    			

静岡理工科大学 菅沼ホーム 目次 索引