情報学部 菅沼ホーム 目次 索引

t 分布

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

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

  1. C++ 

    /****************************/
/* ｔ分布の計算             */
/*      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 t(double, double *, int);
double p_t(int *);
double gamma(double, int *);
double newton(double(*)(double), double(*)(double), double, double,
              double, int, int *);
double t_f(double);
double t_df(double);

double p;     // α％値を計算するとき時α/100を設定
int dof;      // 自由度
              // 上の２つの変数は，必ず，この位置で定義しておく必要がある

int main()
{
	double h, x, x1, x2, y, z;
	int sw;
	char file1[100], file2[100];
	FILE *out1, *out2;

	printf("自由度は？ ");
	scanf("%d", &dof);
	printf("目的とする結果は？ \n");
	printf("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n");
	printf("     =1 : ｐ％値（ 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 = t(x, &z, dof);
			printf("P(X = %f) = %f,  P( X < %f) = %f (自由度 = %d)\n", x, z, x, y, dof);
		}
					// グラフ出力
		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 = t(x, &z, dof);
				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)\n", x, dof);
		else if ((1.0-p) < 1.0e-7)
			printf("%f％値 = -∞ (自由度 = %d)\n", x, dof);
		else {
			y = p_t(&sw);
			printf("%f％値 = %f  sw %d (自由度 = %d)\n", x, y, sw, dof);
		}
	}

	return 0;
}

/**************************************************/
/* ｔ分布の計算（P(X = tt), P(X < tt)）           */
/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
/*      dd : P(X = tt)                            */
/*      df : 自由度                               */
/*      return : P(X < tt)                        */
/**************************************************/
#include <math.h>

double t(double tt, double *dd, int df)
{
	double pi = 4.0 * atan(1.0);
	double p, pp, sign, t2, u, x;
	int ia, i1;

	sign = (tt < 0.0) ? -1.0 : 1.0;
	if (fabs(tt) < 1.0e-10)
		tt = sign * 1.0e-10;
	t2 = tt * tt;
	x  = t2 / (t2 + df);

	if(df%2 != 0) {
		u  = sqrt(x*(1.0-x)) / pi;
		p  = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi;
		ia = 1;
	}

	else {
		u  = sqrt(x) * (1.0 - x) / 2.0;
		p  = sqrt(x);
		ia = 2;
	}

	if (ia != df) {
		for (i1 = ia; i1 <= df-2; i1 += 2) {
			p += 2.0 * u / i1;
			u *= (1.0 + i1) / i1 * (1.0 - x);
		}
	}

	*dd = u / fabs(tt);
	pp  = 0.5 + 0.5 * sign * p;

	return pp;
}

/**************************************************/
/* ｔ分布のｐ％値（P(X > u) = 0.01p）             */
/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
/*      ind : >= 0 : normal（収束回数）           */
/*            = -1 : 収束しなかった               */
/**************************************************/
#include <math.h>

double p_t(int *ind)
{
	double pi = 4.0 * atan(1.0);
	double c, df, df2, e, pis, p2, tt = 0.0, t0, x, yq;

	pis = sqrt(pi);
	df  = (double)dof;
	df2 = 0.5 * df;
					// 自由度＝１
	if (dof == 1)
		tt = tan(pi*(0.5-p));

	else {
					// 自由度＝２
		if (dof == 2) {
			c   = (p > 0.5) ? -1.0 : 1.0;
			p2  = (1.0 - 2.0 * p);
			p2 *= p2;
			tt  = c * sqrt(2.0 * p2 / (1.0 - p2));
		}
					// 自由度＞２
		else {

			yq = p_normal(ind);   // 初期値計算のため

			if (*ind >= 0) {

				x  = 1.0 - 1.0 / (4.0 * df);
				e  = x * x - yq * yq / (2.0 * df);

				if (e > 0.5)
					t0 = yq / sqrt(e);
				else {
					x  = sqrt(df) / (pis * p * df * gamma(df2, ind) / gamma(df2+0.5, ind));
					t0 = pow(x, 1.0/df);
				}
						// ニュートン法
				tt = newton(t_f, t_df, t0, 1.0e-6, 1.0e-10, 100, ind);
			}
		}
	}

	return tt;
}

/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/*      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 t_f(double x)
{
	double y;

	return t(x, &y, dof) - 1.0 + p;
}

/**************************/
/* P(X = x)（関数の微分） */
/**************************/
double t_df(double x)
{
	double y, z;

	z = t(x, &y, dof);

	return y;
}

/*****************************************************/
/* Newton法による非線形方程式(f(x)=0)の解            */
/*      f : f(x)を計算する関数名                     */
/*      df : f(x)の微分を計算する関数名              */
/*      x0 : 初期値                                  */
/*      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   */
/*      eps2 : 終了条件２（｜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(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
/*      eps2 : 終了条件２（｜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

    　　グラフ出力を指定すると，指定された範囲の密度関数および分布関数の値をファイルに出力すると共にグラフも表示しますが，グラフ出力を指定しないと，指定された値における密度関数および分布関数の値，または，％値を出力します．なお，グラフに関しては，「グラフの表示」の表示を参照して下さい． 
    /****************************/
