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

F 分布

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

  プログラムは,グラフ出力を指定すると F 分布の密度関数および分布関数の値を指定した範囲だけファイルに出力します.また,グラフ出力を指定しないと,指定された値における密度関数および分布関数の値,または,%値を出力します.

  1. C++

    /****************************/
    /* F分布の計算             */
    /*      coded by Y.Suganuma */
    /****************************/
    #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;   // 自由度
                      // 上の3つの変数は,必ず,この位置で定義しておく必要がある
    
    int main()
    {
    	double h, x, x1, x2, y, z;
    	int sw;
    	char file1[100], file2[100];
    	FILE *out1, *out2;
    
    	printf("自由度1は? ");
    	scanf("%d", &dof1);
    	printf("自由度2は? ");
    	scanf("%d", &dof2);
    	printf("目的とする結果は? \n");
    	printf("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)\n");
    	printf("     =1 : p%値( P(X > u) = 0.01p となるuの値) ");
    	scanf("%d", &sw);
    
    	if (sw == 0) {
    
    		printf("グラフ出力?(=1: yes,  =0: no) ");
    		scanf("%d", &sw);
    					// 密度関数と分布関数の値
    		if (sw == 0) {
    			printf("   データは? ");
    			scanf("%lf", &x);
    			y = f(x, &z, dof1, dof2);
    			printf("P(X = %f) = %f,  P(X < %f) = %f (自由度 = %d, %d)\n",
                       x, z, x, y, dof1, dof2);
    		}
    					// グラフ出力
    		else {
    			printf("   密度関数のファイル名は? ");
    			scanf("%s", file1);
    			printf("   分布関数のファイル名は? ");
    			scanf("%s", file2);
    			out1 = fopen(file1,"w");
    			out2 = fopen(file2,"w");
    			printf("   データの下限は? ");
    			scanf("%lf", &x1);
    			printf("   データの上限は? ");
    			scanf("%lf", &x2);
    			printf("   刻み幅は? ");
    			scanf("%lf", &h);
    			for (x = x1; x < x2+0.5*h; x += h) {
    				y = f(x, &z, dof1, dof2);
    				fprintf(out1, "%f %f\n", x, z);
    				fprintf(out2, "%f %f\n", x, y);
    			}
    		}
    	}
    					// %値
    	else {
    		printf("%の値は? ");
    		scanf("%lf", &x);
    		p = 0.01 * x;
    		if (p < 1.0e-7)
    			printf("%f%値 = ∞(自由度 = %d, %d)\n", x, dof1, dof2);
    		else if ((1.0-p) < 1.0e-7)
    			printf("%f%値 = 0(自由度 = %d, %d)\n", x, dof1, dof2);
    		else {
    			y = p_f(&sw);
    			printf("%f%値 = %f  sw %d (自由度 = %d, %d)\n",
                       x, y, sw, dof1, dof2);
    		}
    	}
    
    	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)の解            */
    /*      f : f(x)を計算する関数名                     */
    /*      df : 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)の解                  */
    /*      f : 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;
    }
    			

  2. Java

      グラフ出力を指定すると,指定された範囲の密度関数および分布関数の値をファイルに出力すると共にグラフも表示しますが,グラフ出力を指定しないと,指定された値における密度関数および分布関数の値,または,%値を出力します.アプレット版では,任意のデータに対する処理を画面上で実行することが可能であり,結果はテキストエリアに出力されると共に,「確率(複数点)」を選択すればグラフも表示されます.なお,グラフに関しては,「グラフの表示」の表示を参照して下さい.
    /****************************/
    /* F分布の計算             */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.io.*;
    import java.text.*;
    import java.awt.*;
    import javax.swing.*;
    import java.awt.event.*;
    import java.util.*;
    
    public class Test {
    
    	static double p;   // α%値を計算するとき時α/100を設定
    	static int dof1, dof2;    // 自由度
    
    	public static void main(String args[]) throws IOException
    	{
    		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
    		double h, x, to, y;
    		double z[] = new double [1];
    		int sw, sw1[] = new int [1];
    
    		System.out.print("自由度1は? ");
    		dof1 = Integer.parseInt(in.readLine());
    		System.out.print("自由度2は? ");
    		dof2 = Integer.parseInt(in.readLine());
    		System.out.println("目的とする結果は? ");
    		System.out.print("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)\n");
    		System.out.print("     =1 : p%値( P(X > u) = 0.01p となるuの値) ");
    		sw = Integer.parseInt(in.readLine());
    
    		if (sw == 0) {
    
    			System.out.print("グラフ出力?(=1: yes,  =0: no) ");
    			sw = Integer.parseInt(in.readLine());
    					// 密度関数と分布関数の値
    			if (sw == 0) {
    				System.out.print("   データは? ");
    				x = Double.parseDouble(in.readLine());
    				y = App.f(x, z, dof1, dof2);
    				System.out.println("P(X = " + x + ") = " + z[0] + ",  P( X < " + x + ") = " + y +
                                       " (自由度 = " + dof1 + " " + dof2 + ")");
    			}
    					// グラフ出力
    			else {
    				String file1, file2;
    				System.out.print("   密度関数のファイル名は? ");
    				file1 = in.readLine();
    				System.out.print("   分布関数のファイル名は? ");
    				file2 = in.readLine();
    				PrintStream out1 = new PrintStream(new FileOutputStream(file1));
    				PrintStream out2 = new PrintStream(new FileOutputStream(file2));
    				System.out.print("   データの上限は? ");
    				to = Double.parseDouble(in.readLine());
    				System.out.print("   刻み幅は? ");
    				h = Double.parseDouble(in.readLine());
    							// データ取得
    				ArrayList <Double> x1 = new ArrayList <Double> ();
    				ArrayList <Double> y1 = new ArrayList <Double> ();
    				ArrayList <Double> y2 = new ArrayList <Double> ();
    				for (x = 0.0; x < to+0.5*h; x += h) {
    					y = App.f(x, z, dof1, dof2);
    					out1.println(x + " " + z[0]);
    					out2.println(x + " " + y);
    					x1.add(x);
    					y1.add(z[0]);
    					y2.add(y);
    				}
    							// グラフの描画
    				graph(x1, y1, y2);
    			}
    		}
    					// %値
    		else {
    			System.out.print("%の値は? ");
    			x = Double.parseDouble(in.readLine());
    			p = 0.01 * x;
    			if (p < 1.0e-7)
    				System.out.println(x + "%値 = ∞ (自由度 = " + dof1 + " " + dof2 + ")");
    			else if ((1.0-p) < 1.0e-7)
    				System.out.println(x + "%値 = 0 (自由度 = " + dof1 + " " + dof2 + ")");
    			else {
    				y = App.p_f(sw1);
    				System.out.println(x + "%値 = " + y + "  sw " + sw1[0] +
                                       " (自由度 = " + dof1 + " " + dof2 + ")");
    			}
    		}
    	}
    
    	/*************************/
    	/* グラフの描画          */
    	/*      x1 : x座標データ */
    	/*      y1 : 密度関数    */
    	/*      y2 : 分布関数    */
    	/*************************/
    	static void graph(ArrayList <Double> x1, ArrayList <Double> y1, ArrayList <Double> y2)
    	{
    					// 密度関数
    		String title1[];   // グラフ,x軸,及び,y軸のタイトル
    		String g_title1[];   // 凡例(グラフの内容)
    		double x_scale[];   // x軸目盛り
    		double y_scale1[];   // y軸目盛り
    		double data_x[][], data1_y[][];   // データ
    		int n = x1.size();
    							// グラフ,x軸,及び,y軸のタイトル
    		title1 = new String [3];
    		title1[0] = "密度関数(F分布 自由度 = " + dof1 + ", " + dof2 + ")";
    		title1[1] = "x";
    		title1[2] = "f(x)";
    							// 凡例
    		g_title1 = new String [1];
    		g_title1[0] = "密度関数";
    							// x軸目盛り
    		x_scale = new double[3];
    		double x_max = (x1.get(n-1)).doubleValue();
    		double x_min = 0.0;
    		double x_step = (x_max - x_min) / 10;
    		int x_p = 0;
    		boolean ok = true;
    		int k = 0;
    		while (ok) {
    			if (x_step < 1.0) {
    				x_step *= 10;
    				k++;
    			}
    			else if (x_step >= 10.0) {
    				x_step /= 10;
    				k--;
    			}
    			else {
    				ok = false;
    				if (x_step-(int)x_step > 1.0e-5)
    					x_step = (int)x_step + 1;
    				else
    					x_step = (int)x_step;
    				if (k != 0) {
    					x_step = x_step * Math.pow(10, -k);
    					if (k > 0)
    						x_p = k;
    				}
    				double t = x_min;
    				while (t < x_max-0.001*x_step)
    					t += x_step;
    				x_max = t;
    			}
    		}
    		x_scale[0] = x_min;   // 最小値
    		x_scale[1] = x_max;   // 最大値
    		x_scale[2] = x_step;   // 刻み幅
    							// y軸目盛り
    		y_scale1 = new double[3];
    		double y_max = 0.0;
    		for (int i1 = 0; i1 < n; i1++) {
    			if ((y1.get(i1)).doubleValue() > y_max)
    				y_max = (y1.get(i1)).doubleValue();
    		}
    		double y_step = y_max / 5;
    		int y_p1 = 0;
    		ok = true;
    		k  = 0;
    		while (ok) {
    			if (y_step < 1.0) {
    				y_step *= 10;
    				k++;
    			}
    			else if (y_step >= 10.0) {
    				y_step /= 10;
    				k--;
    			}
    			else {
    				ok = false;
    				if (y_step-(int)y_step > 1.0e-5)
    					y_step = (int)y_step + 1;
    				else
    					y_step = (int)y_step;
    				if (k != 0) {
    					y_step = y_step * Math.pow(10, -k);
    					if (k > 0)
    						y_p1 = k;
    				}
    				double t = 0.0;
    				while (t < y_max-0.001*y_step)
    					t += y_step;
    				y_max = t;
    			}
    		}
    		y_scale1[0] = 0.0;   // 最小値
    		y_scale1[1] = y_max;   // 最大値
    		y_scale1[2] = y_step;   // 刻み幅
    							// データ
    		data_x = new double [1][n];
    		data1_y = new double [1][n];
    		for (int i1 = 0; i1 < n; i1++)
    			data_x[0][i1] = (x1.get(i1)).doubleValue();  
    		for (int i1 = 0; i1 < n; i1++)
    			data1_y[0][i1] = (y1.get(i1)).doubleValue();  
    							// 作図
    		LineGraph gp1 = new LineGraph(title1, g_title1, x_scale, x_p, y_scale1, y_p1, data_x, data1_y, true, false);
    					// 分布関数
    		String title2[];   // グラフ,x軸,及び,y軸のタイトル
    		String g_title2[];   // 凡例(グラフの内容)
    		double y_scale2[];   // y軸目盛り
    		double data2_y[][];   // データ
    							// グラフ,x軸,及び,y軸のタイトル
    		title2 = new String [3];
    		title2[0] = "分布関数(F分布 自由度 = " + dof1 + ", " + dof2 + ")";
    		title2[1] = "x";
    		title2[2] = "F(x)";
    							// 凡例
    		g_title2 = new String [1];
    		g_title2[0] = "分布関数";
    							// y軸目盛り
    		y_scale2 = new double[3];
    		int y_p2 = 1;
    		y_scale2[0] = 0.0;   // 最小値
    		y_scale2[1] = 1.0;   // 最大値
    		y_scale2[2] = 0.2;   // 刻み幅
    							// データ
    		data2_y = new double [1][n];
    		for (int i1 = 0; i1 < n; i1++)
    			data2_y[0][i1] = (y2.get(i1)).doubleValue();  
    							// 作図
    		LineGraph gp2 = new LineGraph(title2, g_title2, x_scale, x_p, y_scale2, y_p2, data_x, data2_y, true, false);
    	}
    }
    
    /****************/
    /* 関数値の計算 */
    /****************/
    class Kansu {
    	private int sw;
    					// コンストラクタ
    	Kansu (int s) {sw = s;}
    					// double型関数
    	double snx(double x)
    	{
    
    		double y = 0.0, w[] = new double [1];
    
    		switch (sw) {
    						// 関数値(f(x))の計算(正規分布)
    			case 0:
    				y = 1.0 - Test.p - App.normal(x, w);
    				break;
    						// 関数値(f(x))の計算(F分布)
    			case 1:
    				y = App.f(x, w, Test.dof1, Test.dof2) - 1.0 + Test.p;
    				break;
    						// 関数の微分の計算(F分布)
    			case 2:
    				y = App.f(x, w, Test.dof1, Test.dof2);
    				y = w[0];
    				break;
    		}
    
    		return y;
    	}
    }
    
    /************************/
    /* 科学技術系算用の手法 */
    /************************/
    class App {
    
    	/****************************************/
    	/* f分布の計算(P(X = ff), P(X < ff)) */
    	/*      dd : P(X = ff)                  */
    	/*      df1,df2 : 自由度                */
    	/*      return : P(X < ff)              */
    	/****************************************/
    	static 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 : 収束しなかった     */
    	/****************************************/
    	static 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 * (Test.dof1 - sw);
    			df2 = 0.5 * Test.dof2;
    			a   = 2.0 / (9.0 * (Test.dof1 - sw));
    			a1  = 1.0 - a;
    			b   = 2.0 / (9.0 * Test.dof2);
    			b1  = 1.0 - b;
    
    			yq  = p_normal(ind);
    
    			e   = b1 * b1 - b * yq * yq;
    
    			if (e > 0.8 || (Test.dof1+Test.dof2-sw) <= MAX)
    				sw = -1;
    			else {
    				sw += 1;
    				if ((Test.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)Test.dof2, df2-1.0);
    				y2 = Math.pow((double)Test.dof1, df2);
    				x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
                         2.0 * y1 / y2 / Test.p;
    				f0 = Math.pow(x, 2.0/Test.dof2);
    			}
    	/*
    	     ニュートン法
    	*/
    			Kansu kn1 = new Kansu(1);
    			Kansu kn2 = new Kansu(2);
    
    			ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind, kn1, kn2);
    		}
    
    		return ff;
    	}
    
    	/****************************************/
    	/* Γ(x)の計算(ガンマ関数,近似式) */
    	/*      ier : =0 : normal               */
    	/*            =-1 : x=-n (n=0,1,2,・・・)  */
    	/*      return : 結果                   */
    	/****************************************/
    	static 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)                        */
    	/*************************************************/
    	static 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 : 収束しなかった                               */
    	/******************************************************************/
    	static double p_normal(int ind[])
    	{
    		double u;
    		int sw[] = new int [1];
    
    		Kansu kn = new Kansu(0);
    
    		u        = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, sw, kn);
    		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 : 実際の試行回数                         */
    	/*            (負の時は解を得ることができなかった) */
    	/*      kn1 : 関数を計算するクラスオブジェクト       */
    	/*      kn2 : 関数の微分を計算するクラスオブジェクト */
    	/*      return : 解                                  */
    	/*****************************************************/
    	static double newton(double x1, double eps1, double eps2, int max,
                             int ind[], Kansu kn1, Kansu kn2)
    	{
    		double g, dg, x;
    		int sw;
    
    		x      = x1;
    		ind[0] = 0;
    		sw     = 0;
    
    		while (sw == 0 && ind[0] >= 0) {
    
    			ind[0]++;
    			sw = 1;
    			g  = kn1.snx(x1);
    
    			if (Math.abs(g) > eps2) {
    				if (ind[0] <= max) {
    					dg = kn2.snx(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 : 実際の試行回数                             */
    	/*            (負の時は解を得ることができなかった)     */
    	/*      kn : 関数値を計算するクラスオブジェクト          */
    	/*      return : 解                                      */
    	/*********************************************************/
    	static double bisection(double x1, double x2, double eps1, double eps2, int max, int ind[], Kansu kn)
    	{
    		double f0, f1, f2, x0 = 0.0;
    		int sw;
    
    		f1 = kn.snx(x1);
    		f2 = kn.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      = kn.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;
    	}
    }
    
    /****************************/
    /* 折れ線グラフの描画       */
    /*      coded by Y.Suganuma */
    /****************************/
    class LineGraph extends JFrame {
    
    	Draw_line pn;
    
    	/*********************************************************/
    	/* コンストラクタ(折れ線グラフ2)                      */
    	/*      title_i : グラフ,x軸,及び,y軸のタイトル     */
    	/*      g_title_i : 凡例                                 */
    	/*      x_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_x_i : 小数点以下の桁数(x軸)             */
    	/*      y_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_y_i : 小数点以下の桁数(y軸)             */
    	/*      data_x_i : グラフのデータ(x軸)                */
    	/*      data_y_i : グラフのデータ(y軸)                */
    	/*      d_t_i : タイトル表示の有無                       */
    	/*      d_g_i : 凡例表示の有無                           */
    	/*********************************************************/
    	LineGraph(String title_i[], String g_title_i[], double x_scale_i[],
                  int place_x_i, double y_scale_i[], int place_y_i,
                  double data_x_i[][], double data_y_i[][], boolean d_t_i,
                  boolean d_g_i)
    	{
    					// JFrameクラスのコンストラクタの呼び出し
    		super("折れ線グラフ(2)");
    					// Windowサイズと表示位置を設定
    		int width = 900, height = 600;   // Windowの大きさ(初期サイズ)
    		setSize(width, height);
    		Toolkit tool = getToolkit();
    		Dimension d  = tool.getScreenSize();
    		setLocation(d.width / 2 - width / 2, d.height / 2 - height / 2);
    					// 描画パネル
    		Container cp = getContentPane();
    		pn = new Draw_line(title_i, g_title_i, x_scale_i, place_x_i, y_scale_i, place_y_i, data_x_i, data_y_i, d_t_i, d_g_i, this);
    		cp.add(pn);
    					// ウィンドウを表示
    		setVisible(true);
    					// イベントアダプタ
    		addWindowListener(new WinEnd());
    		addComponentListener(new ComponentResize());
    	}
    
    	/**********************/
    	/* Windowのサイズ変化 */
    	/**********************/
    	class ComponentResize extends ComponentAdapter
    	{
    		public void componentResized(ComponentEvent e)
    		{
    			pn.repaint();
    		}
    	}
    
    	/************/
    	/* 終了処理 */
    	/************/
    	class WinEnd extends WindowAdapter
    	{
    		public void windowClosing(WindowEvent e) {
    			setVisible(false);
    		}
    	}
    }
    
    class Draw_line extends JPanel {
    
    	String title[];   // グラフのタイトル
    	String g_title[];   // 凡例(グラフの内容)
    	String x_title[];   // x軸への表示
    	double x_scale[];   // y軸目盛り
    	double y_scale[];   // y軸目盛り
    	double data_x[][], data_y[][];   // データ
    	boolean d_t;   // タイトル表示の有無
    	boolean d_g;   // 凡例表示の有無
    	boolean type = true;   // 横軸が項目かデータか
    	boolean ver = true;   // 縦か横か
    	int place_x;   // 小数点以下の桁数(x軸)
    	int place_y;   // 小数点以下の桁数(y軸)
    	int width = 900, height = 600;   // Windowの大きさ(初期サイズ)
    	int bx1, bx2, by1, by2;   // 表示切り替えボタンの位置
    	LineGraph line;
    	String change = "横 色";   // 表示切り替えボタン
    	float line_w = 1.0f;   // 折れ線グラフ等の線の太さ
    	boolean line_m = true;   // 折れ線グラフ等にマークを付けるか否か
    	Color cl[] = {Color.black, Color.magenta, Color.blue, Color.orange, Color.cyan,
    	              Color.pink, Color.green, Color.yellow, Color.darkGray, Color.red};   // グラフの色
    	int n_g;   // グラフの数
    
    	/*********************************************************/
    	/* コンストラクタ(折れ線グラフ2)                      */
    	/*      title_i : グラフ,x軸,及び,y軸のタイトル     */
    	/*      g_title_i : 凡例                                 */
    	/*      x_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_x_i : 小数点以下の桁数(x軸)             */
    	/*      y_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_y_i : 小数点以下の桁数(y軸)             */
    	/*      data_x_i : グラフのデータ(x軸)                */
    	/*      data_y_i : グラフのデータ(y軸)                */
    	/*      d_t_i : タイトル表示の有無                       */
    	/*      d_g_i : 凡例表示の有無                           */
    	/*********************************************************/
    	Draw_line(String title_i[], String g_title_i[], double x_scale_i[],
                  int place_x_i, double y_scale_i[], int place_y_i,
                  double data_x_i[][], double data_y_i[][], boolean d_t_i,
                  boolean d_g_i, LineGraph line_i)
    	{
    					// 背景色
    		setBackground(Color.white);
    					// テーブルデータの保存
    		title   = title_i;
    		g_title = g_title_i;
    		x_scale = x_scale_i;
    		place_x = place_x_i;
    		y_scale = y_scale_i;
    		place_y = place_y_i;
    		data_x  = data_x_i;
    		data_y  = data_y_i;
    		d_t     = d_t_i;
    		d_g     = d_g_i;
    		type    = false;
    		line    = line_i;
    					// イベントアダプタ
    		addMouseListener(new ClickMouse(this));
    	}
    
    	/********/
    	/* 描画 */
    	/********/
    	public void paintComponent (Graphics g)
    	{
    		super.paintComponent(g);   // 親クラスの描画(必ず必要)
    
    		double r, x1, y1, sp;
    		int i1, i2, cr, k, k_x, k_y, k1, k2, kx, kx1, ky, ky1, han, len;
    		int x_l, x_r, y_u, y_d;   // 描画領域
    		int f_size;   // フォントサイズ
    		int n_p;   // データの数
    		String s1;
    		Font f;
    		FontMetrics fm;
    		Graphics2D g2 = (Graphics2D)g;
    					//
    					// Windowサイズの取得
    					//
    		Insets insets = line.getInsets();
    		Dimension d = line.getSize();
    		width  = d.width - (insets.left + insets.right);
    		height = d.height - (insets.top + insets.bottom);
    		x_l    = insets.left + 10;
    		x_r    = d.width - insets.right - 10;
    		y_u    = 20;
    		y_d    = d.height - insets.bottom - insets.top;
    					//
    					// グラフタイトルの表示
    					//
    		r      = 0.05;   // タイトル領域の割合
    		f_size = ((y_d - y_u) < (x_r - x_l)) ? (int)((y_d - y_u) * r) : (int)((x_r - x_l) * r);
    		if (f_size < 5)
    			f_size = 5;
    		if (d_t) {
    			f = new Font("TimesRoman", Font.BOLD, f_size);
    			g.setFont(f);
    			fm  = g.getFontMetrics(f);
    			len = fm.stringWidth(title[0]);
    			g.drawString(title[0], (x_l+x_r)/2-len/2, y_d-f_size/2);
    			y_d -= f_size;
    		}
    					//
    					// 表示切り替えボタンの設置
    					//
    		f_size = (int)(0.8 * f_size);
    		if (f_size < 5)
    			f_size = 5;
    		f  = new Font("TimesRoman", Font.PLAIN, f_size);
    		fm = g.getFontMetrics(f);
    		g.setFont(f);
    		g.setColor(Color.yellow);
    		len = fm.stringWidth(change);
    		bx1 = x_r - len - 7 * f_size / 10;
    		by1 = y_u - f_size / 2;
    		bx2 = bx1 + len + f_size / 2;
    		by2 = by1 + 6 * f_size / 5;
    		g.fill3DRect(bx1, by1, len+f_size/2, 6*f_size/5, true);
    		g.setColor(Color.black);
    		g.drawString(change, x_r-len-f_size/2, y_u+f_size/2);
    					//
    					// 凡例の表示
    					//
    		n_g = g_title.length;
    		if (d_g) {
    			han = 0;
    			for (i1 = 0; i1 < n_g; i1++) {
    				len = fm.stringWidth(g_title[i1]);
    				if (len > han)
    					han = len;
    			}
    			han += 15;
    			r    = 0.2;   // 凡例領域の割合
    			k1   = (int)((x_r - x_l) * r);
    			if (han > k1)
    				han = k1;
    			kx = x_r - han;
    			ky = y_u + 3 * f_size / 2;
    			k  = 0;
    			g2.setStroke(new BasicStroke(7.0f));
    			for (i1 = 0; i1 < n_g; i1++) {
    				g.setColor(cl[k]);
    				g.drawLine(kx, ky, kx+10, ky);
    				g.setColor(Color.black);
    				g.drawString(g_title[i1], kx+15, ky+2*f_size/5);
    				k++;
    				if (k >= cl.length)
    					k = 0;
    				ky += f_size;
    			}
    			g2.setStroke(new BasicStroke(1.0f));
    			x_r -= (han + 10);
    		}
    		else
    			x_r -= (int)(0.03 * (x_r - x_l));
    					//
    					// x軸及びy軸のタイトルの表示
    					//
    		if (ver) {   // 縦
    			if (title[1].length() > 0 && !title[1].equals("-")) {
    				len = fm.stringWidth(title[1]);
    				g.drawString(title[1], (x_l+x_r)/2-len/2, y_d-4*f_size/5);
    				y_d -= 7 * f_size / 4;
    			}
    			else
    				y_d -= f_size / 2;
    			if (title[2].length() > 0 && !title[2].equals("-")) {
    				g.drawString(title[2], x_l, y_u+f_size/2);
    				y_u += f_size;
    			}
    		}
    		else {   // 横
    			if (title[2].length() > 0 && !title[2].equals("-")) {
    				len = fm.stringWidth(title[2]);
    				g.drawString(title[2], (x_l+x_r)/2-len/2, y_d-4*f_size/5);
    				y_d -= 7 * f_size / 4;
    			}
    			else
    				y_d -= f_size / 2;
    			if (title[1].length() > 0 && !title[1].equals("-")) {
    				g.drawString(title[1], x_l, y_u+f_size/2);
    				y_u += f_size;
    			}
    		}
    					//
    					// x軸,y軸,及び,各軸の目盛り
    					//
    		f_size = (int)(0.8 * f_size);
    		if (f_size < 5)
    			f_size = 5;
    		f    = new Font("TimesRoman", Font.PLAIN, f_size);
    		fm   = g.getFontMetrics(f);
    		y_d -= 3 * f_size / 2;
    		k_y  = (int)((y_scale[1] - y_scale[0]) / (0.99 * y_scale[2]));
    		k_x  = 0;
    		if (!type)
    			k_x = (int)((x_scale[1] - x_scale[0]) / (0.99 * x_scale[2]));
    		g.setFont(f);
    
    		DecimalFormat df_x, df_y;
    		df_x = new DecimalFormat("#");
    		df_y = new DecimalFormat("#");
    		if (!type) {
    			if (place_x != 0) {
    				s1 = "#.";
    				for (i1 = 0; i1 < place_x; i1++)
    					s1 += "0";
    				df_x = new DecimalFormat(s1);
    			}
    		}
    		if (place_y != 0) {
    			s1 = "#.";
    			for (i1 = 0; i1 < place_y; i1++)
    				s1 += "0";
    			df_y = new DecimalFormat(s1);
    		}
    						// 縦表示
    		if (ver) {
    							// y軸
    			y1  = y_scale[0];
    			len = 0;
    			for (i1 = 0; i1 < k_y+1; i1++) {
    				s1 = df_y.format(y1);
    				k1 = fm.stringWidth(s1);
    				if (k1 > len)
    					len = k1;
    				y1 += y_scale[2];
    			}
    			g.drawLine(x_l+len+5, y_u, x_l+len+5, y_d);
    			g.drawLine(x_r, y_u, x_r, y_d);
    			y1 = y_scale[0];
    			x1 = y_d;
    			sp = (double)(y_d - y_u) / k_y;
    			for (i1 = 0; i1 < k_y+1; i1++) {
    				ky = (int)Math.round(x1);
    				s1 = df_y.format(y1);
    				k1 = fm.stringWidth(s1);
    				g.drawString(s1, x_l+len-k1, ky+f_size/2);
    				g.drawLine(x_l+len+5, ky, x_r, ky);
    				y1 += y_scale[2];
    				x1 -= sp;
    			}
    			x_l += (len + 5);
    							// x軸
    			if (type) {
    				n_p = x_title.length;
    				sp  = (double)(x_r - x_l) / n_p;
    				x1  = x_l + sp / 2.0;
    				for (i1 = 0; i1 < n_p; i1++) {
    					kx = (int)Math.round(x1);
    					k1 = fm.stringWidth(x_title[i1]);
    					g.drawString(x_title[i1], kx-k1/2, y_d+6*f_size/5);
    					g.drawLine(kx, y_d, kx, y_d-5);
    					x1 += sp;
    				}
    			}
    			else {
    				x1 = x_scale[0];
    				y1 = x_l;
    				sp = (double)(x_r - x_l) / k_x;
    				for (i1 = 0; i1 < k_x+1; i1++) {
    					kx = (int)Math.round(y1);
    					s1 = df_x.format(x1);
    					k1 = fm.stringWidth(s1);
    					g.drawString(s1, kx-k1/2, y_d+6*f_size/5);
    					g.drawLine(kx, y_d, kx, y_u);
    					x1 += x_scale[2];
    					y1 += sp;
    				}
    			}
    		}
    						// 横表示
    		else {
    							// y軸
    			if (type) {
    				n_p = x_title.length;
    				len = 0;
    				for (i1 = 0; i1 < n_p; i1++) {
    					k1 = fm.stringWidth(x_title[i1]);
    					if (k1 > len)
    						len = k1;
    				}
    				g.drawLine(x_l+len+5, y_u, x_l+len+5, y_d);
    				g.drawLine(x_r, y_u, x_r, y_d);
    				sp = (double)(y_d - y_u) / n_p;
    				x1 = y_d - sp / 2.0;
    				for (i1 = 0; i1 < n_p; i1++) {
    					ky = (int)Math.round(x1);
    					k1 = fm.stringWidth(x_title[n_p-1-i1]);
    					g.drawString(x_title[n_p-1-i1], x_l+len-k1, ky+f_size/2);
    					g.drawLine(x_l+len+5, ky, x_l+len+10, ky);
    					x1 -= sp;
    				}
    				g.drawLine(x_l+len+5, y_u, x_r, y_u);
    				g.drawLine(x_l+len+5, y_d, x_r, y_d);
    				x_l += (len + 5);
    			}
    			else {
    				y1  = x_scale[0];
    				len = 0;
    				for (i1 = 0; i1 < k_x+1; i1++) {
    					s1 = df_x.format(y1);
    					k1 = fm.stringWidth(s1);
    					if (k1 > len)
    						len = k1;
    					y1 += x_scale[2];
    				}
    				g.drawLine(x_l+len+5, y_u, x_l+len+5, y_d);
    				g.drawLine(x_r, y_u, x_r, y_d);
    				y1 = x_scale[0];
    				x1 = y_d;
    				sp = (double)(y_d - y_u) / k_x;
    				for (i1 = 0; i1 < k_x+1; i1++) {
    					ky = (int)Math.round(x1);
    					s1 = df_x.format(y1);
    					k1 = fm.stringWidth(s1);
    					g.drawString(s1, x_l+len-k1, ky+f_size/2);
    					g.drawLine(x_l+len+5, ky, x_r, ky);
    					y1 += x_scale[2];
    					x1 -= sp;
    				}
    				x_l += (len + 5);
    			}
    							// x軸
    			x1 = y_scale[0];
    			y1 = x_l;
    			sp = (double)(x_r - x_l) / k_y;
    			for (i1 = 0; i1 < k_y+1; i1++) {
    				kx = (int)Math.round(y1);
    				s1 = df_y.format(x1);
    				k1 = fm.stringWidth(s1);
    				g.drawString(s1, kx-k1/2, y_d+6*f_size/5);
    				g.drawLine(kx, y_d, kx, y_u);
    				x1 += y_scale[2];
    				y1 += sp;
    			}
    		}
    					//
    					// グラフの表示
    					//
    		g2.setStroke(new BasicStroke(line_w));
    		cr = (int)line_w + 6;
    						// 縦表示
    		if (ver) {
    			if (type) {
    				n_p = x_title.length;
    				sp  = (double)(x_r - x_l) / n_p;
    				k1  = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					x1  = x_l + sp / 2.0;
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						kx = (int)Math.round(x1);
    						ky = y_d - (int)((y_d - y_u) * (data_y[i1][i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1  = kx;
    						ky1  = ky;
    						x1  += sp;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    			else {
    				n_p = data_x[0].length;
    				k1  = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						kx = x_l + (int)((x_r - x_l) * (data_x[i1][i2] - x_scale[0]) / (x_scale[1] - x_scale[0]));
    						ky = y_d - (int)((y_d - y_u) * (data_y[i1][i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1 = kx;
    						ky1 = ky;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    		}
    						// 横表示
    		else {
    			if (type) {
    				n_p = x_title.length;
    				sp = (double)(y_d - y_u) / n_p;
    				k1 = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					y1  = y_d - sp / 2.0;
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						ky = (int)Math.round(y1);
    						kx = x_l + (int)((x_r - x_l) * (data_y[i1][n_p-1-i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1  = kx;
    						ky1  = ky;
    						y1  -= sp;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    			else {
    				n_p = data_x[0].length;
    				k1  = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						kx = x_l + (int)((x_r - x_l) * (data_y[i1][i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						ky = y_d - (int)((y_d - y_u) * (data_x[i1][i2] - x_scale[0]) / (x_scale[1] - x_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1 = kx;
    						ky1 = ky;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    		}
    		g2.setStroke(new BasicStroke(1.0f));
    	}
    
    	/************************************/
    	/* マウスがクリックされたときの処理 */
    	/************************************/
    	class ClickMouse extends MouseAdapter
    	{
    		Draw_line dr;
    
    		ClickMouse(Draw_line dr1)
    		{
    			dr = dr1;
    		}
    
    		public void mouseClicked(MouseEvent e)
    		{
    			int xp = e.getX();
    			int yp = e.getY();
    					// 縦表示と横表示の変換
    			if (xp > bx1 && xp < bx1+(bx2-bx1)/2 && yp > by1 && yp < by2) {
    				if (ver) {
    					ver    = false;
    					change = "縦 色";
    				}
    				else {
    					ver    = true;
    					change = "横 色";
    				}
    				repaint();
    			}
    					// グラフの色,線の太さ等
    			else if (xp > bx1+(bx2-bx1)/2 && xp < bx2 && yp > by1 && yp < by2) {
    				Modify md = new Modify(dr.line, dr);
    				md.setVisible(true);
    			}
    		}
    	}
    }
    
    /****************************/
    /* 色及び線の太さの変更     */
    /*      coded by Y.Suganuma */
    /****************************/
    class Modify extends JDialog implements ActionListener, TextListener {
    	Draw_line dr;   // 折れ線グラフ
    	JButton bt_dr;
    	TextField rgb[];
    	TextField r[];
    	TextField g[];
    	TextField b[];
    	JTextField tx;
    	JRadioButton r1, r2;
    	float line_w = 1.0f;   // 折れ線グラフ等の線の太さ
    	boolean line_m = true;   // 折れ線グラフ等にマークを付けるか否か
    	Color cl[];   // グラフの色
    	int n_g;   // グラフの数
    	int wd;   // 線の太さを変更するか
    	int mk;   // マークを変更するか
    	int n;
    	JPanel jp[];
    					// 折れ線グラフ
    	Modify(Frame host, Draw_line dr1)
    	{
    		super(host, "色と線の変更", true);
    							// 初期設定
    		dr  = dr1;
    		wd  = 1;
    		mk  = 1;
    		n_g = dr.n_g;
    		if (n_g > 10)
    			n_g = 10;
    		n      = n_g + 3;
    		line_w = dr.line_w;
    		line_m = dr.line_m;
    		cl     = new Color[n_g];
    		for (int i1 = 0; i1 < n_g; i1++)
    			cl[i1] = dr.cl[i1];
    		set();
    							// ボタン
    		Font f = new Font("TimesRoman", Font.BOLD, 20);
    		bt_dr = new JButton("OK");
    		bt_dr.setFont(f);
    		bt_dr.addActionListener(this);
    		jp[n-1].add(bt_dr);
    	}
    					// 設定
    	void set()
    	{
    		setSize(450, 60*(n));
    		Container cp = getContentPane();
    		cp.setBackground(Color.white);
    		cp.setLayout(new GridLayout(n, 1, 5, 5));
    		jp = new JPanel[n];
    		for (int i1 = 0; i1 < n; i1++) {
    			jp[i1] = new JPanel();
    			cp.add(jp[i1]);
    		}
    		Font f = new Font("TimesRoman", Font.BOLD, 20);
    							// 色の変更
    		JLabel lb[][] = new JLabel[n_g][3];
    		rgb = new TextField[n_g];
    		r = new TextField[n_g];
    		g = new TextField[n_g];
    		b = new TextField[n_g];
    		for (int i1 = 0; i1 < n_g; i1++) {
    			rgb[i1] = new TextField(3);
    			rgb[i1].setFont(f);
    			rgb[i1].setBackground(new Color(cl[i1].getRed(), cl[i1].getGreen(), cl[i1].getBlue()));
    			jp[i1].add(rgb[i1]);
    			lb[i1][0] = new JLabel(" 赤");
    			lb[i1][0].setFont(f);
    			jp[i1].add(lb[i1][0]);
    			r[i1] = new TextField(3);
    			r[i1].setFont(f);
    			r[i1].setBackground(Color.white);
    			r[i1].setText(Integer.toString(cl[i1].getRed()));
    			r[i1].addTextListener(this);
    			jp[i1].add(r[i1]);
    			lb[i1][1] = new JLabel("緑");
    			lb[i1][1].setFont(f);
    			jp[i1].add(lb[i1][1]);
    			g[i1] = new TextField(3);
    			g[i1].setFont(f);
    			g[i1].setBackground(Color.white);
    			g[i1].setText(Integer.toString(cl[i1].getGreen()));
    			g[i1].addTextListener(this);
    			jp[i1].add(g[i1]);
    			lb[i1][2] = new JLabel("青");
    			lb[i1][2].setFont(f);
    			jp[i1].add(lb[i1][2]);
    			b[i1] = new TextField(3);
    			b[i1].setFont(f);
    			b[i1].setBackground(Color.white);
    			b[i1].setText(Integer.toString(cl[i1].getBlue()));
    			b[i1].addTextListener(this);
    			jp[i1].add(b[i1]);
    		}
    							// 線の変更
    		if (wd > 0) {
    			JLabel lb1 = new JLabel("線の太さ:");
    			lb1.setFont(f);
    			jp[n_g].add(lb1);
    			tx = new JTextField(2);
    			tx.setFont(f);
    			tx.setBackground(Color.white);
    			tx.setText(Integer.toString((int)line_w));
    			jp[n_g].add(tx);
    		}
    
    		if (mk > 0) {
    			JLabel lb2 = new JLabel("マーク:");
    			lb2.setFont(f);
    			jp[n-2].add(lb2);
    			ButtonGroup gp = new ButtonGroup();
    			r1 = new JRadioButton("付ける");
    			r1.setFont(f);
    			gp.add(r1);
    			jp[n-2].add(r1);
    			r2 = new JRadioButton("付けない");
    			r2.setFont(f);
    			gp.add(r2);
    			jp[n-2].add(r2);
    			if (line_m)
    				r1.doClick();
    			else
    				r2.doClick();
    		}
    	}
    					// TextFieldの内容が変更されたときの処理
    	public void textValueChanged(TextEvent e)
    	{
    		for (int i1 = 0; i1 < n_g; i1++) {
    			if (e.getSource() == r[i1] || e.getSource() == g[i1] || e.getSource() == b[i1]) {
    				String str = r[i1].getText();
    				int rc = str.length()>0 ? Integer.parseInt(str) : 0;
    				str = g[i1].getText();
    				int gc = str.length()>0 ? Integer.parseInt(str) : 0;
    				str = b[i1].getText();
    				int bc = str.length()>0 ? Integer.parseInt(str) : 0;
    				rgb[i1].setBackground(new Color(rc, gc, bc));
    			}
    		}
    	}
    					// 値の設定
    	public void actionPerformed(ActionEvent e)
    	{
    		for (int i1 = 0; i1 < n_g; i1++) {
    			String str = r[i1].getText();
    			int rc = str.length()>0 ? Integer.parseInt(str) : 0;
    			str = g[i1].getText();
    			int gc = str.length()>0 ? Integer.parseInt(str) : 0;
    			str = b[i1].getText();
    			int bc = str.length()>0 ? Integer.parseInt(str) : 0;
    			dr.cl[i1] = new Color(rc, gc, bc);
    		}
    		dr.line_w = Integer.parseInt(tx.getText());
    		if (r1.isSelected())
    			dr.line_m = true;
    		else
    			dr.line_m = false;
    		dr.repaint();
    
    		setVisible(false);
    	}
    }
    			
    アプレット版
    /****************************/
    /* F分布                   */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.io.*;
    import java.text.*;
    import java.awt.*;
    import javax.swing.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.*;
    
    public class F extends Applet implements ActionListener, ItemListener {
    	double p;   // α%値を計算するとき時α/100を設定
    	int dof1, dof2;   // 自由度
    	int sw = 0;
    	Choice ch;
    	Label lb1, lb2, lb3;
    	TextField tx1, tx2, tx3, tx4, tx5;
    	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("自由度1:"));
    		tx1 = new TextField("6", 3);
    		pn1.add(tx1);
    
    		pn1.add(new Label("自由度2:"));
    		tx2 = new TextField("4", 3);
    		pn1.add(tx2);
    
    		pn1.add(new Label("目的:"));
    		ch = new Choice ();
    		ch.addItem("確率(1点)");
    		ch.addItem("確率(複数点)");
    		ch.addItem("p%値");
    		pn1.add(ch);
    		ch.addItemListener(this);
    					// 中央のパネル(データ)
    		Panel pn2 = new Panel();
    		add(pn2, BorderLayout.CENTER);
    
    		lb1 = new Label("x:");
    		pn2.add(lb1);
    		tx3 = new TextField("1.0", 3);
    		pn2.add(tx3);
    
    		lb2 = new Label("上限:");
    		lb2.setVisible(false);
    		pn2.add(lb2);
    		tx4 = new TextField("10.0", 3);
    		tx4.setVisible(false);
    		pn2.add(tx4);
    
    		lb3 = new Label("刻み幅:");
    		lb3.setVisible(false);
    		pn2.add(lb3);
    		tx5 = new TextField("0.1", 3);
    		tx5.setVisible(false);
    		pn2.add(tx5);
    
    		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);
    
    		ta1 = new TextArea(25, 25);
    		pn3.add(ta1);
    
    		ta2 = new TextArea(25, 25);
    		ta2.setVisible(false);
    		pn3.add(ta2);
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    			double z[] = new double [1];
    			int sw1[] = new int [1];
    					// 問題の設定
    			dof1 = Integer.parseInt(tx1.getText());
    			dof2 = Integer.parseInt(tx2.getText());
    					// 計算
    							// 確率(1点)
    			if (sw == 0) {
    				double x = Double.parseDouble(tx3.getText());
    				double y = f(x, z, dof1, dof2);
    				ta1.setText("P(X = " + String.format("%.5f",x) + ") = " + String.format("%.5f",z[0]) + "\n");
    				ta1.append("P(X < " + String.format("%.5f",x) + ") = " + String.format("%.5f",y));
    			}
    							// 確率(複数点)
    			else if (sw == 1) {
    				double from = Double.parseDouble(tx3.getText());
    				double to = Double.parseDouble(tx4.getText());
    				double h = Double.parseDouble(tx5.getText());
    									// データ取得
    				ArrayList <Double> x1 = new ArrayList <Double> ();
    				ArrayList <Double> y1 = new ArrayList <Double> ();
    				ArrayList <Double> y2 = new ArrayList <Double> ();
    				ta1.setText("密度関数\n");
    				ta2.setText("分布関数\n");
    				for (double x = from; x < to+0.5*h; x += h) {
    					double y = f(x, z, dof1, dof2);
    					ta1.append(String.format("%.5f",x) + " " + String.format("%.5f",z[0]) + "\n");
    					ta2.append(String.format("%.5f",x) + " " + String.format("%.5f",y) + "\n");
    					x1.add(x);
    					y1.add(z[0]);
    					y2.add(y);
    				}
    									// グラフの描画
    				graph(x1, y1, y2);
    			}
    							// p%値
    			else {
    				double x = Double.parseDouble(tx3.getText());
    				p = x / 100.0;
    				if (p < 1.0e-7)
    					ta1.setText(String.format("%.5f",x) + "%値 = ∞");
    				else if ((1.0-p) < 1.0e-7)
    					ta1.setText(String.format("%.5f",x) + "%値 = 0");
    				else {
    					double y = p_f(sw1);
    					if (sw1[0] < 0)
    						ta1.setText("収束しませんでした");
    					else
    						ta1.setText(String.format("%.5f",x) + "%値 = " + String.format("%.5f",y));
    				}
    			}
    		}
    	}
    
    	/************************/
    	/* 選択されたときの処理 */
    	/************************/
    	public void itemStateChanged(ItemEvent e)
    	{
    		if (e.getItemSelectable() == ch) {
    			sw = ch.getSelectedIndex();
    					// 確率(1点)
    			if (sw == 0) {
    				lb1.setText("x:");
    				tx3.setText("1.0");
    				tx3.setEditable(true);
    				lb2.setVisible(false);
    				tx4.setVisible(false);
    				lb3.setVisible(false);
    				tx5.setVisible(false);
    				ta2.setVisible(false);
    			}
    					// 確率(複数点)
    			else if (sw == 1) {
    				lb1.setText("下限:");
    				tx3.setText("0.0");
    				tx3.setEditable(false);
    				lb2.setVisible(true);
    				tx4.setVisible(true);
    				lb3.setVisible(true);
    				tx5.setVisible(true);
    				ta2.setVisible(true);
    			}
    					// p%値
    			else {
    				lb1.setText("%:");
    				tx3.setText("2.5");
    				tx3.setEditable(true);
    				lb2.setVisible(false);
    				tx4.setVisible(false);
    				lb3.setVisible(false);
    				tx5.setVisible(false);
    				ta2.setVisible(false);
    			}
    			getParent().validate();
    		}
    	}
    
    	/****************/
    	/* 関数値の計算 */
    	/****************/
    	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))の計算(χ2分布)
    			case 1:
    				y = f(x, w, dof1, dof2) - 1.0 + p;
    				break;
    						// 関数の微分の計算(χ2分布)
    			case 2:
    				y = f(x, w, dof1, dof2);
    				y = w[0];
    				break;
    		}
    
    		return y;
    	}
    
    	/***************************************/
    	/*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-10, 1.0e-15, 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;
    	}
    
    	/*************************/
    	/* グラフの描画          */
    	/*      x1 : x座標データ */
    	/*      y1 : 密度関数    */
    	/*      y2 : 分布関数    */
    	/*************************/
    	void graph(ArrayList <Double> x1, ArrayList <Double> y1, ArrayList <Double> y2)
    	{
    					// 密度関数
    		String title1[];   // グラフ,x軸,及び,y軸のタイトル
    		String g_title1[];   // 凡例(グラフの内容)
    		double x_scale[];   // x軸目盛り
    		double y_scale1[];   // y軸目盛り
    		double data_x[][], data1_y[][];   // データ
    		int n = x1.size();
    							// グラフ,x軸,及び,y軸のタイトル
    		title1 = new String [3];
    		title1[0] = "密度関数(F分布 自由度 = " + dof1 + ", " + dof2 + ")";
    		title1[1] = "x";
    		title1[2] = "f(x)";
    							// 凡例
    		g_title1 = new String [1];
    		g_title1[0] = "密度関数";
    							// x軸目盛り
    		x_scale = new double[3];
    		double x_max = (x1.get(n-1)).doubleValue();
    		double x_min = 0.0;
    		double x_step = (x_max - x_min) / 10;
    		int x_p = 0;
    		boolean ok = true;
    		int k = 0;
    		while (ok) {
    			if (x_step < 1.0) {
    				x_step *= 10;
    				k++;
    			}
    			else if (x_step >= 10.0) {
    				x_step /= 10;
    				k--;
    			}
    			else {
    				ok = false;
    				if (x_step-(int)x_step > 1.0e-5)
    					x_step = (int)x_step + 1;
    				else
    					x_step = (int)x_step;
    				if (k != 0) {
    					x_step = x_step * Math.pow(10, -k);
    					if (k > 0)
    						x_p = k;
    				}
    				double t = x_min;
    				while (t < x_max-0.001*x_step)
    					t += x_step;
    				x_max = t;
    			}
    		}
    		x_scale[0] = x_min;   // 最小値
    		x_scale[1] = x_max;   // 最大値
    		x_scale[2] = x_step;   // 刻み幅
    							// y軸目盛り
    		y_scale1 = new double[3];
    		double y_max = 0.0;
    		for (int i1 = 0; i1 < n; i1++) {
    			if ((y1.get(i1)).doubleValue() > y_max)
    				y_max = (y1.get(i1)).doubleValue();
    		}
    		double y_step = y_max / 5;
    		int y_p1 = 0;
    		ok = true;
    		k  = 0;
    		while (ok) {
    			if (y_step < 1.0) {
    				y_step *= 10;
    				k++;
    			}
    			else if (y_step >= 10.0) {
    				y_step /= 10;
    				k--;
    			}
    			else {
    				ok = false;
    				if (y_step-(int)y_step > 1.0e-5)
    					y_step = (int)y_step + 1;
    				else
    					y_step = (int)y_step;
    				if (k != 0) {
    					y_step = y_step * Math.pow(10, -k);
    					if (k > 0)
    						y_p1 = k;
    				}
    				double t = 0.0;
    				while (t < y_max-0.001*y_step)
    					t += y_step;
    				y_max = t;
    			}
    		}
    		y_scale1[0] = 0.0;   // 最小値
    		y_scale1[1] = y_max;   // 最大値
    		y_scale1[2] = y_step;   // 刻み幅
    							// データ
    		data_x = new double [1][n];
    		data1_y = new double [1][n];
    		for (int i1 = 0; i1 < n; i1++)
    			data_x[0][i1] = (x1.get(i1)).doubleValue();  
    		for (int i1 = 0; i1 < n; i1++)
    			data1_y[0][i1] = (y1.get(i1)).doubleValue();  
    							// 作図
    		LineGraph gp1 = new LineGraph(title1, g_title1, x_scale, x_p, y_scale1, y_p1, data_x, data1_y, true, false);
    					// 分布関数
    		String title2[];   // グラフ,x軸,及び,y軸のタイトル
    		String g_title2[];   // 凡例(グラフの内容)
    		double y_scale2[];   // y軸目盛り
    		double data2_y[][];   // データ
    							// グラフ,x軸,及び,y軸のタイトル
    		title2 = new String [3];
    		title2[0] = "分布関数(F分布 自由度 = " + dof1 + ", " + dof2 + ")";
    		title2[1] = "x";
    		title2[2] = "F(x)";
    							// 凡例
    		g_title2 = new String [1];
    		g_title2[0] = "分布関数";
    							// y軸目盛り
    		y_scale2 = new double[3];
    		int y_p2 = 1;
    		y_scale2[0] = 0.0;   // 最小値
    		y_scale2[1] = 1.0;   // 最大値
    		y_scale2[2] = 0.2;   // 刻み幅
    							// データ
    		data2_y = new double [1][n];
    		for (int i1 = 0; i1 < n; i1++)
    			data2_y[0][i1] = (y2.get(i1)).doubleValue();  
    							// 作図
    		LineGraph gp2 = new LineGraph(title2, g_title2, x_scale, x_p, y_scale2, y_p2, data_x, data2_y, true, false);
    	}
    }
    
    /****************************/
    /* 折れ線グラフの描画       */
    /*      coded by Y.Suganuma */
    /****************************/
    class LineGraph extends JFrame {
    
    	Draw_line pn;
    
    	/*********************************************************/
    	/* コンストラクタ(折れ線グラフ2)                      */
    	/*      title_i : グラフ,x軸,及び,y軸のタイトル     */
    	/*      g_title_i : 凡例                                 */
    	/*      x_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_x_i : 小数点以下の桁数(x軸)             */
    	/*      y_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_y_i : 小数点以下の桁数(y軸)             */
    	/*      data_x_i : グラフのデータ(x軸)                */
    	/*      data_y_i : グラフのデータ(y軸)                */
    	/*      d_t_i : タイトル表示の有無                       */
    	/*      d_g_i : 凡例表示の有無                           */
    	/*********************************************************/
    	LineGraph(String title_i[], String g_title_i[], double x_scale_i[],
                  int place_x_i, double y_scale_i[], int place_y_i,
                  double data_x_i[][], double data_y_i[][], boolean d_t_i,
                  boolean d_g_i)
    	{
    					// JFrameクラスのコンストラクタの呼び出し
    		super("折れ線グラフ(2)");
    					// Windowサイズと表示位置を設定
    		int width = 900, height = 600;   // Windowの大きさ(初期サイズ)
    		setSize(width, height);
    		Toolkit tool = getToolkit();
    		Dimension d  = tool.getScreenSize();
    		setLocation(d.width / 2 - width / 2, d.height / 2 - height / 2);
    					// 描画パネル
    		Container cp = getContentPane();
    		pn = new Draw_line(title_i, g_title_i, x_scale_i, place_x_i, y_scale_i, place_y_i, data_x_i, data_y_i, d_t_i, d_g_i, this);
    		cp.add(pn);
    					// ウィンドウを表示
    		setVisible(true);
    					// イベントアダプタ
    		addWindowListener(new WinEnd());
    		addComponentListener(new ComponentResize());
    	}
    
    	/**********************/
    	/* Windowのサイズ変化 */
    	/**********************/
    	class ComponentResize extends ComponentAdapter
    	{
    		public void componentResized(ComponentEvent e)
    		{
    			pn.repaint();
    		}
    	}
    
    	/************/
    	/* 終了処理 */
    	/************/
    	class WinEnd extends WindowAdapter
    	{
    		public void windowClosing(WindowEvent e) {
    			setVisible(false);
    		}
    	}
    }
    
    class Draw_line extends JPanel {
    
    	String title[];   // グラフのタイトル
    	String g_title[];   // 凡例(グラフの内容)
    	String x_title[];   // x軸への表示
    	double x_scale[];   // y軸目盛り
    	double y_scale[];   // y軸目盛り
    	double data_x[][], data_y[][];   // データ
    	boolean d_t;   // タイトル表示の有無
    	boolean d_g;   // 凡例表示の有無
    	boolean type = true;   // 横軸が項目かデータか
    	boolean ver = true;   // 縦か横か
    	int place_x;   // 小数点以下の桁数(x軸)
    	int place_y;   // 小数点以下の桁数(y軸)
    	int width = 900, height = 600;   // Windowの大きさ(初期サイズ)
    	int bx1, bx2, by1, by2;   // 表示切り替えボタンの位置
    	LineGraph line;
    	String change = "横 色";   // 表示切り替えボタン
    	float line_w = 1.0f;   // 折れ線グラフ等の線の太さ
    	boolean line_m = true;   // 折れ線グラフ等にマークを付けるか否か
    	Color cl[] = {Color.black, Color.magenta, Color.blue, Color.orange, Color.cyan,
    	              Color.pink, Color.green, Color.yellow, Color.darkGray, Color.red};   // グラフの色
    	int n_g;   // グラフの数
    
    	/*********************************************************/
    	/* コンストラクタ(折れ線グラフ2)                      */
    	/*      title_i : グラフ,x軸,及び,y軸のタイトル     */
    	/*      g_title_i : 凡例                                 */
    	/*      x_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_x_i : 小数点以下の桁数(x軸)             */
    	/*      y_scale_i : データの最小値,最大値,目盛幅(y) */
    	/*      place_y_i : 小数点以下の桁数(y軸)             */
    	/*      data_x_i : グラフのデータ(x軸)                */
    	/*      data_y_i : グラフのデータ(y軸)                */
    	/*      d_t_i : タイトル表示の有無                       */
    	/*      d_g_i : 凡例表示の有無                           */
    	/*********************************************************/
    	Draw_line(String title_i[], String g_title_i[], double x_scale_i[],
                  int place_x_i, double y_scale_i[], int place_y_i,
                  double data_x_i[][], double data_y_i[][], boolean d_t_i,
                  boolean d_g_i, LineGraph line_i)
    	{
    					// 背景色
    		setBackground(Color.white);
    					// テーブルデータの保存
    		title   = title_i;
    		g_title = g_title_i;
    		x_scale = x_scale_i;
    		place_x = place_x_i;
    		y_scale = y_scale_i;
    		place_y = place_y_i;
    		data_x  = data_x_i;
    		data_y  = data_y_i;
    		d_t     = d_t_i;
    		d_g     = d_g_i;
    		type    = false;
    		line    = line_i;
    					// イベントアダプタ
    		addMouseListener(new ClickMouse(this));
    	}
    
    	/********/
    	/* 描画 */
    	/********/
    	public void paintComponent (Graphics g)
    	{
    		super.paintComponent(g);   // 親クラスの描画(必ず必要)
    
    		double r, x1, y1, sp;
    		int i1, i2, cr, k, k_x, k_y, k1, k2, kx, kx1, ky, ky1, han, len;
    		int x_l, x_r, y_u, y_d;   // 描画領域
    		int f_size;   // フォントサイズ
    		int n_p;   // データの数
    		String s1;
    		Font f;
    		FontMetrics fm;
    		Graphics2D g2 = (Graphics2D)g;
    					//
    					// Windowサイズの取得
    					//
    		Insets insets = line.getInsets();
    		Dimension d = line.getSize();
    		width  = d.width - (insets.left + insets.right);
    		height = d.height - (insets.top + insets.bottom);
    		x_l    = insets.left + 10;
    		x_r    = d.width - insets.right - 10;
    		y_u    = 20;
    		y_d    = d.height - insets.bottom - insets.top;
    					//
    					// グラフタイトルの表示
    					//
    		r      = 0.05;   // タイトル領域の割合
    		f_size = ((y_d - y_u) < (x_r - x_l)) ? (int)((y_d - y_u) * r) : (int)((x_r - x_l) * r);
    		if (f_size < 5)
    			f_size = 5;
    		if (d_t) {
    			f = new Font("TimesRoman", Font.BOLD, f_size);
    			g.setFont(f);
    			fm  = g.getFontMetrics(f);
    			len = fm.stringWidth(title[0]);
    			g.drawString(title[0], (x_l+x_r)/2-len/2, y_d-f_size/2);
    			y_d -= f_size;
    		}
    					//
    					// 表示切り替えボタンの設置
    					//
    		f_size = (int)(0.8 * f_size);
    		if (f_size < 5)
    			f_size = 5;
    		f  = new Font("TimesRoman", Font.PLAIN, f_size);
    		fm = g.getFontMetrics(f);
    		g.setFont(f);
    		g.setColor(Color.yellow);
    		len = fm.stringWidth(change);
    		bx1 = x_r - len - 7 * f_size / 10;
    		by1 = y_u - f_size / 2;
    		bx2 = bx1 + len + f_size / 2;
    		by2 = by1 + 6 * f_size / 5;
    		g.fill3DRect(bx1, by1, len+f_size/2, 6*f_size/5, true);
    		g.setColor(Color.black);
    		g.drawString(change, x_r-len-f_size/2, y_u+f_size/2);
    					//
    					// 凡例の表示
    					//
    		n_g = g_title.length;
    		if (d_g) {
    			han = 0;
    			for (i1 = 0; i1 < n_g; i1++) {
    				len = fm.stringWidth(g_title[i1]);
    				if (len > han)
    					han = len;
    			}
    			han += 15;
    			r    = 0.2;   // 凡例領域の割合
    			k1   = (int)((x_r - x_l) * r);
    			if (han > k1)
    				han = k1;
    			kx = x_r - han;
    			ky = y_u + 3 * f_size / 2;
    			k  = 0;
    			g2.setStroke(new BasicStroke(7.0f));
    			for (i1 = 0; i1 < n_g; i1++) {
    				g.setColor(cl[k]);
    				g.drawLine(kx, ky, kx+10, ky);
    				g.setColor(Color.black);
    				g.drawString(g_title[i1], kx+15, ky+2*f_size/5);
    				k++;
    				if (k >= cl.length)
    					k = 0;
    				ky += f_size;
    			}
    			g2.setStroke(new BasicStroke(1.0f));
    			x_r -= (han + 10);
    		}
    		else
    			x_r -= (int)(0.03 * (x_r - x_l));
    					//
    					// x軸及びy軸のタイトルの表示
    					//
    		if (ver) {   // 縦
    			if (title[1].length() > 0 && !title[1].equals("-")) {
    				len = fm.stringWidth(title[1]);
    				g.drawString(title[1], (x_l+x_r)/2-len/2, y_d-4*f_size/5);
    				y_d -= 7 * f_size / 4;
    			}
    			else
    				y_d -= f_size / 2;
    			if (title[2].length() > 0 && !title[2].equals("-")) {
    				g.drawString(title[2], x_l, y_u+f_size/2);
    				y_u += f_size;
    			}
    		}
    		else {   // 横
    			if (title[2].length() > 0 && !title[2].equals("-")) {
    				len = fm.stringWidth(title[2]);
    				g.drawString(title[2], (x_l+x_r)/2-len/2, y_d-4*f_size/5);
    				y_d -= 7 * f_size / 4;
    			}
    			else
    				y_d -= f_size / 2;
    			if (title[1].length() > 0 && !title[1].equals("-")) {
    				g.drawString(title[1], x_l, y_u+f_size/2);
    				y_u += f_size;
    			}
    		}
    					//
    					// x軸,y軸,及び,各軸の目盛り
    					//
    		f_size = (int)(0.8 * f_size);
    		if (f_size < 5)
    			f_size = 5;
    		f    = new Font("TimesRoman", Font.PLAIN, f_size);
    		fm   = g.getFontMetrics(f);
    		y_d -= 3 * f_size / 2;
    		k_y  = (int)((y_scale[1] - y_scale[0]) / (0.99 * y_scale[2]));
    		k_x  = 0;
    		if (!type)
    			k_x = (int)((x_scale[1] - x_scale[0]) / (0.99 * x_scale[2]));
    		g.setFont(f);
    
    		DecimalFormat df_x, df_y;
    		df_x = new DecimalFormat("#");
    		df_y = new DecimalFormat("#");
    		if (!type) {
    			if (place_x != 0) {
    				s1 = "#.";
    				for (i1 = 0; i1 < place_x; i1++)
    					s1 += "0";
    				df_x = new DecimalFormat(s1);
    			}
    		}
    		if (place_y != 0) {
    			s1 = "#.";
    			for (i1 = 0; i1 < place_y; i1++)
    				s1 += "0";
    			df_y = new DecimalFormat(s1);
    		}
    						// 縦表示
    		if (ver) {
    							// y軸
    			y1  = y_scale[0];
    			len = 0;
    			for (i1 = 0; i1 < k_y+1; i1++) {
    				s1 = df_y.format(y1);
    				k1 = fm.stringWidth(s1);
    				if (k1 > len)
    					len = k1;
    				y1 += y_scale[2];
    			}
    			g.drawLine(x_l+len+5, y_u, x_l+len+5, y_d);
    			g.drawLine(x_r, y_u, x_r, y_d);
    			y1 = y_scale[0];
    			x1 = y_d;
    			sp = (double)(y_d - y_u) / k_y;
    			for (i1 = 0; i1 < k_y+1; i1++) {
    				ky = (int)Math.round(x1);
    				s1 = df_y.format(y1);
    				k1 = fm.stringWidth(s1);
    				g.drawString(s1, x_l+len-k1, ky+f_size/2);
    				g.drawLine(x_l+len+5, ky, x_r, ky);
    				y1 += y_scale[2];
    				x1 -= sp;
    			}
    			x_l += (len + 5);
    							// x軸
    			if (type) {
    				n_p = x_title.length;
    				sp  = (double)(x_r - x_l) / n_p;
    				x1  = x_l + sp / 2.0;
    				for (i1 = 0; i1 < n_p; i1++) {
    					kx = (int)Math.round(x1);
    					k1 = fm.stringWidth(x_title[i1]);
    					g.drawString(x_title[i1], kx-k1/2, y_d+6*f_size/5);
    					g.drawLine(kx, y_d, kx, y_d-5);
    					x1 += sp;
    				}
    			}
    			else {
    				x1 = x_scale[0];
    				y1 = x_l;
    				sp = (double)(x_r - x_l) / k_x;
    				for (i1 = 0; i1 < k_x+1; i1++) {
    					kx = (int)Math.round(y1);
    					s1 = df_x.format(x1);
    					k1 = fm.stringWidth(s1);
    					g.drawString(s1, kx-k1/2, y_d+6*f_size/5);
    					g.drawLine(kx, y_d, kx, y_u);
    					x1 += x_scale[2];
    					y1 += sp;
    				}
    			}
    		}
    						// 横表示
    		else {
    							// y軸
    			if (type) {
    				n_p = x_title.length;
    				len = 0;
    				for (i1 = 0; i1 < n_p; i1++) {
    					k1 = fm.stringWidth(x_title[i1]);
    					if (k1 > len)
    						len = k1;
    				}
    				g.drawLine(x_l+len+5, y_u, x_l+len+5, y_d);
    				g.drawLine(x_r, y_u, x_r, y_d);
    				sp = (double)(y_d - y_u) / n_p;
    				x1 = y_d - sp / 2.0;
    				for (i1 = 0; i1 < n_p; i1++) {
    					ky = (int)Math.round(x1);
    					k1 = fm.stringWidth(x_title[n_p-1-i1]);
    					g.drawString(x_title[n_p-1-i1], x_l+len-k1, ky+f_size/2);
    					g.drawLine(x_l+len+5, ky, x_l+len+10, ky);
    					x1 -= sp;
    				}
    				g.drawLine(x_l+len+5, y_u, x_r, y_u);
    				g.drawLine(x_l+len+5, y_d, x_r, y_d);
    				x_l += (len + 5);
    			}
    			else {
    				y1  = x_scale[0];
    				len = 0;
    				for (i1 = 0; i1 < k_x+1; i1++) {
    					s1 = df_x.format(y1);
    					k1 = fm.stringWidth(s1);
    					if (k1 > len)
    						len = k1;
    					y1 += x_scale[2];
    				}
    				g.drawLine(x_l+len+5, y_u, x_l+len+5, y_d);
    				g.drawLine(x_r, y_u, x_r, y_d);
    				y1 = x_scale[0];
    				x1 = y_d;
    				sp = (double)(y_d - y_u) / k_x;
    				for (i1 = 0; i1 < k_x+1; i1++) {
    					ky = (int)Math.round(x1);
    					s1 = df_x.format(y1);
    					k1 = fm.stringWidth(s1);
    					g.drawString(s1, x_l+len-k1, ky+f_size/2);
    					g.drawLine(x_l+len+5, ky, x_r, ky);
    					y1 += x_scale[2];
    					x1 -= sp;
    				}
    				x_l += (len + 5);
    			}
    							// x軸
    			x1 = y_scale[0];
    			y1 = x_l;
    			sp = (double)(x_r - x_l) / k_y;
    			for (i1 = 0; i1 < k_y+1; i1++) {
    				kx = (int)Math.round(y1);
    				s1 = df_y.format(x1);
    				k1 = fm.stringWidth(s1);
    				g.drawString(s1, kx-k1/2, y_d+6*f_size/5);
    				g.drawLine(kx, y_d, kx, y_u);
    				x1 += y_scale[2];
    				y1 += sp;
    			}
    		}
    					//
    					// グラフの表示
    					//
    		g2.setStroke(new BasicStroke(line_w));
    		cr = (int)line_w + 6;
    						// 縦表示
    		if (ver) {
    			if (type) {
    				n_p = x_title.length;
    				sp  = (double)(x_r - x_l) / n_p;
    				k1  = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					x1  = x_l + sp / 2.0;
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						kx = (int)Math.round(x1);
    						ky = y_d - (int)((y_d - y_u) * (data_y[i1][i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1  = kx;
    						ky1  = ky;
    						x1  += sp;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    			else {
    				n_p = data_x[0].length;
    				k1  = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						kx = x_l + (int)((x_r - x_l) * (data_x[i1][i2] - x_scale[0]) / (x_scale[1] - x_scale[0]));
    						ky = y_d - (int)((y_d - y_u) * (data_y[i1][i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1 = kx;
    						ky1 = ky;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    		}
    						// 横表示
    		else {
    			if (type) {
    				n_p = x_title.length;
    				sp = (double)(y_d - y_u) / n_p;
    				k1 = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					y1  = y_d - sp / 2.0;
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						ky = (int)Math.round(y1);
    						kx = x_l + (int)((x_r - x_l) * (data_y[i1][n_p-1-i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1  = kx;
    						ky1  = ky;
    						y1  -= sp;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    			else {
    				n_p = data_x[0].length;
    				k1  = 0;
    				for (i1 = 0; i1 < n_g; i1++) {
    					g.setColor(cl[k1]);
    					kx1 = 0;
    					ky1 = 0;
    					for (i2 = 0; i2 < n_p; i2++) {
    						kx = x_l + (int)((x_r - x_l) * (data_y[i1][i2] - y_scale[0]) / (y_scale[1] - y_scale[0]));
    						ky = y_d - (int)((y_d - y_u) * (data_x[i1][i2] - x_scale[0]) / (x_scale[1] - x_scale[0]));
    						if (line_m)
    							g.fillOval(kx-cr/2, ky-cr/2, cr, cr);
    						if (i2 > 0)
    							g.drawLine(kx1, ky1, kx, ky);
    						kx1 = kx;
    						ky1 = ky;
    					}
    					k1++;
    					if (k1 >= cl.length)
    						k1 = 0;
    				}
    			}
    		}
    		g2.setStroke(new BasicStroke(1.0f));
    	}
    
    	/************************************/
    	/* マウスがクリックされたときの処理 */
    	/************************************/
    	class ClickMouse extends MouseAdapter
    	{
    		Draw_line dr;
    
    		ClickMouse(Draw_line dr1)
    		{
    			dr = dr1;
    		}
    
    		public void mouseClicked(MouseEvent e)
    		{
    			int xp = e.getX();
    			int yp = e.getY();
    					// 縦表示と横表示の変換
    			if (xp > bx1 && xp < bx1+(bx2-bx1)/2 && yp > by1 && yp < by2) {
    				if (ver) {
    					ver    = false;
    					change = "縦 色";
    				}
    				else {
    					ver    = true;
    					change = "横 色";
    				}
    				repaint();
    			}
    					// グラフの色,線の太さ等
    			else if (xp > bx1+(bx2-bx1)/2 && xp < bx2 && yp > by1 && yp < by2) {
    				Modify md = new Modify(dr.line, dr);
    				md.setVisible(true);
    			}
    		}
    	}
    }
    
    /****************************/
    /* 色及び線の太さの変更     */
    /*      coded by Y.Suganuma */
    /****************************/
    class Modify extends JDialog implements ActionListener, TextListener {
    	Draw_line dr;   // 折れ線グラフ
    	JButton bt_dr;
    	TextField rgb[];
    	TextField r[];
    	TextField g[];
    	TextField b[];
    	JTextField tx;
    	JRadioButton r1, r2;
    	float line_w = 1.0f;   // 折れ線グラフ等の線の太さ
    	boolean line_m = true;   // 折れ線グラフ等にマークを付けるか否か
    	Color cl[];   // グラフの色
    	int n_g;   // グラフの数
    	int wd;   // 線の太さを変更するか
    	int mk;   // マークを変更するか
    	int n;
    	JPanel jp[];
    					// 折れ線グラフ
    	Modify(Frame host, Draw_line dr1)
    	{
    		super(host, "色と線の変更", true);
    							// 初期設定
    		dr  = dr1;
    		wd  = 1;
    		mk  = 1;
    		n_g = dr.n_g;
    		if (n_g > 10)
    			n_g = 10;
    		n      = n_g + 3;
    		line_w = dr.line_w;
    		line_m = dr.line_m;
    		cl     = new Color[n_g];
    		for (int i1 = 0; i1 < n_g; i1++)
    			cl[i1] = dr.cl[i1];
    		set();
    							// ボタン
    		Font f = new Font("TimesRoman", Font.BOLD, 20);
    		bt_dr = new JButton("OK");
    		bt_dr.setFont(f);
    		bt_dr.addActionListener(this);
    		jp[n-1].add(bt_dr);
    	}
    					// 設定
    	void set()
    	{
    		setSize(450, 60*(n));
    		Container cp = getContentPane();
    		cp.setBackground(Color.white);
    		cp.setLayout(new GridLayout(n, 1, 5, 5));
    		jp = new JPanel[n];
    		for (int i1 = 0; i1 < n; i1++) {
    			jp[i1] = new JPanel();
    			cp.add(jp[i1]);
    		}
    		Font f = new Font("TimesRoman", Font.BOLD, 20);
    							// 色の変更
    		JLabel lb[][] = new JLabel[n_g][3];
    		rgb = new TextField[n_g];
    		r = new TextField[n_g];
    		g = new TextField[n_g];
    		b = new TextField[n_g];
    		for (int i1 = 0; i1 < n_g; i1++) {
    			rgb[i1] = new TextField(3);
    			rgb[i1].setFont(f);
    			rgb[i1].setBackground(new Color(cl[i1].getRed(), cl[i1].getGreen(), cl[i1].getBlue()));
    			jp[i1].add(rgb[i1]);
    			lb[i1][0] = new JLabel(" 赤");
    			lb[i1][0].setFont(f);
    			jp[i1].add(lb[i1][0]);
    			r[i1] = new TextField(3);
    			r[i1].setFont(f);
    			r[i1].setBackground(Color.white);
    			r[i1].setText(Integer.toString(cl[i1].getRed()));
    			r[i1].addTextListener(this);
    			jp[i1].add(r[i1]);
    			lb[i1][1] = new JLabel("緑");
    			lb[i1][1].setFont(f);
    			jp[i1].add(lb[i1][1]);
    			g[i1] = new TextField(3);
    			g[i1].setFont(f);
    			g[i1].setBackground(Color.white);
    			g[i1].setText(Integer.toString(cl[i1].getGreen()));
    			g[i1].addTextListener(this);
    			jp[i1].add(g[i1]);
    			lb[i1][2] = new JLabel("青");
    			lb[i1][2].setFont(f);
    			jp[i1].add(lb[i1][2]);
    			b[i1] = new TextField(3);
    			b[i1].setFont(f);
    			b[i1].setBackground(Color.white);
    			b[i1].setText(Integer.toString(cl[i1].getBlue()));
    			b[i1].addTextListener(this);
    			jp[i1].add(b[i1]);
    		}
    							// 線の変更
    		if (wd > 0) {
    			JLabel lb1 = new JLabel("線の太さ:");
    			lb1.setFont(f);
    			jp[n_g].add(lb1);
    			tx = new JTextField(2);
    			tx.setFont(f);
    			tx.setBackground(Color.white);
    			tx.setText(Integer.toString((int)line_w));
    			jp[n_g].add(tx);
    		}
    
    		if (mk > 0) {
    			JLabel lb2 = new JLabel("マーク:");
    			lb2.setFont(f);
    			jp[n-2].add(lb2);
    			ButtonGroup gp = new ButtonGroup();
    			r1 = new JRadioButton("付ける");
    			r1.setFont(f);
    			gp.add(r1);
    			jp[n-2].add(r1);
    			r2 = new JRadioButton("付けない");
    			r2.setFont(f);
    			gp.add(r2);
    			jp[n-2].add(r2);
    			if (line_m)
    				r1.doClick();
    			else
    				r2.doClick();
    		}
    	}
    					// TextFieldの内容が変更されたときの処理
    	public void textValueChanged(TextEvent e)
    	{
    		for (int i1 = 0; i1 < n_g; i1++) {
    			if (e.getSource() == r[i1] || e.getSource() == g[i1] || e.getSource() == b[i1]) {
    				String str = r[i1].getText();
    				int rc = str.length()>0 ? Integer.parseInt(str) : 0;
    				str = g[i1].getText();
    				int gc = str.length()>0 ? Integer.parseInt(str) : 0;
    				str = b[i1].getText();
    				int bc = str.length()>0 ? Integer.parseInt(str) : 0;
    				rgb[i1].setBackground(new Color(rc, gc, bc));
    			}
    		}
    	}
    					// 値の設定
    	public void actionPerformed(ActionEvent e)
    	{
    		for (int i1 = 0; i1 < n_g; i1++) {
    			String str = r[i1].getText();
    			int rc = str.length()>0 ? Integer.parseInt(str) : 0;
    			str = g[i1].getText();
    			int gc = str.length()>0 ? Integer.parseInt(str) : 0;
    			str = b[i1].getText();
    			int bc = str.length()>0 ? Integer.parseInt(str) : 0;
    			dr.cl[i1] = new Color(rc, gc, bc);
    		}
    		dr.line_w = Integer.parseInt(tx.getText());
    		if (r1.isSelected())
    			dr.line_m = true;
    		else
    			dr.line_m = false;
    		dr.repaint();
    
    		setVisible(false);
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,任意のデータに対する処理を画面上で実行することが可能であり,結果はテキストエリアに出力されると共に,「確率(複数点)」を選択すればグラフも表示されます.グラフを表示する場合,クリックしたページ( test.html )とは異なるページにグラフを表示するため,test.html に入力されたデータを graph_js.htm に送ってグラフを表示しています.なお,graph_js.htm において利用している JavaScrip のプログラムファイル graph.js に関しては,「グラフの表示」を参照して下さい.
    test.html(データの設定)
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>F分布</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()
    		{
    			dof1 = parseInt(document.getElementById("dof1").value);
    			dof2 = parseInt(document.getElementById("dof2").value);
    			let z = new Array();
    			let sw1 = new Array();
    					// 1点
    			if (method == 0) {
    				let x = parseFloat(document.getElementById("x").value);
    				let y = f(x, z, dof1, dof2);
    				let str = "P(X = " + x + ") = " + z[0] + "\n";
    				str = str + "P(X < " + x + ") = " + y;
    				document.getElementById("xx").value = str;
    			}
    					// 複数点
    			else if (method == 1) {
    							// データ取得
    				let from = parseFloat(document.getElementById("from").value);
    				let to = parseFloat(document.getElementById("to").value);
    				let h  = parseFloat(document.getElementById("h").value);
    				let x1 = new Array();
    				let y1 = new Array();
    				let y2 = new Array();
    				let mx = 0;
    				let n = 0;
    				for (let x = from; x < to+0.5*h; x += h) {
    					let y = f(x, z, dof1, dof2);
    					x1[n] = x;
    					y1[n] = z[0];
    					y2[n] = y;
    					if (z[0] > mx)
    						mx = z[0];
    					n++;
    				}
    				let str1 = "密度関数\n";
    				let str2 = "分布関数\n";
    				for (let i1 = 0; i1 < n; i1++) {
    					str1 = str1 + x1[i1] + " " + y1[i1] + "\n";
    					str2 = str2 + x1[i1] + " " + y2[i1] + "\n";
    				}
    				document.getElementById("xx").value = str1;
    				document.getElementById("yy").value = str2;
    							// グラフの描画
    				let gp1 = "2,F分布(密度関数),x,f(x),1,グラフ1," + from + "," + to;
    				let sp1 = Math.floor((to - from) / 5 * 10);
    				if (sp1 % 5 > 0)
    					sp1 = sp1 + 5 - sp1 % 5;
    				sp1 /= 10;
    				let sp2 = Math.floor(mx / 5 * 100);
    				if (sp2 % 5 > 0)
    					sp2 = sp2 + 5 - sp2 % 5;
    				sp2 /= 100;
    				mx = sp2 * 5;
    				gp1 = gp1 + "," + sp1 + ",1,0.0," + mx + "," + sp2 + ",2," + n;
    				for (let i1 = 0; i1 < n; i1++)
    					gp1 = gp1 + "," + Math.round(100 * x1[i1]) / 100;
    				for (let i1 = 0; i1 < n; i1++)
    					gp1 = gp1 + "," + Math.round(1000 * y1[i1]) / 1000;
    				gp1 = gp1 + ",1,0";
    				let str = "graph_js.htm?gp=" + gp1;
    				open(str, "density", "width=950, height=700");
    
    				let gp2 = "2,F分布(分布関数),x,F(x),1,グラフ1," + from + "," + to;
    				sp2 = 0.2;
    				mx = 1.0;
    				gp2 = gp2 + "," + sp1 + ",1,0.0," + mx + "," + sp2 + ",1," + n;
    				for (let i1 = 0; i1 < n; i1++)
    					gp2 = gp2 + "," + Math.round(100 * x1[i1]) / 100;
    				for (let i1 = 0; i1 < n; i1++)
    					gp2 = gp2 + "," + Math.round(1000 * y2[i1]) / 1000;
    				gp2 = gp2 + ",1,0";
    				str = "graph_js.htm?gp=" + gp2;
    				open(str, "distribution", "width=950, height=700");
    			}
    					// %値
    			else {
    				let x = parseFloat(document.getElementById("p").value);
    				p = x / 100.0;
    				if (p < 1.0e-7) {
    					str = x + "%値 = ∞";
    					document.getElementById("xx").value = str;
    				}
    				else if ((1.0-p) < 1.0e-7) {
    					str = x + "%値 = 0";
    					document.getElementById("xx").value = str;
    				}
    				else {
    					let y = p_f(sw1);
    					if (sw1[0] < 0)
    						document.getElementById("xx").value = "収束しませんでした";
    					else {
    						str = x + "%値 = " + y;
    						document.getElementById("xx").value = str;
    					}
    				}
    			}
    		}
    
    		/****************/
    		/* 関数値の計算 */
    		/****************/
    		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;
    		}
    
    		/***************************************/
    		/*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-10, 1.0e-15, 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 < 3; i1++) {
    				if (sel.options[i1].selected) {
    					method = i1;
    					break;
    				}
    			}
    							// 1点
    			if (method == 0) {
    				document.getElementById("x_t").style.display = "";
    				document.getElementById("p_t").style.display = "none";
    				document.getElementById("from_t").style.display = "none";
    				document.getElementById("to_t").style.display = "none";
    				document.getElementById("h_t").style.display = "none";
    				document.getElementById("yy").style.display = "none";
    			}
    							// 複数点
    			else if (method == 1) {
    				document.getElementById("x_t").style.display = "none";
    				document.getElementById("p_t").style.display = "none";
    				document.getElementById("from_t").style.display = "";
    				document.getElementById("to_t").style.display = "";
    				document.getElementById("h_t").style.display = "";
    				document.getElementById("yy").style.display = "";
    			}
    							// %値
    			else {
    				document.getElementById("x_t").style.display = "none";
    				document.getElementById("p_t").style.display = "";
    				document.getElementById("from_t").style.display = "none";
    				document.getElementById("to_t").style.display = "none";
    				document.getElementById("h_t").style.display = "none";
    				document.getElementById("yy").style.display = "none";
    			}
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; text-align:center; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>F分布</B></H2>
    
    	自由度1:<INPUT ID="dof1" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="6"> 
    	自由度2:<INPUT ID="dof2" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="4"> 
    	<SPAN ID="x_t">x:<INPUT ID="x" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="1.0"> </SPAN>
    	<SPAN ID="p_t" STYLE="display: none">%値:<INPUT ID="p" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="2.5"> </SPAN>
    	<SPAN ID="from_t" STYLE="display: none">下限:<INPUT ID="from" DISABLED STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="0.0"> </SPAN>
    	<SPAN ID="to_t" STYLE="display: none">上限:<INPUT ID="to" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="10.0"> </SPAN>
    	<SPAN ID="h_t" STYLE="display: none">刻み幅:<INPUT ID="h" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="0.1"> </SPAN>
    	目的:<SELECT ID="method" onChange="set_m()" STYLE="font-size:100%">
    		<OPTION SELECTED>確率(1点)
    		<OPTION>確率(複数点)
    		<OPTION>p%値
    	</SELECT> 
    	<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
    	<TEXTAREA ID="xx" COLS="40" ROWS="30" STYLE="font-size: 100%"></TEXTAREA>
    	<TEXTAREA ID="yy" COLS="40" ROWS="30" STYLE="font-size: 100%; display: none"></TEXTAREA>
    </BODY>
    
    </HTML>
    			
    graph_js.htm(グラフの表示)
    <!DOCTYPE HTML>
    <HTML>
    <HEAD>
    	<TITLE>グラフの表示</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<LINK REL="stylesheet" TYPE="text/css" HREF="../master.css">
    	<SCRIPT TYPE="text/javascript" SRC="graph.js"></SCRIPT>
    	<SCRIPT TYPE="text/javascript">
    		function GetParameter()
    		{
    			let result = new Array();
    			if(1 < window.location.search.length) {
    							// 最初の1文字 (?記号) を除いた文字列を取得する
    				let str = window.location.search.substring(1);
    							// 区切り記号 (&) で文字列を配列に分割する
    				let param = str.split('&');
    				for(let i1 = 0; i1 < param.length; i1++ ) {
    							// パラメータ名とパラメータ値に分割する
    					let element = param[i1].split('=');
    					let Name = decodeURIComponent(element[0]);
    					let Value = decodeURIComponent(element[1]);
    							// パラメータ名をキーとして連想配列に追加する
    					result[Name] = Value;
    				}
    			}
    			return result;
    		}
    	</SCRIPT>
    </HEAD>
    <BODY CLASS="white" STYLE="text-align: center">
    	<DIV ID="cl_line" STYLE="text-align: center; display: none">
    		<FORM>
    			<DIV ID="c0">
    				<INPUT ID="rgb0" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r0" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(0)"> 
    				緑<INPUT ID="g0" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(0)"> 
    				青<INPUT ID="b0" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(0)">
    			</DIV>
    			<DIV ID="c1">
    				<INPUT ID="rgb1" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r1" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(1)"> 
    				緑<INPUT ID="g1" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(1)"> 
    				青<INPUT ID="b1" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(1)">
    			</DIV>
    			<DIV ID="c2">
    				<INPUT ID="rgb2" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r2" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(2)"> 
    				緑<INPUT ID="g2" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(2)"> 
    				青<INPUT ID="b2" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(2)">
    			</DIV>
    			<DIV ID="c3">
    				<INPUT ID="rgb3" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r3" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(3)"> 
    				緑<INPUT ID="g3" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(3)"> 
    				青<INPUT ID="b3" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(3)">
    			</DIV>
    			<DIV ID="c4">
    				<INPUT ID="rgb4" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r4" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(4)"> 
    				緑<INPUT ID="g4" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(4)"> 
    				青<INPUT ID="b4" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(4)">
    			</DIV>
    			<DIV ID="c5">
    				<INPUT ID="rgb5" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r5" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(5)"> 
    				緑<INPUT ID="g5" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(5)"> 
    				青<INPUT ID="b5" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(5)">
    			</DIV>
    			<DIV ID="c6">
    				<INPUT ID="rgb6" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r6" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(6)"> 
    				緑<INPUT ID="g6" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(6)"> 
    				青<INPUT ID="b6" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(6)">
    			</DIV>
    			<DIV ID="c7">
    				<INPUT ID="rgb7" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r7" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(7)"> 
    				緑<INPUT ID="g7" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(7)"> 
    				青<INPUT ID="b7" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(7)">
    			</DIV>
    			<DIV ID="c8">
    				<INPUT ID="rgb8" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r8" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(8)"> 
    				緑<INPUT ID="g8" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(8)"> 
    				青<INPUT ID="b8" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(8)">
    			</DIV>
    			<DIV ID="c9">
    				<INPUT ID="rgb9" TYPE="text" SIZE="3" STYLE="font-size: 90%"> 
    				赤<INPUT ID="r9" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(9)"> 
    				緑<INPUT ID="g9" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(9)"> 
    				青<INPUT ID="b9" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(9)">
    			</DIV>
    			<DIV ID="line_m">
    				マーク:<INPUT ID="l_m1" TYPE="radio" NAME="mark" STYLE="font-size: 90%" CHECKED>付ける 
    				<INPUT ID="l_m2" TYPE="radio" NAME="mark" STYLE="font-size: 90%">付けない
    			</DIV>
    			<DIV ID="line_w">
    				線の太さ:<INPUT ID="l_w" TYPE="text" SIZE="3" STYLE="font-size: 90%" VALUE="2">
    			</DIV>
    			<DIV>
    				<SPAN STYLE="background-color: pink; font-size: 100%" onClick="D_Change()">OK</SPAN>
    			</DIV>
    		</FORM>
    	</DIV>
    	<BR>
    	<DIV STYLE="text-align: center">
    		<CANVAS ID="canvas_e" STYLE="background-color: #eeffee;" WIDTH="900" HEIGHT="600" onClick="Click(event)"></CANVAS>
    	</DIV>
    	<SCRIPT TYPE="text/javascript">
    		let result = GetParameter();
    		graph(result['gp']);
    	</SCRIPT>
    </BODY>
    </HTML>
    			

  4. PHP

    <?php
    
    /****************************/
    /* F分布の計算             */
    /*      coded by Y.Suganuma */
    /****************************/
    
    /*********************************************************/
    /* 二分法による非線形方程式(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;
    }
    
    /********/
    /* main */
    /********/
    
    	printf("自由度1は? ");
    	fscanf(STDIN, "%d", $dof1);
    	printf("自由度2は? ");
    	fscanf(STDIN, "%d", $dof2);
    	printf("目的とする結果は? \n");
    	printf("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)\n");
    	printf("     =1 : p%値( P(X > u) = 0.01p となるuの値) ");
    	fscanf(STDIN, "%d", $sw);
    
    	if ($sw == 0) {
    
    		printf("グラフ出力?(=1: yes,  =0: no) ");
    		fscanf(STDIN, "%d", $sw);
    					// 密度関数と分布関数の値
    		if ($sw == 0) {
    			printf("   データは? ");
    			fscanf(STDIN, "%lf", $x);
    			$y = f($x, $z, $dof1, $dof2);
    			printf("P(X = %f) = %f,  P(X < %f) = %f (自由度 = %d, %d)\n", $x, $z, $x, $y, $dof1, $dof2);
    		}
    					// グラフ出力
    		else {
    			printf("   密度関数のファイル名は? ");
    			fscanf(STDIN, "%s", $file1);
    			printf("   分布関数のファイル名は? ");
    			fscanf(STDIN, "%s", $file2);
    			$out1 = fopen($file1, "wb");
    			$out2 = fopen($file2, "wb");
    			printf("   データの下限は? ");
    			fscanf(STDIN, "%lf", $x1);
    			printf("   データの上限は? ");
    			fscanf(STDIN, "%lf", $x2);
    			printf("   刻み幅は? ");
    			fscanf(STDIN, "%lf", $h);
    			for ($x = $x1; $x < $x2+0.5*$h; $x += $h) {
    				$y = f($x, $z, $dof1, $dof2);
    				fwrite($out1, $x." ".$z."\n");
    				fwrite($out2, $x." ".$y."\n");
    			}
    		}
    	}
    					// %値
    	else {
    		printf("%の値は? ");
    		fscanf(STDIN, "%lf", $x);
    		$p = 0.01 * $x;
    		if ($p < 1.0e-7)
    			printf("%f%値 = ∞(自由度 = %d, %d)\n", $x, $dof1, $dof2);
    		else if ((1.0-$p) < 1.0e-7)
    			printf("%f%値 = 0(自由度 = %d, %d)\n", $x, $dof1, $dof2);
    		else {
    			$y = p_f($sw);
    			printf("%f%値 = %f  sw %d (自由度 = %d, %d)\n", $x, $y, $sw, $dof1, $dof2);
    		}
    	}
    
    ?>
    			

  5. Ruby

    ############################
    # f分布の計算
    #      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
    
    			# 密度関数と分布関数の値
    print("自由度1は? ")
    $dof1 = Integer(gets())
    print("自由度2は? ")
    $dof2 = Integer(gets())
    print("目的とする結果は? \n")
    print("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)\n")
    print("     =1 : p%値( P(X > u) = 0.01p となるuの値) ")
    sw = Integer(gets())
    z  = Array.new(1)
    
    if sw == 0
    	print("グラフ出力?(=1: yes,  =0: no) ")
    	sw = Integer(gets())
    	if sw == 0
    			# 密度関数と分布関数の値
    		print("   データは? ")
    		x = Float(gets())
    		y = f(x, z, $dof1, $dof2)
    		print("P(X = " + String(x) + ") = " + String(z[0]) + ",  P( X < " + String(x) + ") = " + String(y) + "(自由度 = " + String($dof1) + ", " + String($dof2) + ")\n")
    			# グラフ出力
    	else
    		print("   密度関数のファイル名は? ")
    		file1 = gets().strip()
    		print("   分布関数のファイル名は? ")
    		file2 = gets().strip()
    		print("   データの下限は? ")
    		x1 = Float(gets())
    		print("   データの上限は? ")
    		x2 = Float(gets())
    		print("   刻み幅は? ")
    		h    = Float(gets())
    		out1 = open(file1, "w")
    		out2 = open(file2, "w")
    		x    = x1
    		while x < x2+0.5*h
    			y = f(x, z, $dof1, $dof2)
    			out1.print(String(x) + " " + String(z[0]) + "\n")
    			out2.print(String(x) + " " + String(y) + "\n")
    			x += h
    		end
    		out1.close()
    		out2.close()
    	end
    			# %値
    else
    	print("%の値は? ")
    	x  = Float(gets())
    	$p = 0.01 * x
    	if $p < 1.0e-7
    		print(String(x) + "%値 = ∞ (自由度 = " + String($dof1) + ", " + String($dof2) + ")\n")
    	elsif (1.0-$p)< 1.0e-7
    		print(String(x) + "%値 = 0 (自由度 = " + String($dof1) + ", " + String($dof2) + ")\n")
    	else
    		ind = Array.new(1)
    		y   = p_f(ind)
    		print(String(x) + "%値 = " + String(y) + "  収束 " + String(ind[0]) + " (自由度 = " + String($dof1) + ", " + String($dof2) + ")\n")
    	end
    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
    
    ############################
    # f分布の計算
    #      coded by Y.Suganuma
    ############################
    
    ##########################################
    # 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
    
    			# 密度関数と分布関数の値
    s    = input("自由度1は? ")
    dof1 = int(s)
    s    = input("自由度2は? ")
    dof2 = int(s)
    print("目的とする結果は? ")
    print("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)")
    s  = input("     =1 : p%値( P(X > u) = 0.01p となるuの値) ")
    sw = int(s)
    z  = np.empty(1, np.float)
    
    if sw == 0 :
    	s  = input("グラフ出力?(=1: yes,  =0: no) ")
    	sw = int(s)
    	if sw == 0 :
    			# 密度関数と分布関数の値
    		s  = input("   データは? ")
    		x  = float(s)
    		y  = f(x, z, dof1, dof2)
    		print("P(X = " + str(x) + ") = " + str(z[0]) + ",  P( X < " + str(x) + ") = " + str(y) + "(自由度 = " + str(dof1) + ", " + str(dof2) + ")")
    			# グラフ出力
    	else :
    		file1 = input("   密度関数のファイル名は? ")
    		file2 = input("   分布関数のファイル名は? ")
    		s     = input("   データの下限は? ")
    		x1    = float(s)
    		s     = input("   データの上限は? ")
    		x2    = float(s)
    		s     = input("   刻み幅は? ")
    		h     = float(s)
    		out1  = open(file1, "w")
    		out2  = open(file2, "w")
    		x     = x1
    		while x < x2+0.5*h :
    			y = f(x, z, dof1, dof2)
    			out1.write(str(x) + " " + str(z[0]) + "\n")
    			out2.write(str(x) + " " + str(y) + "\n")
    			x += h
    		out1.close()
    		out2.close()
    			# %値
    else :
    	s = input("%の値は? ")
    	x = float(s)
    	p = 0.01 * x
    	if p < 1.0e-7 :
    		print(str(x) + "%値 = ∞ (自由度 = " + str(dof1) + ", " + str(dof2) + ")")
    	elif (1.0-p)< 1.0e-7 :
    		print(str(x) + "%値 = 0 (自由度 = " + str(dof1) + ", " + str(dof2) + ")")
    	else :
    		ind = np.empty(1, np.int)
    		y   = p_f(ind)
    		print(str(x) + "%値 = " + str(y) + "  収束 " + str(ind[0]) + " (自由度 = " + str(dof1) + ", " + str(dof2) + ")")
    			

  7. C#

    /****************************/
    /* F分布の計算             */
    /*      coded by Y.Suganuma */
    /****************************/
    using System;
    using System.IO;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	double p;     // α%値を計算するとき時α/100を設定
    	int dof1, dof2;      // 自由度
    
    	public Test1()
    	{
    		Console.Write("自由度1は? ");
    		dof1 = int.Parse(Console.ReadLine());
    		Console.Write("自由度2は? ");
    		dof2 = int.Parse(Console.ReadLine());
    		Console.WriteLine("目的とする結果は? ");
    		Console.WriteLine("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)");
    		Console.Write("     =1 : p%値( P(X > u) = 0.01p となるuの値) ");
    		int sw = int.Parse(Console.ReadLine());
    
    		if (sw == 0) {
    
    			Console.Write("グラフ出力?(=1: yes,  =0: no) ");
    			sw = int.Parse(Console.ReadLine());
    
    			double z = 0.0;
    					// 密度関数と分布関数の値
    			if (sw == 0) {
    				Console.Write("   データは? ");
    				double x = double.Parse(Console.ReadLine());
    				double y = f(x, ref z, dof1, dof2);
    				Console.WriteLine("P(X = " + x + ") = " + z + ",  P( X < " + x + ") = " + y +
                                      " (自由度 = " + dof1 + " " + dof2 + ")");
    			}
    					// グラフ出力
    			else {
    				Console.Write("   密度関数のファイル名は? ");
    				String file1 = Console.ReadLine();
    				Console.Write("   分布関数のファイル名は? ");
    				String file2 = Console.ReadLine();
    				Console.Write("   データの上限は? ");
    				double to = double.Parse(Console.ReadLine());
    				Console.Write("   刻み幅は? ");
    				double h = double.Parse(Console.ReadLine());
    				StreamWriter out1 = new StreamWriter(file1);
    				StreamWriter out2 = new StreamWriter(file2);
    							// データ取得
    				for (double x = 0.0; x < to+0.5*h; x += h) {
    					double y = f(x, ref z, dof1, dof2);
    					out1.WriteLine(x + " " + z);
    					out2.WriteLine(x + " " + y);
    				}
    				out1.Close();
    				out2.Close();
    			}
    		}
    					// %値
    		else {
    			Console.Write("%の値は? ");
    			double x = double.Parse(Console.ReadLine());
    			p        = 0.01 * x;
    			if (p < 1.0e-7)
    				Console.WriteLine(x + "%値 = ∞ (自由度 = " + dof1 + " " + dof2 + ")");
    			else if ((1.0-p) < 1.0e-7)
    				Console.WriteLine(x + "%値 = 0 (自由度 = " + dof1 + " " + dof2 + ")");
    			else {
    				int sw1  = 0;
    				double y = p_f(ref sw1);
    				Console.WriteLine(x + "%値 = " + y + "  sw " + sw1 +
                                      " (自由度 = " + dof1 + " " + dof2 + ")");
    			}
    		}
    	}
    
    	/****************************************/
    	/* 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;
    	}
    }
    			

  8. VB

    '**************************'
    ' F分布の計算             '
    '      coded by Y.Suganuma '
    '**************************'
    Imports System.IO
    
    Module Test
    
    	Dim p As Double     ' α%値を計算するとき時α/100を設定
    	Dim dof1 As Integer     ' 自由度
    	Dim dof2 As Integer     ' 自由度
    
    	Sub Main()
    
    		Console.Write("自由度1は? ")
    		dof1 = Integer.Parse(Console.ReadLine())
    		Console.Write("自由度2は? ")
    		dof2 = Integer.Parse(Console.ReadLine())
    		Console.WriteLine("目的とする結果は? ")
    		Console.WriteLine("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)")
    		Console.Write("     =1 : p%値( P(X > u) = 0.01p となるuの値) ")
    		Dim sw As Integer = Integer.Parse(Console.ReadLine())
    
    		If sw = 0
    
    			Console.Write("グラフ出力?(=1: yes,  =0: no) ")
    			sw = Integer.Parse(Console.ReadLine())
    
    			Dim z As Double = 0.0
    					' 密度関数と分布関数の値
    			If sw = 0
    				Console.Write("   データは? ")
    				Dim x As Double = double.Parse(Console.ReadLine())
    				Dim y As Double = f(x, z, dof1, dof2)
    				Console.WriteLine("P(X = " & x & ") = " & z & ",  P( X < " & x & ") = " & y &
                                      " (自由度 = " & dof1 & " " & dof2 & ")")
    					' グラフ出力
    			Else
    				Console.Write("   密度関数のファイル名は? ")
    				Dim file1 As String = Console.ReadLine().Trim()
    				Console.Write("   分布関数のファイル名は? ")
    				Dim file2 As String = Console.ReadLine().Trim()
    				Console.Write("   データの上限は? ")
    				Dim to1 As Double = Double.Parse(Console.ReadLine())
    				Console.Write("   刻み幅は? ")
    				Dim h As Double = Double.Parse(Console.ReadLine())
    				Dim out1 As StreamWriter = new StreamWriter(file1)
    				Dim out2 As StreamWriter = new StreamWriter(file2)
    							' データ取得
    				Dim x As Double = 0.0
    				Do While x < to1+0.5*h
    					Dim y As Double = f(x, z, dof1, dof2)
    					out1.WriteLine(x & " " & z)
    					out2.WriteLine(x & " " & y)
    					x += h
    				Loop
    				out1.Close()
    				out2.Close()
    			End If
    					' %値
    		Else
    			Console.Write("%の値は? ")
    			Dim x As Double = Double.Parse(Console.ReadLine())
    			p               = 0.01 * x
    			If (p < 1.0e-7)
    				Console.WriteLine(x & "%値 = ∞ (自由度 = " & dof1 & " " & dof2 & ")")
    			Else If ((1.0-p) < 1.0e-7)
    				Console.WriteLine(x & "%値 = 0 (自由度 = " & dof1 & " " & dof2 & ")")
    			Else
    						' 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 sw1 As Integer = 0
    				Dim y As Double    = p_f(sw1, f_f, f_df, normal_f)
    				Console.WriteLine(x & "%値 = " & y & "  sw " & sw1 &
    				                  " (自由度 = " & dof1 & " " & dof2 & ")")
    			End If
    		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
    			

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