| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/****************************/
/* χ2分布の計算 */
/* 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 chi(double, double *, int);
double p_chi(int *);
double newton(double(*)(double), double(*)(double), double, double,
double, int, int *);
double chi_f(double);
double chi_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 = chi(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 = chi(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％値 = 0 (自由度 = %d)\n", x, dof);
else {
y = p_chi(&sw);
printf("%f％値 = %f sw %d (自由度 = %d)\n", x, y, sw, dof);
}
}
return 0;
}
/********************************************/
/* χ２乗分布の計算（P(X = ch), P(X < ch)） */
/* dd : P(X = ch) */
/* df : 自由度 */
/* return : P(X < ch) */
/********************************************/
#include <math.h>
double chi(double ch, double *dd, int df)
{
double pi = 4.0 * atan(1.0);
double chs, pis, pp, x, y, u;
int ia, i1;
if (ch < 1.0e-10)
ch = 1.0e-10;
pis = sqrt(pi);
chs = sqrt(ch);
x = 0.5 * ch;
if(df%2 != 0) {
u = sqrt(x) * exp(-x) / pis;
pp = 2.0 * (normal(chs, &y) - 0.5);
ia = 1;
}
else {
u = x * exp(-x);
pp = 1.0 - exp(-x);
ia = 2;
}
if (ia != df) {
for (i1 = ia; i1 <= df-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= (ch / i1);
}
}
*dd = u / ch;
return pp;
}
/******************************************/
/* χ２乗分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >=0 : normal（収束回数） */
/* =-1 : 収束しなかった */
/******************************************/
#include <math.h>
double p_chi(int *ind)
{
double po, x, xx = 0.0, x0, w;
// 自由度＝１（正規分布を使用）
if (dof == 1) {
po = p;
p *= 0.5;
x = p_normal(ind);
xx = x * x;
p = po;
}
else {
// 自由度＝２
if (dof == 2)
xx = -2.0 * log(p);
// 自由度＞２
else {
x = p_normal(ind); // 初期値計算のため
if (*ind >= 0) {
w = 2.0 / (9.0 * dof);
x = 1.0 - w + x * sqrt(w);
x0 = pow(x, 3.0) * dof; // ニュートン法の初期値
// ニュートン法
xx = newton(chi_f, chi_df, x0, 1.0e-6, 1.0e-10, 100, ind);
}
}
}
return xx;
}
/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
double chi_f(double x)
{
double y;
return chi(x, &y, dof) - 1.0 + p;
}
/**************************/
/* P(X = x)（関数の微分） */
/**************************/
double chi_df(double x)
{
double y,z;
z = chi(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分布の計算 */
/* 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, 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.chi(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(" データの上限は？ ");
to = Double.parseDouble(in.readLine());
System.out.print(" 刻み幅は？ ");
h = Double.parseDouble(in.readLine());
// データ取得
ArrayList <Double> x1 = new ArrayList <Double> ();
ArrayList <Double> y1 = new ArrayList <Double> ();
ArrayList <Double> y2 = new ArrayList <Double> ();
for (x = 0.0; x < to+0.5*h; x += h) {
y = App.chi(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 + "％値 = 0 (自由度 = " + dof + ")");
else {
y = App.p_chi(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] = "密度関数（χ2分布 自由度 = " + 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 = 0.0;
double x_step = (x_max - x_min) / 10;
int x_p = 0;
boolean ok = true;
int k = 0;
while (ok) {
if (x_step < 1.0) {
x_step *= 10;
k++;
}
else if (x_step >= 10.0) {
x_step /= 10;
k--;
}
else {
ok = false;
if (x_step-(int)x_step > 1.0e-5)
x_step = (int)x_step + 1;
else
x_step = (int)x_step;
if (k != 0) {
x_step = x_step * Math.pow(10, -k);
if (k > 0)
x_p = k;
}
double t = x_min;
while (t < x_max-0.001*x_step)
t += x_step;
x_max = t;
}
}
x_scale[0] = x_min; // 最小値
x_scale[1] = x_max; // 最大値
x_scale[2] = x_step; // 刻み幅
// ｙ軸目盛り
y_scale1 = new double[3];
double y_max = 1.0;
double y_step = 0.2;
int y_p1 = 1;
if (dof > 1) {
y_max = 0.0;
for (int i1 = 0; i1 < n; i1++) {
if ((y1.get(i1)).doubleValue() > y_max)
y_max = (y1.get(i1)).doubleValue();
}
y_step = y_max / 5;
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] = "分布関数（χ2分布 自由度 = " + 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)）の計算（χ2分布）
case 1:
y = App.chi(x, w, Test.dof) - 1.0 + Test.p;
break;
// 関数の微分の計算（χ2分布）
case 2:
y = App.chi(x, w, Test.dof);
y = w[0];
break;
}
return y;
}
}
/************************/
/* 科学技術系算用の手法 */
/************************/
class App {
/********************************************/
/* χ２乗分布の計算（P(X = ch), P(X < ch)） */
/* dd : P(X = ch) */
/* df : 自由度 */
/* return : P(X < ch) */
/********************************************/
static double chi(double ch, double dd[], int df)
{
double pi = Math.PI;
double chs, pis, pp, x, u, y[] = new double [1];
int ia, i1;
if (ch < 1.0e-10)
ch = 1.0e-10;
pis = Math.sqrt(pi);
chs = Math.sqrt(ch);
x = 0.5 * ch;
if(df%2 != 0) {
u = Math.sqrt(x) * Math.exp(-x) / pis;
pp = 2.0 * (normal(chs, y) - 0.5);
ia = 1;
}
else {
u = x * Math.exp(-x);
pp = 1.0 - Math.exp(-x);
ia = 2;
}
if (ia != df) {
for (i1 = ia; i1 <= df-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= (ch / i1);
}
}
dd[0] = u / ch;
return pp;
}
/******************************************/
/* χ２乗分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >=0 : normal（収束回数） */
/* =-1 : 収束しなかった */
/******************************************/
static double p_chi(int ind[])
{
double po, x, xx = 0.0, x0, w;
// 自由度＝１（正規分布を使用）
if (Test.dof == 1) {
po = Test.p;
Test.p *= 0.5;
x = p_normal(ind);
xx = x * x;
Test.p = po;
}
else {
// 自由度＝２
if (Test.dof == 2)
xx = -2.0 * Math.log(Test.p);
// 自由度＞２
else {
x = p_normal(ind); // 初期値計算のため
if (ind[0] >= 0) {
// ニュートン法
w = 2.0 / (9.0 * Test.dof);
x = 1.0 - w + x * Math.sqrt(w);
x0 = Math.pow(x, 3.0) * Test.dof;
Kansu kn1 = new Kansu(1);
Kansu kn2 = new Kansu(2);
xx = newton(x0, 1.0e-6, 1.0e-10, 100, ind, kn1, kn2);
}
}
}
return xx;
}
/*************************************************/
/* 標準正規分布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);
}
}
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>χ<SUP>2</SUP>分布</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 = chi(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 = chi(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,χ2分布（密度関数）,x,f(x),1,グラフ１," + from + "," + to;
let sp1 = Math.floor((to - from) / 5 * 10);
if (sp1 % 5 > 0)
sp1 = sp1 + 5 - sp1 % 5;
sp1 /= 10;
let sp2 = Math.floor(mx / 5 * 100);
if (sp2 % 5 > 0)
sp2 = sp2 + 5 - sp2 % 5;
sp2 /= 100;
mx = sp2 * 5;
gp1 = gp1 + "," + sp1 + ",1,0.0," + mx + "," + sp2 + ",2," + n;
for (let i1 = 0; i1 < n; i1++)
gp1 = gp1 + "," + Math.round(100 * x1[i1]) / 100;
for (let i1 = 0; i1 < n; i1++)
gp1 = gp1 + "," + Math.round(1000 * y1[i1]) / 1000;
gp1 = gp1 + ",1,0";
let str = "graph_js.htm?gp=" + gp1;
open(str, "density", "width=950, height=700");
let gp2 = "2,χ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 + "％値 = 0";
document.getElementById("xx").value = str;
}
else {
let y = p_chi(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)）の計算（χ2分布）
case 1:
y = chi(x, w, dof) - 1.0 + p;
break;
// 関数の微分の計算（χ2分布）
case 2:
y = chi(x, w, dof);
y = w[0];
break;
}
return y;
}
/********************************************/
/* χ２乗分布の計算（P(X = ch), P(X < ch)） */
/* dd : P(X = ch) */
/* df : 自由度 */
/* return : P(X < ch) */
/********************************************/
function chi(ch, dd, df)
{
let pi = Math.PI;
let chs;
let pis;
let pp;
let x;
let u;
let y = new Array();
let ia;
let i1;
if (ch < 1.0e-10)
ch = 1.0e-10;
pis = Math.sqrt(pi);
chs = Math.sqrt(ch);
x = 0.5 * ch;
if(df%2 != 0) {
u = Math.sqrt(x) * Math.exp(-x) / pis;
pp = 2.0 * (normal(chs, y) - 0.5);
ia = 1;
}
else {
u = x * Math.exp(-x);
pp = 1.0 - Math.exp(-x);
ia = 2;
}
if (ia != df) {
for (i1 = ia; i1 <= df-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= (ch / i1);
}
}
dd[0] = u / ch;
return pp;
}
/******************************************/
/* χ２乗分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >=0 : normal（収束回数） */
/* =-1 : 収束しなかった */
/******************************************/
function p_chi(ind)
{
let po;
let x;
let xx = 0.0;
let x0;
let w;
// 自由度＝１（正規分布を使用）
if (dof == 1) {
po = p;
p *= 0.5;
x = p_normal(ind);
xx = x * x;
p = po;
}
else {
// 自由度＝２
if (dof == 2)
xx = -2.0 * Math.log(p);
// 自由度＞２
else {
x = p_normal(ind); // 初期値計算のため
if (ind[0] >= 0) {
// ニュートン法
w = 2.0 / (9.0 * dof);
x = 1.0 - w + x * Math.sqrt(w);
x0 = Math.pow(x, 3.0) * dof;
xx = newton(x0, 1.0e-10, 1.0e-15, 100, ind);
}
}
}
return xx;
}
/*************************************************/
/* 標準正規分布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-6, 1.0e-10, 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>χ<SUP>2</SUP>分布</B></H2>
自由度：<INPUT ID="dof" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="5">
<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" DISABLED STYLE="font-size: 100%;" TYPE="text" SIZE="5" VALUE="0.0"> </SPAN>
<SPAN ID="to_t" STYLE="display: none">上限：<INPUT ID="to" STYLE="font-size: 100%;" TYPE="text" SIZE="5.0" VALUE="15.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>
<!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>
<?php
/****************************/
/* χ2分布の計算 */
/* 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;
}
/********************************************/
/* χ２乗分布の計算（P(X = ch), P(X < ch)） */
/* dd : P(X = ch) */
/* df : 自由度 */
/* return : P(X < ch) */
/********************************************/
function chi($ch, &$dd, $df)
{
$pi = 4.0 * atan(1.0);
if ($ch < 1.0e-10)
$ch = 1.0e-10;
$pis = sqrt($pi);
$chs = sqrt($ch);
$x = 0.5 * $ch;
if($df%2 != 0) {
$u = sqrt($x) * exp(-$x) / $pis;
$pp = 2.0 * (normal($chs, $y) - 0.5);
$ia = 1;
}
else {
$u = $x * exp(-$x);
$pp = 1.0 - exp(-$x);
$ia = 2;
}
if ($ia != $df) {
for ($i1 = $ia; $i1 <= $df-2; $i1 += 2) {
$pp -= 2.0 * $u / $i1;
$u *= ($ch / $i1);
}
}
$dd = $u / $ch;
return $pp;
}
/******************************************/
/* χ２乗分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >=0 : normal（収束回数） */
/* =-1 : 収束しなかった */
/******************************************/
function p_chi(&$ind)
{
global $p, $dof;
$xx = 0.0;
// 自由度＝１（正規分布を使用）
if ($dof == 1) {
$po = $p;
$p *= 0.5;
$x = p_normal($ind);
$xx = $x * $x;
$p = $po;
}
else {
// 自由度＝２
if ($dof == 2)
$xx = -2.0 * log($p);
// 自由度＞２
else {
$x = p_normal($ind); // 初期値計算のため
if ($ind >= 0) {
$w = 2.0 / (9.0 * $dof);
$x = 1.0 - $w + $x * sqrt($w);
$x0 = pow($x, 3.0) * $dof; // ニュートン法の初期値
// ニュートン法
$xx = newton("chi_f", "chi_df", $x0, 1.0e-6, 1.0e-10, 100, $ind);
}
}
}
return $xx;
}
/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
function chi_f($x)
{
global $p, $dof;
return chi($x, $y, $dof) - 1.0 + $p;
}
/**************************/
/* P(X = x)（関数の微分） */
/**************************/
function chi_df($x)
{
global $dof;
$z = chi($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 = chi($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,"w");
$out2 = fopen($file2,"w");
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 = chi($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％値 = 0 (自由度 = %d)\n", $x, $dof);
else {
$y = p_chi($sw);
printf("%f％値 = %f sw %d (自由度 = %d)\n", $x, $y, $sw, $dof);
}
}
?>
############################
# χ2分布の計算
# coded by Y.Suganuma
############################
################################################
# 標準正規分布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 = ch), P(X < ch)）
# dd : P(X = ch)
# df : 自由度
# return : P(X < ch)
##########################################
def chi(ch, dd, df)
if ch < 1.0e-10
ch = 1.0e-10
end
pis = Math.sqrt(Math::PI)
chs = Math.sqrt(ch)
x = 0.5 * ch
y = Array.new(1)
if df%2 != 0
u = Math.sqrt(x) * Math.exp(-x) / pis
pp = 2.0 * (normal(chs, y) - 0.5)
ia = 1
else
u = x * Math.exp(-x)
pp = 1.0 - Math.exp(-x)
ia = 2
end
if ia != df
i1 = ia
while i1 < df-1
pp -= 2.0 * u / i1
u *= (ch / i1)
i1 += 2
end
end
dd[0] = u / ch
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）
# ind : >=0 : normal（収束回数）
# =-1 : 収束しなかった
##########################################
def p_chi(ind)
chi_f = Proc.new { |sw, x|
y = Array.new(1)
z = chi(x, y, $dof)
if sw == 0
z - 1.0 + $p
else
z = y[0]
end
}
xx = 0.0
# 自由度＝１（正規分布を使用）
if $dof == 1
po = $p
$p *= 0.5
x = p_normal(ind)
xx = x * x
$p = po
else
# 自由度＝２
if $dof == 2
xx = -2.0 * Math.log($p)
# 自由度＞２
else
x = p_normal(ind) # 初期値計算のため
if ind[0] >= 0
w = 2.0 / (9.0 * $dof)
x = 1.0 - w + x * Math.sqrt(w)
x0 = (x ** 3.0) * $dof # ニュートン法の初期値
xx = newton(x0, 1.0e-6, 1.0e-10, 100, ind, &chi_f)
end
end
end
return xx
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 = chi(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 = chi(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) + "％値 = 0 (自由度 = " + String($dof) + ")\n")
else
ind = Array.new(1)
y = p_chi(ind)
print(String(x) + "％値 = " + String(y) + " 収束 " + String(ind[0]) + " (自由度 = " + String($dof) + ")\n")
end
end
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
################################################
# 標準正規分布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 = ch), P(X < ch)）
# dd : P(X = ch)
# df : 自由度
# return : P(X < ch)
##########################################
def chi(ch, dd, df) :
if ch < 1.0e-10 :
ch = 1.0e-10
pis = sqrt(pi)
chs = sqrt(ch)
x = 0.5 * ch
y = np.empty(1, np.float)
if df%2 != 0 :
u = sqrt(x) * exp(-x) / pis
pp = 2.0 * (normal(chs, y) - 0.5)
ia = 1
else :
u = x * exp(-x)
pp = 1.0 - exp(-x)
ia = 2
if ia != df :
for i1 in range(ia, df-1, 2) :
pp -= 2.0 * u / i1
u *= (ch / i1)
dd[0] = u / ch
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
############################
# χ2分布の計算
# 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)（関数値, χ２乗分布）
##########################################
def chi_f(x) :
y = np.empty(1, np.float)
return (chi(x, y, dof) - 1.0 + p)
####################################
# P(X = x)（関数の微分, χ２乗分布）
####################################
def chi_df(x) :
y = np.empty(1, np.float)
z = chi(x, y, dof)
return y[0]
##########################################
# χ２乗分布のｐ％値（P(X > u) = 0.01p）
# ind : >=0 : normal（収束回数）
# =-1 : 収束しなかった
##########################################
def p_chi(ind) :
xx = 0.0
# 自由度＝１（正規分布を使用）
if dof == 1 :
global p # この宣言がないとエラー
po = p
p *= 0.5
x = p_normal(ind)
xx = x * x
p = po
else :
# 自由度＝２
if dof == 2 :
xx = -2.0 * log(p)
# 自由度＞２
else :
x = p_normal(ind) # 初期値計算のため
if ind[0] >= 0 :
w = 2.0 / (9.0 * dof)
x = 1.0 - w + x * sqrt(w)
x0 = (x ** 3.0) * dof # ニュートン法の初期値
xx = newton(chi_f, chi_df, x0, 1.0e-6, 1.0e-10, 100, ind)
return xx
# 密度関数と分布関数の値
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 = chi(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 = chi(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) + "％値 = 0 (自由度 = " + str(dof) + ")")
else :
ind = np.empty(1, np.int)
y = p_chi(ind)
print(str(x) + "％値 = " + str(y) + " 収束 " + str(ind[0]) + " (自由度 = " + str(dof) + ")")
/****************************/
/* χ2分布の計算 */
/* 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 = chi(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 to = double.Parse(Console.ReadLine());
Console.Write(" 刻み幅は？ ");
double h = double.Parse(Console.ReadLine());
StreamWriter out1 = new StreamWriter(file1);
StreamWriter out2 = new StreamWriter(file2);
// データ取得
for (double x = 0.0; x < to+0.5*h; x += h) {
double y = chi(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 + "％値 = 0 (自由度 = " + dof + ")");
else {
int sw1 = 0;
double y = p_chi(ref sw1);
Console.WriteLine(x + "％値 = " + y + " sw " + sw1 + " (自由度 = " + dof + ")");
}
}
}
/********************************************/
/* χ２乗分布の計算（P(X = ch), P(X < ch)） */
/* dd : P(X = ch) */
/* df : 自由度 */
/* return : P(X < ch) */
/********************************************/
double chi(double ch, ref double dd, int df)
{
double pp, y = 0.0, u;
int ia = 0;
if (ch < 1.0e-10)
ch = 1.0e-10;
double pis = Math.Sqrt(Math.PI);
double chs = Math.Sqrt(ch);
double x = 0.5 * ch;
if(df%2 != 0) {
u = Math.Sqrt(x) * Math.Exp(-x) / pis;
pp = 2.0 * (normal(chs, ref y) - 0.5);
ia = 1;
}
else {
u = x * Math.Exp(-x);
pp = 1.0 - Math.Exp(-x);
ia = 2;
}
if (ia != df) {
for (int i1 = ia; i1 <= df-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= (ch / i1);
}
}
dd = u / ch;
return pp;
}
/******************************************/
/* χ２乗分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >=0 : normal（収束回数） */
/* =-1 : 収束しなかった */
/******************************************/
double p_chi(ref int ind)
{
double xx = 0.0;
// 自由度＝１（正規分布を使用）
if (dof == 1) {
double po = p;
p *= 0.5;
double x = p_normal(ref ind);
xx = x * x;
p = po;
}
else {
// 自由度＝２
if (dof == 2)
xx = -2.0 * Math.Log(p);
// 自由度＞２
else {
double x = p_normal(ref ind); // 初期値計算のため
if (ind >= 0) {
double w = 2.0 / (9.0 * dof);
x = 1.0 - w + x * Math.Sqrt(w);
double x0 = Math.Pow(x, 3.0) * dof; // ニュートン法の初期値
// ニュートン法
xx = newton(x0, 1.0e-6, 1.0e-10, 100, ref ind, chi_f, chi_df);
}
}
}
return xx;
}
/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
double chi_f(double x)
{
double y = 0.0;
return chi(x, ref y, dof) - 1.0 + p;
}
/**************************/
/* P(X = x)（関数の微分） */
/**************************/
double chi_df(double x)
{
double y = 0.0;
double z = chi(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;
}
}
'**************************'
' χ2分布の計算 '
' 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 = chi(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 to1 As Double = Double.Parse(Console.ReadLine())
Console.Write(" 刻み幅は？ ")
Dim h As Double = Double.Parse(Console.ReadLine())
Dim out1 As StreamWriter = new StreamWriter(file1)
Dim out2 As StreamWriter = new StreamWriter(file2)
' データ取得
Dim x As Double = 0.0
Do While x < to1+0.5*h
Dim y As Double = chi(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 & "％値 = 0 (自由度 = " & dof & ")")
Else
' P(X > x) - 1 + p（関数値）（ラムダ式）
Dim chi_f = Function(v1) As Double
Dim v2 As Double = 0.0
Return chi(v1, v2, dof) - 1.0 + p
End Function
' P(X = x)（関数の微分）（ラムダ式）
Dim chi_df = Function(v1) As Double
Dim v2 As Double = 0.0
Dim v3 As Double = chi(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_chi(sw1, chi_f, chi_df, normal_f)
Console.WriteLine(x & "％値 = " & y & " sw " & sw1 & " (自由度 = " & dof & ")")
End If
End If
End Sub
'******************************************'
' χ２乗分布の計算（P(X = ch), P(X < ch)） '
' dd : P(X = ch) '
' df : 自由度 '
' return : P(X < ch) '
'******************************************'
Function chi(ch As Double, ByRef dd As Double, df As Integer)
Dim pp As Double
Dim y As Double = 0.0
Dim u As Double
Dim ia As Integer = 0
If ch < 1.0e-10
ch = 1.0e-10
End If
Dim pis As Double = Math.Sqrt(Math.PI)
Dim chs As Double = Math.Sqrt(ch)
Dim x As Double = 0.5 * ch
If (df MOd 2) <> 0
u = Math.Sqrt(x) * Math.Exp(-x) / pis
pp = 2.0 * (normal(chs, y) - 0.5)
ia = 1
Else
u = x * Math.Exp(-x)
pp = 1.0 - Math.Exp(-x)
ia = 2
End If
If ia <> df
For i1 As Integer = ia To df-2 Step 2
pp -= 2.0 * u / i1
u *= (ch / i1)
Next
End If
dd = u / ch
Return pp
End Function
'****************************************'
' χ２乗分布のｐ％値（P(X > u) = 0.01p） '
' ind : >=0 : normal（収束回数） '
' =-1 : 収束しなかった '
'****************************************'
Function p_chi(ByRef ind As Integer, chi_f As Func(Of Double, Double),
chi_df As Func(Of Double, Double), normal_f As Func(Of Double, Double))
Dim xx As Double = 0.0
' 自由度＝１（正規分布を使用）
If dof = 1
Dim po As Double = p
p *= 0.5
Dim x As Double = p_normal(ind, normal_f)
xx = x * x
p = po
Else
' 自由度＝２
If dof = 2
xx = -2.0 * Math.Log(p)
' 自由度＞２
Else
Dim x As Double = p_normal(ind, normal_f) ' 初期値計算のため
If ind >= 0
Dim w As Double = 2.0 / (9.0 * dof)
x = 1.0 - w + x * Math.Sqrt(w)
Dim x0 As Double = Math.Pow(x, 3.0) * dof ' ニュートン法の初期値
' ニュートン法
xx = newton(x0, 1.0e-6, 1.0e-10, 100, ind, chi_f, chi_df)
End If
End If
End If
Return xx
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
| 情報学部 | 菅沼ホーム | 目次 | 索引 |