/* ｔ分布の計算             */
/*      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 dof;    // 自由度

	public static void main(String args[]) throws IOException
	{
		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
		double h, x, from, to, y;
		double z[] = new double [1];
		int sw, sw1[] = new int [1];

		System.out.print("自由度は？ ");
		dof = Integer.parseInt(in.readLine());
		System.out.println("目的とする結果は？ ");
		System.out.print("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n");
		System.out.print("     =1 : ｐ％値（ 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.t(x, z, dof);
				System.out.println("P(X = " + x + ") = " + z[0] + ",  P( X < " + x + ") = " + y +
                                   " (自由度 = " + dof + ")");
			}
					// グラフ出力
			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("   データの下限は？ ");
				from = Double.parseDouble(in.readLine());
				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 = from; x < to+0.5*h; x += h) {
					y = App.t(x, z, dof);
					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 + "％値 = ∞ (自由度 = " + dof + ")");
			else if ((1.0-p) < 1.0e-7)
				System.out.println(x + "％値 = -∞ (自由度 = " + dof + ")");
			else {
				y = App.p_t(sw1);
				System.out.println(x + "％値 = " + y + "  sw " + sw1[0] + " (自由度 = " + dof + ")");
			}
		}
	}

	/*************************/
	/* グラフの描画          */
	/*      x1 : x座標データ */
	/*      y1 : 密度関数    */
	/*      y2 : 分布関数    */
	/*************************/
	static void graph(ArrayList <Double> x1, ArrayList <Double> y1, ArrayList <Double> y2)
	{
					// 密度関数
		String title1[];   // グラフ，ｘ軸，及び，ｙ軸のタイトル
		String g_title1[];   // 凡例（グラフの内容）
		double x_scale[];   // ｘ軸目盛り
		double y_scale1[];   // ｙ軸目盛り
		double data_x[][], data1_y[][];   // データ
		int n = x1.size();
							// グラフ，ｘ軸，及び，ｙ軸のタイトル
		title1 = new String [3];
		title1[0] = "密度関数（ｔ分布　自由度 = " + dof + "）";
		title1[1] = "x";
		title1[2] = "f(x)";
							// 凡例
		g_title1 = new String [1];
		g_title1[0] = "密度関数";
							// ｘ軸目盛り
		x_scale = new double[3];
		double x_max = (x1.get(n-1)).doubleValue();
		double x_min = (x1.get(0)).doubleValue();
		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_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[];   // グラフ，ｘ軸，及び，ｙ軸のタイトル
		String g_title2[];   // 凡例（グラフの内容）
		double y_scale2[];   // ｙ軸目盛り
		double data2_y[][];   // データ
							// グラフ，ｘ軸，及び，ｙ軸のタイトル
		title2 = new String [3];
		title2[0] = "分布関数（ｔ分布　自由度 = " + dof + "）";
		title2[1] = "x";
		title2[2] = "F(x)";
							// 凡例
		g_title2 = new String [1];
		g_title2[0] = "分布関数";
							// ｙ軸目盛り
		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)）の計算（ｔ分布）
			case 1:
				y = App.t(x, w, Test.dof) - 1.0 + Test.p;
				break;
						// 関数の微分の計算（ｔ分布）
			case 2:
				y = App.t(x, w, Test.dof);
				y = w[0];
				break;
		}

		return y;
	}
}

/************************/
/* 科学技術系算用の手法 */
/************************/
class App {

	/**************************************************/
	/* ｔ分布の計算（P(X = tt), P(X < tt)）           */
	/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
	/*      dd : P(X = tt)                            */
	/*      df : 自由度                               */
	/*      return : P(X < tt)                        */
	/**************************************************/
	static double t(double tt, double dd[], int df)
	{
		double pi = Math.PI;
		double p, pp, sign, t2, u, x;
		int ia, i1;

		sign = (tt < 0.0) ? -1.0 : 1.0;
		if (Math.abs(tt) < 1.0e-10)
			tt = sign * 1.0e-10;
		t2 = tt * tt;
		x  = t2 / (t2 + df);

		if(df%2 != 0) {
			u  = Math.sqrt(x*(1.0-x)) / pi;
			p  = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
			ia = 1;
		}

		else {
			u  = Math.sqrt(x) * (1.0 - x) / 2.0;
			p  = Math.sqrt(x);
			ia = 2;
		}

		if (ia != df) {
			for (i1 = ia; i1 <= df-2; i1 += 2) {
				p += 2.0 * u / i1;
				u *= (1.0 + i1) / i1 * (1.0 - x);
			}
		}

		dd[0] = u / Math.abs(tt);
		pp  = 0.5 + 0.5 * sign * p;

		return pp;
	}

	/**************************************************/
	/* ｔ分布のｐ％値（P(X > u) = 0.01p）             */
	/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
	/*      ind : >= 0 : normal（収束回数）           */
	/*            = -1 : 収束しなかった               */
	/**************************************************/
	static double p_t(int ind[])
	{
		double pi = Math.PI;
		double c, df, df2, e, pis, p2, tt = 0.0, t0, x, yq;

		pis = Math.sqrt(pi);
		df  = (double)Test.dof;
		df2 = 0.5 * df;
					// 自由度＝１
		if (Test.dof == 1)
			tt = Math.tan(pi*(0.5-Test.p));

		else {
					// 自由度＝２
			if (Test.dof == 2) {
				c   = (Test.p > 0.5) ? -1.0 : 1.0;
				p2  = (1.0 - 2.0 * Test.p);
				p2 *= p2;
				tt  = c * Math.sqrt(2.0 * p2 / (1.0 - p2));
			}
					// 自由度＞２
			else {

				yq = p_normal(ind);   // 初期値計算のため

				if (ind[0] >= 0) {

					x  = 1.0 - 1.0 / (4.0 * df);
					e  = x * x - yq * yq / (2.0 * df);

					if (e > 0.5)
						t0 = yq / Math.sqrt(e);
					else {
						x  = Math.sqrt(df) / (pis * Test.p * df * gamma(df2, ind) / gamma(df2+0.5, ind));
						t0 = Math.pow(x, 1.0/df);
					}
						// ニュートン法
					Kansu kn1 = new Kansu(1);
					Kansu kn2 = new Kansu(2);

					tt = newton(t0, 1.0e-6, 1.0e-10, 100, ind, kn1, kn2);
				}
			}
		}

		return tt;
	}

	/****************************************/
	/* Γ（ｘ）の計算（ガンマ関数，近似式） */
	/*      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(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   */
	/*      eps2 : 終了条件２（｜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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
	/*      eps2 : 終了条件２（｜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;

	/*********************************************************/
	/* コンストラクタ（折れ線グラフ２）                      */
	/*      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 : 凡例表示の有無                           */
	/*********************************************************/
	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("折れ線グラフ（２）");
					// 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[];   // ｘ軸への表示
	double x_scale[];   // ｙ軸目盛り
	double y_scale[];   // ｙ軸目盛り
	double data_x[][], data_y[][];   // データ
	boolean d_t;   // タイトル表示の有無
	boolean d_g;   // 凡例表示の有無
	boolean type = true;   // 横軸が項目かデータか
	boolean ver = true;   // 縦か横か
	int place_x;   // 小数点以下の桁数（ｘ軸）
	int place_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;   // グラフの数

	/*********************************************************/
	/* コンストラクタ（折れ線グラフ２）                      */
	/*      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 : 凡例表示の有無                           */
	/*********************************************************/
	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));
					//
					// ｘ軸及びｙ軸のタイトルの表示
					//
		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;
			}
		}
					//
					// ｘ軸，ｙ軸，及び，各軸の目盛り
					//
		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) {
							// ｙ軸
			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);
							// ｘ軸
			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 {
							// ｙ軸
			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);
			}
							// ｘ軸
			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>ｔ分布</TITLE>
	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
	<SCRIPT TYPE="text/javascript">
		method = 0;
		dof = 0;
		p = 0.0;
		function main()
		{
			dof = parseInt(document.getElementById("dof").value);
			let z = new Array();
			let sw1 = new Array();
					// 1点
			if (method == 0) {
				let x = parseFloat(document.getElementById("x").value);
				let y = t(x, z, dof);
				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 = t(x, z, dof);
					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,ｔ分布（密度関数）,x,f(x),1,グラフ１," + from + "," + to;
				let sp1 = (to - from) / 6;
				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,ｔ分布（分布関数）,x,F(x),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 + "％値 = -∞";
					document.getElementById("xx").value = str;
				}
				else {
					let y = p_t(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)）の計算（ｔ分布）
				case 1:
					y = t(x, w, dof) - 1.0 + p;
					break;
							// 関数の微分の計算（ｔ分布）
				case 2:
					y = t(x, w, dof);
					y = w[0];
					break;
			}

			return y;
		}

		/**************************************************/
		/* ｔ分布の計算（P(X = tt), P(X < tt)）           */
		/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
		/*      dd : P(X = tt)                            */
		/*      df : 自由度                               */
		/*      return : P(X < tt)                        */
		/**************************************************/
		function t(tt, dd, df)
		{
			let pi = Math.PI;
			let p;
			let pp;
			let sign;
			let t2;
			let u;
			let x;
			let ia;
			let i1;

			sign = (tt < 0.0) ? -1.0 : 1.0;
			if (Math.abs(tt) < 1.0e-10)
				tt = sign * 1.0e-10;
			t2 = tt * tt;
			x  = t2 / (t2 + df);

			if(df%2 != 0) {
				u  = Math.sqrt(x*(1.0-x)) / pi;
				p  = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
				ia = 1;
			}

			else {
				u  = Math.sqrt(x) * (1.0 - x) / 2.0;
				p  = Math.sqrt(x);
				ia = 2;
			}

			if (ia != df) {
				for (i1 = ia; i1 <= df-2; i1 += 2) {
					p += 2.0 * u / i1;
					u *= (1.0 + i1) / i1 * (1.0 - x);
				}
			}

			dd[0] = u / Math.abs(tt);
			pp  = 0.5 + 0.5 * sign * p;

			return pp;
		}

		/**************************************************/
		/* ｔ分布のｐ％値（P(X > u) = 0.01p）             */
		/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
		/*      ind : >= 0 : normal（収束回数）           */
		/*            = -1 : 収束しなかった               */
		/**************************************************/
		function p_t(ind)
		{
			let pi = Math.PI;
			let c;
			let df;
			let df2;
			let e;
			let pis;
			let p2;
			let tt = 0.0;
			let t0;
			let x;
			let yq;

			pis = Math.sqrt(pi);
			df  = dof;
			df2 = 0.5 * df;
						// 自由度＝１
			if (dof == 1)
				tt = Math.tan(pi*(0.5-p));

			else {
						// 自由度＝２
				if (dof == 2) {
					c   = (p > 0.5) ? -1.0 : 1.0;
					p2  = (1.0 - 2.0 * p);
					p2 *= p2;
					tt  = c * Math.sqrt(2.0 * p2 / (1.0 - p2));
				}
						// 自由度＞２
				else {

					yq = p_normal(ind);   // 初期値計算のため

					if (ind[0] >= 0) {

						x  = 1.0 - 1.0 / (4.0 * df);
						e  = x * x - yq * yq / (2.0 * df);

						if (e > 0.5)
							t0 = yq / Math.sqrt(e);
						else {
							x  = Math.sqrt(df) / (pis * p * df * gamma(df2, ind) / gamma(df2+0.5, ind));
							t0 = Math.pow(x, 1.0/df);
						}
							// ニュートン法
						tt = newton(t0, 1.0e-10, 1.0e-15, 100, ind);
					}
				}
			}

			return tt;
		}

		/****************************************/
		/* Γ（ｘ）の計算（ガンマ関数，近似式） */
		/*      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(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   */
		/*      eps2 : 終了条件２（｜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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   */
		/*      eps2 : 終了条件２（｜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>ｔ分布</B></H2>

	自由度：<INPUT ID="dof" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="10"> 
	<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="5.0"> </SPAN>
	<SPAN ID="from_t" STYLE="display: none">下限：<INPUT ID="from" STYLE="font-size: 100%;" TYPE="text" SIZE="5" VALUE="-3.0"> </SPAN>
	<SPAN ID="to_t" STYLE="display: none">上限：<INPUT ID="to" STYLE="font-size: 100%;" TYPE="text" SIZE="5.0" VALUE="3.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>ｐ％値
	</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

/****************************/
/* ｔ分布の計算             */
/*      coded by Y.Suganuma */
/****************************/


/*********************************************************/
/* 二分法による非線形方程式(f(x)=0)の解                  */
/*      f : f(x)を計算する関数名                         */
/*      x1,x2 : 初期値                                   */
/*      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
/*      eps2 : 終了条件２（｜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(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   */
/*      eps2 : 終了条件２（｜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;
}

/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/*      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;
}

/**************************************************/
/* ｔ分布の計算（P(X = tt), P(X < tt)）           */
/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
/*      dd : P(X = tt)                            */
/*      df : 自由度                               */
/*      return : P(X < tt)                        */
/**************************************************/
function t($tt, &$dd, $df)
{
	$pi   = 4.0 * atan(1.0);
	$sign = ($tt < 0.0) ? -1.0 : 1.0;

	if (abs($tt) < 1.0e-10)
		$tt = $sign * 1.0e-10;
	$t2 = $tt * $tt;
	$x  = $t2 / ($t2 + $df);

	if($df%2 != 0) {
		$u  = sqrt($x*(1.0-$x)) / $pi;
		$p  = 1.0 - 2.0 * atan2(sqrt(1.0-$x), sqrt($x)) / $pi;
		$ia = 1;
	}

	else {
		$u  = sqrt($x) * (1.0 - $x) / 2.0;
		$p  = sqrt($x);
		$ia = 2;
	}

	if ($ia != $df) {
		for ($i1 = $ia; $i1 <= $df-2; $i1 += 2) {
			$p += 2.0 * $u / $i1;
			$u *= (1.0 + $i1) / $i1 * (1.0 - $x);
		}
	}

	$dd = $u / abs($tt);
	$pp = 0.5 + 0.5 * $sign * $p;

	return $pp;
}

/**************************************************/
/* ｔ分布のｐ％値（P(X > u) = 0.01p）             */
/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
/*      ind : >= 0 : normal（収束回数）           */
/*            = -1 : 収束しなかった               */
/**************************************************/
function p_t(&$ind)
{
	global $p, $dof;

	$pi  = 4.0 * atan(1.0);
	$tt  = 0.0;
	$pis = sqrt($pi);
	$df  = floatval($dof);
	$df2 = 0.5 * $df;
					// 自由度＝１
	if ($dof == 1)
		$tt = tan($pi*(0.5-$p));

	else {
					// 自由度＝２
		if ($dof == 2) {
			$c   = ($p > 0.5) ? -1.0 : 1.0;
			$p2  = (1.0 - 2.0 * $p);
			$p2 *= $p2;
			$tt  = $c * sqrt(2.0 * $p2 / (1.0 - $p2));
		}
					// 自由度＞２
		else {

			$yq = p_normal($ind);   // 初期値計算のため

			if ($ind >= 0) {

				$x = 1.0 - 1.0 / (4.0 * $df);
				$e = $x * $x - $yq * $yq / (2.0 * $df);

				if ($e > 0.5)
					$t0 = $yq / sqrt($e);
				else {
					$x  = sqrt($df) / ($pis * $p * $df * gamma($df2, $ind) / gamma($df2+0.5, $ind));
					$t0 = pow($x, 1.0/$df);
				}
						// ニュートン法
				$tt = newton("t_f", "t_df", $t0, 1.0e-6, 1.0e-10, 100, $ind);
			}
		}
	}

	return $tt;
}

/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
function t_f($x)
{
	global $p, $dof;
	return t($x, $y, $dof) - 1.0 + $p;
}

/**************************/
/* P(X = x)（関数の微分） */
/**************************/
function t_df($x)
{
	global $dof;

	$z = t($x, $y, $dof);

	return $y;
}

/********/
/* main */
/********/

	printf("自由度は？ ");
	fscanf(STDIN, "%d", $dof);
	printf("目的とする結果は？ \n");
	printf("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n");
	printf("     =1 : ｐ％値（ 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 = t($x, $z, $dof);
			printf("P(X = %f) = %f,  P( X < %f) = %f (自由度 = %d)\n", $x, $z, $x, $y, $dof);
		}
					// グラフ出力
		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 = t($x, $z, $dof);
				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)\n", $x, $dof);
		else if ((1.0-$p) < 1.0e-7)
			printf("%f％値 = -∞ (自由度 = %d)\n", $x, $dof);
		else {
			$y = p_t($sw);
			printf("%f％値 = %f  sw %d (自由度 = %d)\n", $x, $y, $sw, $dof);
		}
	}

?>

  5. Ruby 

    ############################
# ｔ分布の計算
#      coded by Y.Suganuma
############################

############################################
# Γ（ｘ）の計算（ガンマ関数，近似式）
#      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

################################################
# ｔ分布の計算（P(X = tt), P(X < tt)）
# （自由度が∞の時の値はN(0,1)を利用して下さい）
#      dd : P(X = tt)
#      df : 自由度
#      return : P(X < tt)
################################################

def t(tt, dd, df)

	if tt < 0.0
		sign = -1.0
	else
		sign = 1.0
	end
	if tt.abs() < 1.0e-10
		tt = sign * 1.0e-10
	end
	t2 = tt * tt
	x  = t2 / (t2 + df)

	if df%2 != 0
		u  = Math.sqrt(x*(1.0-x)) / Math::PI
		p  = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / Math::PI
		ia = 1
	else
		u  = Math.sqrt(x) * (1.0 - x) / 2.0
		p  = Math.sqrt(x)
		ia = 2
	end

	if ia != df
		i1 = ia
		while i1 < df-1
			p  += 2.0 * u / i1
			u  *= ((1.0 + i1) / i1 * (1.0 - x))
			i1 += 2
		end
	end

	dd[0] = u / tt.abs()
	pp    = 0.5 + 0.5 * sign * p

	return pp
end

############################################
# 二分法による非線形方程式(f(x)=0)の解
#      x1,x2 : 初期値
#      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
#      eps2 : 終了条件２（｜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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
#      eps2 : 終了条件２（｜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(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

#################################################
# ｔ分布のｐ％値（P(X > u) = 0.01p）
# （自由度が∞の時の値はN(0,1)を利用して下さい）
#      ind : >= 0 : normal（収束回数）
#            = -1 : 収束しなかった
#################################################
def p_t(ind)

	t_snx = Proc.new { |sw, x|
		y = Array.new(1)
		z = t(x, y, $dof)
		if sw == 0
			z - 1.0 + $p
		else
			z = y[0]
		end
	}

	tt  = 0.0
	pis = Math.sqrt(Math::PI)
	df  = Float($dof)
	df2 = 0.5 * df
			# 自由度＝１
	if $dof == 1
		tt = Math.tan(Math::PI*(0.5-$p))

	else
			# 自由度＝２
		if $dof == 2
			if $p > 0.5
				c = -1.0
			else
				c = 1.0
			end
			p2  = (1.0 - 2.0 * $p)
			p2 *= p2
			tt  = c * Math.sqrt(2.0 * p2 / (1.0 - p2))
			# 自由度＞２
		else

			yq = p_normal(ind)   # 初期値計算のため

			if ind[0] >= 0

				x  = 1.0 - 1.0 / (4.0 * df)
				e  = x * x - yq * yq / (2.0 * df)

				if e > 0.5
					t0 = yq / Math.sqrt(e)
				else
					x  = Math.sqrt(df) / (pis * $p * df * gamma(df2, ind) / gamma(df2+0.5, ind))
					t0 = x ** (1.0/df)
				end
					# ニュートン法
				tt = newton(t0, 1.0e-6, 1.0e-10, 100, ind, &t_snx)
			end
		end
	end

	return tt
end

			# 密度関数と分布関数の値
print("自由度は？ ")
$dof = Integer(gets())
print("目的とする結果は？ \n")
print("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n")
print("     =1 : ｐ％値（ 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 = t(x, z, $dof)
		print("P(X = " + String(x) + ") = " + String(z[0]) + ",  P( X < " + String(x) + ") = " + String(y) + "(自由度 = " + String($dof) + ")\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 = t(x, z, $dof)
			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($dof) + ")\n")
	elsif (1.0-$p)< 1.0e-7
		print(String(x) + "％値 = -∞ (自由度 = " + String($dof) + ")\n")
	else
		ind = Array.new(1)
		y   = p_t(ind)
		print(String(x) + "％値 = " + String(y) + "  収束 " + String(ind[0]) + " (自由度 = " + String($dof) + ")\n")
	end
end

  6. Python 

    # -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *

############################################
# Γ（ｘ）の計算（ガンマ関数，近似式）
#      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

################################################
# ｔ分布の計算（P(X = tt), P(X < tt)）
# （自由度が∞の時の値はN(0,1)を利用して下さい）
#      dd : P(X = tt)
#      df : 自由度
#      return : P(X < tt)
################################################

def t(tt, dd, df) :

	if tt < 0.0 :
		sign = -1.0
	else :
		sign = 1.0
	if abs(tt) < 1.0e-10 :
		tt = sign * 1.0e-10
	t2 = tt * tt
	x  = t2 / (t2 + df)

	if df%2 != 0 :
		u  = sqrt(x*(1.0-x)) / pi
		p  = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi
		ia = 1

	else :
		u  = sqrt(x) * (1.0 - x) / 2.0
		p  = sqrt(x)
		ia = 2

	if ia != df :
		for i1 in range(ia, df-1, 2) :
			p += 2.0 * u / i1
			u *= ((1.0 + i1) / i1 * (1.0 - x))

	dd[0] = u / abs(tt)
	pp    = 0.5 + 0.5 * sign * p

	return pp

############################################
# 二分法による非線形方程式(f(x)=0)の解
#      fn : f(x)を計算する関数名
#      x1,x2 : 初期値
#      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
#      eps2 : 終了条件２（｜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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
#      eps2 : 終了条件２（｜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

############################
# ｔ分布の計算
#      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(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)（関数値, t 分布）
#######################################
def t_f(x) :
	y = np.empty(1, np.float)
	return t(x, y, dof) - 1.0 + p

#################################
# P(X = x)（関数の微分, t 分布）
#################################
def t_df(x) :
	y = np.empty(1, np.float)
	z = t(x, y, dof)
	return y[0]

#################################################
# ｔ分布のｐ％値（P(X > u) = 0.01p）
# （自由度が∞の時の値はN(0,1)を利用して下さい）
#      ind : >= 0 : normal（収束回数）
#            = -1 : 収束しなかった
#################################################
def p_t(ind) :

	tt  = 0.0
	pis = sqrt(pi)
	df  = float(dof)
	df2 = 0.5 * df
			# 自由度＝１
	if dof == 1 :
		tt = tan(pi*(0.5-p))

	else :
			# 自由度＝２
		if dof == 2 :
			if p > 0.5 :
				c = -1.0
			else :
				c = 1.0
			p2  = (1.0 - 2.0 * p)
			p2 *= p2
			tt  = c * sqrt(2.0 * p2 / (1.0 - p2))
			# 自由度＞２
		else :

			yq = p_normal(ind)   # 初期値計算のため

			if ind[0] >= 0 :

				x  = 1.0 - 1.0 / (4.0 * df)
				e  = x * x - yq * yq / (2.0 * df)

				if e > 0.5 :
					t0 = yq / sqrt(e)
				else :
					x  = sqrt(df) / (pis * p * df * gamma(df2, ind) / gamma(df2+0.5, ind))
					t0 = x ** (1.0/df)
					# ニュートン法
				tt = newton(t_f, t_df, t0, 1.0e-6, 1.0e-10, 100, ind)

	return tt

			# 密度関数と分布関数の値
s   = input("自由度は？ ")
dof = int(s)
print("目的とする結果は？ ")
print("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）")
s  = input("     =1 : ｐ％値（ 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  = t(x, z, dof)
		print("P(X = " + str(x) + ") = " + str(z[0]) + ",  P( X < " + str(x) + ") = " + str(y) + "(自由度 = " + str(dof) + ")")
			# グラフ出力
	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 = t(x, z, dof)
			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(dof) + ")")
	elif (1.0-p)< 1.0e-7 :
		print(str(x) + "％値 = -∞ (自由度 = " + str(dof) + ")")
	else :
		ind = np.empty(1, np.int)
		y   = p_t(ind)
		print(str(x) + "％値 = " + str(y) + "  収束 " + str(ind[0]) + " (自由度 = " + str(dof) + ")")

  7. C# 

    /****************************/
/* ｔ分布の計算             */
/*      coded by Y.Suganuma */
/****************************/
using System;
using System.IO;

class Program
{
	static void Main()
	{
		Test1 ts = new Test1();
	}
}

class Test1
{
	double p;     // α％値を計算するとき時α/100を設定
	int dof;      // 自由度

	public Test1()
	{
		Console.Write("自由度は？ ");
		dof = int.Parse(Console.ReadLine());
		Console.WriteLine("目的とする結果は？ ");
		Console.WriteLine("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）");
		Console.Write("     =1 : ｐ％値（ 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 = t(x, ref z, dof);
				Console.WriteLine("P(X = " + x + ") = " + z + ",  P( X < " + x + ") = " + y +
                                  " (自由度 = " + dof + ")");
			}
					// グラフ出力
			else {
				Console.Write("   密度関数のファイル名は？ ");
				String file1 = Console.ReadLine();
				Console.Write("   分布関数のファイル名は？ ");
				String file2 = Console.ReadLine();
				Console.Write("   データの下限は？ ");
				double from = double.Parse(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 = from; x < to+0.5*h; x += h) {
					double y = t(x, ref z, dof);
					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 + "％値 = ∞ (自由度 = " + dof + ")");
			else if ((1.0-p) < 1.0e-7)
				Console.WriteLine(x + "％値 = -∞ (自由度 = " + dof + ")");
			else {
				int sw1  = 0;
				double y = p_t(ref sw1);
				Console.WriteLine(x + "％値 = " + y + "  sw " + sw1 + " (自由度 = " + dof + ")");
			}
		}
	}

	/**************************************************/
	/* ｔ分布の計算（P(X = tt), P(X < tt)）           */
	/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
	/*      dd : P(X = tt)                            */
	/*      df : 自由度                               */
	/*      return : P(X < tt)                        */
	/**************************************************/
	double t(double tt, ref double dd, int df)
	{
		double q, u;
		int ia = 0;

		double sign = (tt < 0.0) ? -1.0 : 1.0;
		if (Math.Abs(tt) < 1.0e-10)
			tt = sign * 1.0e-10;
		double t2 = tt * tt;
		double x  = t2 / (t2 + df);

		if(df%2 != 0) {
			u  = Math.Sqrt(x*(1.0-x)) / Math.PI;
			q  = 1.0 - 2.0 * Math.Atan2(Math.Sqrt(1.0-x), Math.Sqrt(x)) / Math.PI;
			ia = 1;
		}

		else {
			u  = Math.Sqrt(x) * (1.0 - x) / 2.0;
			q  = Math.Sqrt(x);
			ia = 2;
		}

		if (ia != df) {
			for (int i1 = ia; i1 <= df-2; i1 += 2) {
				q += 2.0 * u / i1;
				u *= (1.0 + i1) / i1 * (1.0 - x);
			}
		}

		dd        = u / Math.Abs(tt);
		double qq = 0.5 + 0.5 * sign * q;

		return qq;
	}

	/**************************************************/
	/* ｔ分布のｐ％値（P(X > u) = 0.01p）             */
	/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
	/*      ind : >= 0 : normal（収束回数）           */
	/*            = -1 : 収束しなかった               */
	/**************************************************/
	double p_t(ref int ind)
	{
		double pis = Math.Sqrt(Math.PI);
		double df  = (double)dof;
		double df2 = 0.5 * df;
		double tt  = 0.0;
					// 自由度＝１
		if (dof == 1)
			tt = Math.Tan(Math.PI*(0.5-p));

		else {
					// 自由度＝２
			if (dof == 2) {
				double c   = (p > 0.5) ? -1.0 : 1.0;
				double p2  = (1.0 - 2.0 * p);
				p2        *= p2;
				tt         = c * Math.Sqrt(2.0 * p2 / (1.0 - p2));
			}
					// 自由度＞２
			else {

				double yq = p_normal(ref ind);   // 初期値計算のため

				if (ind >= 0) {

					double x  = 1.0 - 1.0 / (4.0 * df);
					double e  = x * x - yq * yq / (2.0 * df);

					double t0 = 0.0;
					if (e > 0.5)
						t0 = yq / Math.Sqrt(e);
					else {
						x  = Math.Sqrt(df) / (pis * p * df * gamma(df2, ref ind) / gamma(df2+0.5, ref ind));
						t0 = Math.Pow(x, 1.0/df);
					}
						// ニュートン法
					tt = newton(t0, 1.0e-6, 1.0e-10, 100, ref ind, t_f, t_df);
				}
			}
		}

		return tt;
	}

	/****************************************/
	/* Γ（ｘ）の計算（ガンマ関数，近似式） */
	/*      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 t_f(double x)
	{
		double y = 0.0;
		return t(x, ref y, dof) - 1.0 + p;
	}

	/**************************/
	/* P(X = x)（関数の微分） */
	/**************************/
	double t_df(double x)
	{
		double y = 0.0;
		double z = t(x, ref y, dof);

		return y;
	}

	/*****************************************************/
	/* Newton法による非線形方程式(f(x)=0)の解            */
	/*      x1 : 初期値                                  */
	/*      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   */
	/*      eps2 : 終了条件２（｜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(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
	/*      eps2 : 終了条件２（｜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 

    '**************************'
' ｔ分布の計算             '
'      coded by Y.Suganuma '
'**************************'
Imports System.IO

Module Test

	Dim p As Double     ' α％値を計算するとき時α/100を設定
	Dim dof As Integer     ' 自由度

	Sub Main()

		Console.Write("自由度は？ ")
		dof = Integer.Parse(Console.ReadLine())
		Console.WriteLine("目的とする結果は？ ")
		Console.WriteLine("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）")
		Console.Write("     =1 : ｐ％値（ 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 = t(x, z, dof)
				Console.WriteLine("P(X = " & x & ") = " & z & ",  P( X < " & x & ") = " & y &
                                  " (自由度 = " & dof & ")")
					' グラフ出力
			Else
				Console.Write("   密度関数のファイル名は？ ")
				Dim file1 As String = Console.ReadLine().Trim()
				Console.Write("   分布関数のファイル名は？ ")
				Dim file2 As String = Console.ReadLine().Trim()
				Console.Write("   データの下限は？ ")
				Dim from1 As Double = Double.Parse(Console.ReadLine())
				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 = from1
				Do While x < to1+0.5*h
					Dim y As Double = t(x, z, dof)
					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 & "％値 = ∞ (自由度 = " & dof & ")")
			ElseIf (1.0-p) < 1.0e-7
				Console.WriteLine(x & "％値 = -∞ (自由度 = " & dof & ")")
			Else
						' P(X > x) - 1 + p（関数値）（ラムダ式）
				Dim t_f = Function(v1) As Double
					Dim v2 As Double = 0.0
					Return t(v1, v2, dof) - 1.0 + p
				End Function
						' P(X = x)（関数の微分）（ラムダ式）
				Dim t_df = Function(v1) As Double
					Dim v2 As Double = 0.0
					Dim v3 As Double = t(v1, v2, dof)
					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_t(sw1, t_f, t_df, normal_f)
				Console.WriteLine(x & "％値 = " & y & "  sw " & sw1 & " (自由度 = " & dof & ")")
			End If
		End If

	End Sub

	'************************************************'
	' ｔ分布の計算（P(X = tt), P(X < tt)）           '
	' （自由度が∞の時の値はN(0,1)を利用して下さい） '
	'      dd : P(X = tt)                            '
	'      df : 自由度                               '
	'      return : P(X < tt)                        '
	'************************************************'
	Function t(tt As Double, ByRef dd As Double, df As Integer)

		Dim q As Double
		Dim u As Double
		Dim ia As Integer = 0

		Dim sign As Double
		If tt < 0.0
			sign = -1.0
		Else
			sign = 1.0
		End If
		If Math.Abs(tt) < 1.0e-10
			tt = sign * 1.0e-10
		End If
		Dim t2 As Double = tt * tt
		Dim x As Double  = t2 / (t2 + df)

		If (df Mod 2) <> 0
			u  = Math.Sqrt(x*(1.0-x)) / Math.PI
			q  = 1.0 - 2.0 * Math.Atan2(Math.Sqrt(1.0-x), Math.Sqrt(x)) / Math.PI
			ia = 1

		Else
			u  = Math.Sqrt(x) * (1.0 - x) / 2.0
			q  = Math.Sqrt(x)
			ia = 2
		End If

		If ia <> df
			For i1 As Integer = ia To df-2 Step 2
				q += 2.0 * u / i1
				u *= (1.0 + i1) / i1 * (1.0 - x)
			Next
		End If

		dd               = u / Math.Abs(tt)
		Dim qq As Double = 0.5 + 0.5 * sign * q

		Return qq

	End Function

	'************************************************'
	' ｔ分布のｐ％値（P(X > u) = 0.01p）             '
	' （自由度が∞の時の値はN(0,1)を利用して下さい） '
	'      ind : >= 0 : normal（収束回数）           '
	'            = -1 : 収束しなかった               '
	'************************************************'
	Function p_t(ByRef ind As Integer, t_f As Func(Of Double, Double),
	             t_df As Func(Of Double, Double), normal_f As Func(Of Double, Double))

		Dim pis As Double = Math.Sqrt(Math.PI)
		Dim df As Double  = dof
		Dim df2 As Double = 0.5 * df
		Dim tt As Double  = 0.0
					' 自由度＝１
		If dof = 1
			tt = Math.Tan(Math.PI*(0.5-p))

		Else
					' 自由度＝２
			If dof = 2
				Dim c As Double
				If p > 0.5
					c = -1.0
				Else
					c = 1.0
				End If
				Dim p2 As Double = (1.0 - 2.0 * p)
				p2              *= p2
				tt               = c * Math.Sqrt(2.0 * p2 / (1.0 - p2))
					' 自由度＞２
			Else

				Dim yq As Double = p_normal(ind, normal_f)   ' 初期値計算のため

				If ind >= 0

					Dim x As Double  = 1.0 - 1.0 / (4.0 * df)
					Dim e As Double  = x * x - yq * yq / (2.0 * df)

					Dim t0 As Double = 0.0
					If e > 0.5
						t0 = yq / Math.Sqrt(e)
					Else
						x  = Math.Sqrt(df) / (pis * p * df * gamma(df2, ind) / gamma(df2+0.5, ind))
						t0 = Math.Pow(x, 1.0/df)
					End If
						' ニュートン法
					tt = newton(t0, 1.0e-6, 1.0e-10, 100, ind, t_f, t_df)
				End If
			End If
		End If

		Return tt

	End Function

	'**************************************'
	' Γ（ｘ）の計算（ガンマ関数，近似式） '
	'      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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   '
	'      eps2 : 終了条件２（｜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(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   '
	'      eps2 : 終了条件２（｜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

情報学部 菅沼ホーム 目次 索引