| 情報学部 | 菅沼ホーム | 目次 | 索引 |
#include <stdio.h>
double normal(double, double *);
double p_normal(int *);
double bisection(double(*)(double), double, double, double, double, int, int *);
double normal_f(double);
double f(double, double *, int, int);
double p_f(int *);
double gamma(double, int *);
double newton(double(*)(double), double(*)(double), double, double,
double, int, int *);
double f_f(double);
double f_df(double);
double p; // α％値を計算するとき時α/100を設定
int dof1, dof2; // 自由度
/**************************************/
/* 分散分析（一元配置法と二元配置法） */
/* method : 1 or 2 */
/* Np[0] : 因子1の水準数 */
/* [1] : 因子2の水準数 */
/* N : データ数 */
/* name[0] : 因子1の名前 */
/* [1] : 因子2の名前 */
/* a : a ％値 */
/* x : データ */
/* coded by Y.Suganuma */
/**************************************/
void aov(int method, int *Np, int N, char name[][100], int a, double ***x)
{
double *xi, *xj, **xij, xa, SP, SQ, SI, SE, VP, VQ, VI, VE, FP, FQ, FI;
int i1, i2, i3, sw;
// 一元配置法
if (method == 1) {
xi = new double [Np[0]];
for (i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (i2 = 0; i2 < N; i2++)
xi[i1] += x[i1][0][i2];
xi[i1] /= N;
}
xa = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < N; i2++)
xa += x[i1][0][i2];
}
xa /= (Np[0] * N);
SP = 0.0;
for (i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= N;
SE = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < N; i2++)
SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
}
VP = SP / (Np[0] - 1);
VE = SE / (Np[0] * (N - 1));
FP = VP / VE;
p = 0.01 * a;
dof2 = Np[0] * (N - 1);
printf("全変動： 平方和 %.2f 自由度 %d\n", SP+SE, Np[0]*N-1);
dof1 = Np[0] - 1;
printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, p_f(&sw));
printf("水準内： 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*(N-1), VE);
}
// 二元配置法
else {
xi = new double [Np[0]];
for (i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
xi[i1] += x[i1][i2][i3];
}
xi[i1] /= (Np[1] * N);
}
xj = new double [Np[1]];
for (i1 = 0; i1 < Np[1]; i1++) {
xj[i1] = 0.0;
for (i2 = 0; i2 < Np[0]; i2++) {
for (i3 = 0; i3 < N; i3++)
xj[i1] += x[i2][i1][i3];
}
xj[i1] /= (Np[0] * N);
}
xij = new double * [Np[0]];
for (i1 = 0; i1 < Np[0]; i1++) {
xij[i1] = new double [Np[1]];
for (i2 = 0; i2 < Np[1]; i2++) {
xij[i1][i2] = 0.0;
for (i3 = 0; i3 < N; i3++)
xij[i1][i2] += x[i1][i2][i3];
xij[i1][i2] /= N;
}
}
xa = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
xa += x[i1][i2][i3];
}
}
xa /= (Np[0] * Np[1] * N);
SP = 0.0;
for (i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= (Np[1] * N);
SQ = 0.0;
for (i1 = 0; i1 < Np[1]; i1++)
SQ += (xj[i1] - xa) * (xj[i1] - xa);
SQ *= (Np[0] * N);
SI = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++)
SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
}
SI *= N;
SE = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
}
}
VP = SP / (Np[0] - 1);
VQ = SQ / (Np[1] - 1);
VI = SI / ((Np[0] - 1) * (Np[1] - 1));
VE = SE / (Np[0] * Np[1] * (N - 1));
FP = VP / VE;
FQ = VQ / VE;
FI = VI / VE;
p = 0.01 * a;
dof2 = Np[0] * Np[1] * (N - 1);
printf("全変動： 平方和 %.2f 自由度 %d\n", SP+SQ+SI+SE, Np[0]*Np[1]*N-1);
dof1 = Np[0] - 1;
printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, p_f(&sw));
dof1 = Np[1] - 1;
printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", name[1], SQ, Np[1]-1, VQ, FQ, a, p_f(&sw));
dof1 = (Np[0] - 1) * (Np[1] - 1);
printf("相互作用： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", SI, (Np[0]-1)*(Np[1]-1), VI, FI, a, p_f(&sw));
printf("水準内： 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*Np[1]*(N-1), VE);
}
}
int main()
{
double ***x;
int i1, i2, i3, method, N, Np[2] = {1, 1}, a;
char name[2][100];
scanf("%d %d %d", &method, &N, &a);
if (method == 1 || method == 2) {
for (i1 = 0; i1 < method; i1++)
scanf("%s %d", name[i1], &Np[i1]);
x = new double ** [Np[0]];
for (i1 = 0; i1 < Np[0]; i1++) {
x[i1] = new double * [Np[1]];
for (i2 = 0; i2 < Np[1]; i2++)
x[i1][i2] = new double [N];
}
for (i1 = 0; i1 < Np[1]; i1++) {
for (i2 = 0; i2 < N; i2++) {
for (i3 = 0; i3 < Np[0]; i3++)
scanf("%lf", &x[i3][i1][i2]);
}
}
aov(method, Np, N, name, a, x);
}
else
printf("一元配置法，または，二元配置法だけです！\n");
return 0;
}
/****************************************/
/* ｆ分布の計算（P(X = ff), P(X < ff)） */
/* dd : P(X = ff) */
/* df1,df2 : 自由度 */
/* return : P(X < ff) */
/****************************************/
#include <math.h>
double f(double ff, double *dd, int df1, int df2)
{
double pi = 4.0 * atan(1.0);
double pp, u, x;
int ia, ib, i1;
if (ff < 1.0e-10)
ff = 1.0e-10;
x = ff * df1 / (ff * df1 + df2);
if(df1%2 == 0) {
if (df2%2 == 0) {
u = x * (1.0 - x);
pp = x;
ia = 2;
ib = 2;
}
else {
u = 0.5 * x * sqrt(1.0-x);
pp = 1.0 - sqrt(1.0-x);
ia = 2;
ib = 1;
}
}
else {
if (df2%2 == 0) {
u = 0.5 * sqrt(x) * (1.0 - x);
pp = sqrt(x);
ia = 1;
ib = 2;
}
else {
u = sqrt(x*(1.0-x)) / pi;
pp = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi;
ia = 1;
ib = 1;
}
}
if (ia != df1) {
for (i1 = ia; i1 <= df1-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= x * (i1 + ib) / i1;
}
}
if (ib != df2) {
for (i1 = ib; i1 <= df2-2; i1 += 2) {
pp += 2.0 * u / i1;
u *= (1.0 - x) * (i1 + df1) / i1;
}
}
*dd = u / ff;
return pp;
}
/****************************************/
/* ｆ分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/****************************************/
#include <math.h>
#define MAX 340
double p_f(int *ind)
{
double a, a1, b, b1, df1, df2, e, ff = 0.0, f0, x, y1, y2, yq;
int sw = 0;
/*
初期値計算の準備
while文は，大きな自由度によるガンマ関数の
オーバーフローを避けるため
*/
while (sw >= 0) {
df1 = 0.5 * (dof1 - sw);
df2 = 0.5 * dof2;
a = 2.0 / (9.0 * (dof1 - sw));
a1 = 1.0 - a;
b = 2.0 / (9.0 * dof2);
b1 = 1.0 - b;
yq = p_normal(ind);
e = b1 * b1 - b * yq * yq;
if (e > 0.8 || (dof1+dof2-sw) <= MAX)
sw = -1;
else {
sw += 1;
if ((dof1-sw) == 0)
sw = -2;
}
}
if (sw == -2)
*ind = -1;
else {
/*
f0 : 初期値
*/
if (e > 0.8) {
x = (a1 * b1 + yq * sqrt(a1*a1*b+a*e)) / e;
f0 = pow(x, 3.0);
}
else {
y1 = pow((double)dof2, df2-1.0);
y2 = pow((double)dof1, df2);
x = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
2.0 * y1 / y2 / p;
f0 = pow(x, 2.0/dof2);
}
/*
ニュートン法
*/
ff = newton(f_f, f_df, f0, 1.0e-6, 1.0e-10, 100, ind);
}
return ff;
}
/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/* ier : =0 : normal */
/* =-1 : x=-n (n=0,1,2,･･･) */
/* return : 結果 */
/****************************************/
#include <math.h>
double gamma(double x, int *ier)
{
double err, g, s, t, v, w, y;
long k;
*ier = 0;
if (x > 5.0) {
v = 1.0 / x;
s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
0.00078403922172) * v - 0.000229472093621) * v -
0.00268132716049) * v + 0.00347222222222) * v +
0.0833333333333) * v + 1.0;
g = 2.506628274631001 * exp(-x) * pow(x,x-0.5) * s;
}
else {
err = 1.0e-20;
w = x;
t = 1.0;
if (x < 1.5) {
if (x < err) {
k = (long)x;
y = (double)k - x;
if (fabs(y) < err || fabs(1.0-y) < err)
*ier = -1;
}
if (*ier == 0) {
while (w < 1.5) {
t /= w;
w += 1.0;
}
}
}
else {
if (w > 2.5) {
while (w > 2.5) {
w -= 1.0;
t *= w;
}
}
}
w -= 2.0;
g = (((((((0.0021385778 * w - 0.0034961289) * w +
0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
0.999999926;
g *= t;
}
return g;
}
/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
double f_f(double x)
{
double y;
return f(x, &y, dof1, dof2) - 1.0 + p;
}
/**************************/
/* P(X = x)（関数の微分） */
/**************************/
double f_df(double x)
{
double y, z;
z = f(x, &y, dof1, dof2);
return y;
}
/*****************************************************/
/* Newton法による非線形方程式(f(x)=0)の解 */
/* fn : f(x)を計算する関数名 */
/* dfn : f(x)の微分を計算する関数名 */
/* x0 : 初期値 */
/* eps1 : 終了条件１（｜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)の解 */
/* fn : 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;
}
/*
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 // 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 // 因子の名前 水準の数
3.1 4.7 5.1 // 各水準に対する1番目のデータ
4.1 5.6 3.7 // 各水準に対する2番目のデータ
3.3 4.3 4.5 // 各水準に対する3番目のデータ
3.9 5.9 6.0 // 各水準に対する4番目のデータ
3.7 6.1 3.9 // 各水準に対する5番目のデータ
2.4 4.2 5.4 // 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 // 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 // 1番目の因子の名前 その水準の数
品種 2 // 2番目の因子の名前 その水準の数
3 4 12 -4 -4 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ
*/
/****************************/
/* 分散分析 */
/* coded by Y.Suganuma */
/****************************/
import java.io.*;
import java.util.StringTokenizer;
public class Test {
public static void main(String args[]) throws IOException
{
double x[][][];
int i1, i2, i3, method, N, a, Np[] = new int [2];
String line, name[] = new String [2];
StringTokenizer str;
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
Np[0] = 1;
Np[1] = 1;
line = in.readLine();
str = new StringTokenizer(line, " \r\n");
method = Integer.parseInt(str.nextToken());
N = Integer.parseInt(str.nextToken());
a = Integer.parseInt(str.nextToken());
if (method == 1 || method == 2) {
for (i1 = 0; i1 < method; i1++) {
line = in.readLine();
str = new StringTokenizer(line, " \r\n");
name[i1] = str.nextToken();
Np[i1] = Integer.parseInt(str.nextToken());
}
x = new double [Np[0]][Np[1]][N];
for (i1 = 0; i1 < Np[1]; i1++) {
for (i2 = 0; i2 < N; i2++) {
line = in.readLine();
str = new StringTokenizer(line, " \r\n");
for (i3 = 0; i3 < Np[0]; i3++) {
x[i3][i1][i2] = Double.parseDouble(str.nextToken());
}
}
}
aov(method, Np, N, name, a, x);
}
else
System.out.println("一元配置法，または，二元配置法だけです！");
}
/**************************************/
/* 分散分析（一元配置法と二元配置法） */
/* method : 1 or 2 */
/* Np[0] : 因子1の水準数 */
/* [1] : 因子2の水準数 */
/* N : データ数 */
/* name[0] : 因子1の名前 */
/* [1] : 因子2の名前 */
/* a : a ％値 */
/* x : データ */
/* coded by Y.Suganuma */
/**************************************/
static void aov(int method, int Np[], int N, String name[], int a, double x[][][])
{
double xi[], xj[], xij[][], xa, SP, SQ, SI, SE, VP, VQ, VI, VE, FP, FQ, FI, p;
int i1, i2, i3, sw[] = new int [1], dof1, dof2;
// 一元配置法
if (method == 1) {
xi = new double [Np[0]];
for (i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (i2 = 0; i2 < N; i2++)
xi[i1] += x[i1][0][i2];
xi[i1] /= N;
}
xa = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < N; i2++)
xa += x[i1][0][i2];
}
xa /= (Np[0] * N);
SP = 0.0;
for (i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= N;
SE = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < N; i2++)
SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
}
VP = SP / (Np[0] - 1);
VE = SE / (Np[0] * (N - 1));
FP = VP / VE;
p = 0.01 * a;
dof2 = Np[0] * (N - 1);
System.out.printf("全変動： 平方和 %.2f 自由度 %d\n", SP+SE, Np[0]*N-1);
dof1 = Np[0] - 1;
F kn = new F(dof1, dof2, p);
System.out.printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, kn.p_f(sw));
System.out.printf("水準内： 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*(N-1), VE);
}
// 二元配置法
else {
xi = new double [Np[0]];
for (i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
xi[i1] += x[i1][i2][i3];
}
xi[i1] /= (Np[1] * N);
}
xj = new double [Np[1]];
for (i1 = 0; i1 < Np[1]; i1++) {
xj[i1] = 0.0;
for (i2 = 0; i2 < Np[0]; i2++) {
for (i3 = 0; i3 < N; i3++)
xj[i1] += x[i2][i1][i3];
}
xj[i1] /= (Np[0] * N);
}
xij = new double [Np[0]][Np[1]];
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++) {
xij[i1][i2] = 0.0;
for (i3 = 0; i3 < N; i3++)
xij[i1][i2] += x[i1][i2][i3];
xij[i1][i2] /= N;
}
}
xa = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
xa += x[i1][i2][i3];
}
}
xa /= (Np[0] * Np[1] * N);
SP = 0.0;
for (i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= (Np[1] * N);
SQ = 0.0;
for (i1 = 0; i1 < Np[1]; i1++)
SQ += (xj[i1] - xa) * (xj[i1] - xa);
SQ *= (Np[0] * N);
SI = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++)
SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
}
SI *= N;
SE = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
}
}
VP = SP / (Np[0] - 1);
VQ = SQ / (Np[1] - 1);
VI = SI / ((Np[0] - 1) * (Np[1] - 1));
VE = SE / (Np[0] * Np[1] * (N - 1));
FP = VP / VE;
FQ = VQ / VE;
FI = VI / VE;
p = 0.01 * a;
dof2 = Np[0] * Np[1] * (N - 1);
System.out.printf("全変動： 平方和 %.2f 自由度 %d\n", SP+SQ+SI+SE, Np[0]*Np[1]*N-1);
dof1 = Np[0] - 1;
F kn = new F(dof1, dof2, p);
System.out.printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", name[0], SP, Np[0]-1, VP, FP, a, kn.p_f(sw));
dof1 = Np[1] - 1;
kn = new F(dof1, dof2, p);
System.out.printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", name[1], SQ, Np[1]-1, VQ, FQ, a, kn.p_f(sw));
dof1 = (Np[0] - 1) * (Np[1] - 1);
kn = new F(dof1, dof2, p);
System.out.printf("相互作用： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", SI, (Np[0]-1)*(Np[1]-1), VI, FI, a, kn.p_f(sw));
System.out.printf("水準内： 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", SE, Np[0]*Np[1]*(N-1), VE);
}
}
}
/****************/
/* 関数値の計算 */
/****************/
class F extends Newton_Bisection_Gamma
{
int sw;
int dof1, dof2; // 自由度
double p; // α％値を計算するとき時α/100を設定
// コンストラクタ
F(int dof1, int dof2, double p)
{
this.dof1 = dof1;
this.dof2 = dof2;
this.p = p;
}
// 関数値の計算
double snx(double x)
{
double y = 0.0, w[] = new double [1];
switch (sw) {
// 関数値（f(x)）の計算（正規分布）
case 0:
y = 1.0 - p - normal(x, w);
break;
// 関数値（f(x)）の計算（Ｆ分布）
case 1:
y = f(x, w) - 1.0 + p;
break;
}
return y;
}
// 関数の微分の計算（Ｆ分布）
double dsnx(double x)
{
double y = 0.0, w[] = new double [1];
y = f(x, w);
y = w[0];
return y;
}
/****************************************/
/* ｆ分布の計算（P(X = ff), P(X < ff)） */
/* dd : P(X = ff) */
/* return : P(X < ff) */
/****************************************/
double f(double ff, double dd[])
{
double pi = Math.PI;
double pp, u, x;
int ia, ib, i1;
if (ff < 1.0e-10)
ff = 1.0e-10;
x = ff * dof1 / (ff * dof1 + dof2);
if(dof1%2 == 0) {
if (dof2%2 == 0) {
u = x * (1.0 - x);
pp = x;
ia = 2;
ib = 2;
}
else {
u = 0.5 * x * Math.sqrt(1.0-x);
pp = 1.0 - Math.sqrt(1.0-x);
ia = 2;
ib = 1;
}
}
else {
if (dof2%2 == 0) {
u = 0.5 * Math.sqrt(x) * (1.0 - x);
pp = Math.sqrt(x);
ia = 1;
ib = 2;
}
else {
u = Math.sqrt(x*(1.0-x)) / pi;
pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
ia = 1;
ib = 1;
}
}
if (ia != dof1) {
for (i1 = ia; i1 <= dof1-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= x * (i1 + ib) / i1;
}
}
if (ib != dof2) {
for (i1 = ib; i1 <= dof2-2; i1 += 2) {
pp += 2.0 * u / i1;
u *= (1.0 - x) * (i1 + dof1) / i1;
}
}
dd[0] = u / ff;
return pp;
}
/****************************************/
/* ｆ分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/****************************************/
double p_f(int ind[])
{
double a = 0.0, a1 = 0.0, b = 0.0, b1 = 0.0, df1 = 0.0, df2 = 0.0,
e = 0.0, ff = 0.0, f0, x, y1, y2, yq = 0.0;
int sww = 0, MAX = 340;
/*
初期値計算の準備
while文は，大きな自由度によるガンマ関数の
オーバーフローを避けるため
*/
while (sww >= 0) {
df1 = 0.5 * (dof1 - sww);
df2 = 0.5 * dof2;
a = 2.0 / (9.0 * (dof1 - sww));
a1 = 1.0 - a;
b = 2.0 / (9.0 * dof2);
b1 = 1.0 - b;
yq = p_normal(ind);
e = b1 * b1 - b * yq * yq;
if (e > 0.8 || (dof1+dof2-sww) <= MAX)
sww = -1;
else {
sww += 1;
if ((dof1-sww) == 0)
sww = -2;
}
}
if (sww == -2)
ind[0] = -1;
/*
f0 : 初期値
*/
else {
if (e > 0.8) {
x = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e;
f0 = Math.pow(x, 3.0);
}
else {
y1 = Math.pow((double)dof2, df2-1.0);
y2 = Math.pow((double)dof1, df2);
x = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
2.0 * y1 / y2 / p;
f0 = Math.pow(x, 2.0/dof2);
}
/*
ニュートン法
*/
sw = 1;
ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind);
}
return ff;
}
/*************************************************/
/* 標準正規分布N(0,1)の計算（P(X = x), P(X < x)）*/
/* w : P(X = x) */
/* return : P(X < x) */
/*************************************************/
double normal(double x, double w[])
{
double y, z, P;
/*
確率密度関数（定義式）
*/
w[0] = Math.exp(-0.5 * x * x) / Math.sqrt(2.0*Math.PI);
/*
確率分布関数（近似式を使用）
*/
y = 0.70710678118654 * Math.abs(x);
z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
y * 0.0000430638)))));
P = 1.0 - Math.pow(z, -16.0);
if (x < 0.0)
P = 0.5 - 0.5 * P;
else
P = 0.5 + 0.5 * P;
return P;
}
/******************************************************************/
/* 標準正規分布N(0,1)のｐ％値（P(X > u) = 0.01p）（二分法を使用） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/******************************************************************/
double p_normal(int ind[])
{
double u;
int sww[] = new int [1];
sw = 0;
u = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, sww);
ind[0] = sww[0];
return u;
}
}
abstract class Newton_Bisection_Gamma
{
abstract double snx(double x);
abstract double dsnx(double x);
/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/* ier : =0 : normal */
/* =-1 : x=-n (n=0,1,2,･･･) */
/* return : 結果 */
/****************************************/
double gamma(double x, int ier[])
{
double err, g, s, t, v, w, y;
int k;
ier[0] = 0;
if (x > 5.0) {
v = 1.0 / x;
s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
0.00078403922172) * v - 0.000229472093621) * v -
0.00268132716049) * v + 0.00347222222222) * v +
0.0833333333333) * v + 1.0;
g = 2.506628274631001 * Math.exp(-x) * Math.pow(x,x-0.5) * s;
}
else {
err = 1.0e-20;
w = x;
t = 1.0;
if (x < 1.5) {
if (x < err) {
k = (int)x;
y = (double)k - x;
if (Math.abs(y) < err || Math.abs(1.0-y) < err)
ier[0] = -1;
}
if (ier[0] == 0) {
while (w < 1.5) {
t /= w;
w += 1.0;
}
}
}
else {
if (w > 2.5) {
while (w > 2.5) {
w -= 1.0;
t *= w;
}
}
}
w -= 2.0;
g = (((((((0.0021385778 * w - 0.0034961289) * w +
0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
0.999999926;
g *= t;
}
return g;
}
/*****************************************************/
/* Newton法による非線形方程式(f(x)=0)の解 */
/* x1 : 初期値 */
/* eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜f(x(k))｜＜eps2） */
/* max : 最大試行回数 */
/* ind : 実際の試行回数 */
/* （負の時は解を得ることができなかった） */
/* return : 解 */
/*****************************************************/
double newton(double x1, double eps1, double eps2, int max, int ind[])
{
double g, dg, x;
int sw;
x = x1;
ind[0] = 0;
sw = 0;
while (sw == 0 && ind[0] >= 0) {
ind[0]++;
sw = 1;
g = snx(x1);
if (Math.abs(g) > eps2) {
if (ind[0] <= max) {
dg = dsnx(x1);
if (Math.abs(dg) > eps2) {
x = x1 - g / dg;
if (Math.abs(x-x1) > eps1 && Math.abs(x-x1) > eps1*Math.abs(x)) {
x1 = x;
sw = 0;
}
}
else
ind[0] = -1;
}
else
ind[0] = -1;
}
}
return x;
}
/*********************************************************/
/* 二分法による非線形方程式(f(x)=0)の解 */
/* x1,x2 : 初期値 */
/* eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜f(x(k))｜＜eps2） */
/* max : 最大試行回数 */
/* ind : 実際の試行回数 */
/* （負の時は解を得ることができなかった） */
/* return : 解 */
/*********************************************************/
double bisection(double x1, double x2, double eps1, double eps2, int max, int ind[])
{
double f0, f1, f2, x0 = 0.0;
int sw;
f1 = snx(x1);
f2 = snx(x2);
if (f1*f2 > 0.0)
ind[0] = -1;
else {
ind[0] = 0;
if (f1*f2 == 0.0)
x0 = (f1 == 0.0) ? x1 : x2;
else {
sw = 0;
while (sw == 0 && ind[0] >= 0) {
sw = 1;
ind[0] += 1;
x0 = 0.5 * (x1 + x2);
f0 = snx(x0);
if (Math.abs(f0) > eps2) {
if (ind[0] <= max) {
if (Math.abs(x1-x2) > eps1 && Math.abs(x1-x2) > eps1*Math.abs(x2)) {
sw = 0;
if (f0*f1 < 0.0) {
x2 = x0;
f2 = f0;
}
else {
x1 = x0;
f1 = f0;
}
}
}
else
ind[0] = -1;
}
}
}
}
return x0;
}
}
/*
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 // 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 // 因子の名前 水準の数
3.1 4.7 5.1 // 各水準に対する1番目のデータ
4.1 5.6 3.7 // 各水準に対する2番目のデータ
3.3 4.3 4.5 // 各水準に対する3番目のデータ
3.9 5.9 6.0 // 各水準に対する4番目のデータ
3.7 6.1 3.9 // 各水準に対する5番目のデータ
2.4 4.2 5.4 // 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 // 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 // 1番目の因子の名前 その水準の数
品種 2 // 2番目の因子の名前 その水準の数
3 4 12 -4 -4 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ
*/
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>分散分析</TITLE>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
<SCRIPT TYPE="text/javascript">
method = 0;
dof1 = 0;
dof2 = 0;
p = 0.0;
function main()
{
// データ
let Np = new Array();
Np[0] = 1;
Np[1] = 1;
let name = new Array();
let N = parseInt(document.getElementById("N").value);
let a = parseFloat(document.getElementById("a").value);
name[0] = document.getElementById("name1").value;
Np[0] = parseInt(document.getElementById("np_1").value);
let meth = method + 1;
if (meth == 2) {
name[1] = document.getElementById("name2").value;
Np[1] = parseInt(document.getElementById("np_2").value);
}
let x = new Array();
for (let i1 = 0; i1 < Np[0]; i1++) {
x[i1] = new Array();
for (let i2 = 0; i2 < Np[1]; i2++)
x[i1][i2] = new Array();
}
let ss = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
let k = 0;
for (let i1 = 0; i1 < Np[1]; i1++) {
for (let i2 = 0; i2 < N; i2++) {
for (let i3 = 0; i3 < Np[0]; i3++) {
x[i3][i1][i2] = parseFloat(ss[k]);
k++;
}
}
}
// 計算
aov(meth, Np, N, name, a, x);
}
/****************/
/* 関数値の計算 */
/****************/
function snx(sw, x)
{
let y = 0.0;
let w = new Array();
switch (sw) {
// 関数値（f(x)）の計算（正規分布）
case 0:
y = 1.0 - p - normal(x, w);
break;
// 関数値（f(x)）の計算（Ｆ分布）
case 1:
y = f(x, w, dof1, dof2) - 1.0 + p;
break;
// 関数の微分の計算（Ｆ分布）
case 2:
y = f(x, w, dof1, dof2);
y = w[0];
break;
}
return y;
}
/**************************************/
/* 分散分析（一元配置法と二元配置法） */
/* method : 1 or 2 */
/* Np[0] : 因子1の水準数 */
/* [1] : 因子2の水準数 */
/* N : データ数 */
/* name[0] : 因子1の名前 */
/* [1] : 因子2の名前 */
/* a : a ％値 */
/* x : データ */
/* coded by Y.Suganuma */
/**************************************/
function aov(method, Np, N, name, a, x)
{
let xa;
let SP;
let SQ;
let SI;
let SE;
let VP;
let VQ;
let VI;
let VE;
let FP;
let FQ;
let FI;
let i1;
let i2;
let i3;
let sw = new Array();
// 一元配置法
if (method == 1) {
let xi = new Array();
for (i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (i2 = 0; i2 < N; i2++)
xi[i1] += x[i1][0][i2];
xi[i1] /= N;
}
xa = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < N; i2++)
xa += x[i1][0][i2];
}
xa /= (Np[0] * N);
SP = 0.0;
for (i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= N;
SE = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < N; i2++)
SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
}
VP = SP / (Np[0] - 1);
VE = SE / (Np[0] * (N - 1));
FP = VP / VE;
p = 0.01 * a;
dof2 = Np[0] * (N - 1);
str = "";
str = str + "全変動： 平方和 " + (SP+SE) + " 自由度 " + (Np[0]*N-1) + "\n";
dof1 = Np[0] - 1;
str = str + name[0] + "の水準間： 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "％値 " + p_f(sw) + "\n";
str = str + "水準内： 平方和 " + SE + " 自由度 " + (Np[0]*(N-1)) + " 不偏分散 " + VE + "\n";
document.getElementById("ans").value = str;
}
// 二元配置法
else {
let xi = new Array();
for (i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
xi[i1] += x[i1][i2][i3];
}
xi[i1] /= (Np[1] * N);
}
let xj = new Array();
for (i1 = 0; i1 < Np[1]; i1++) {
xj[i1] = 0.0;
for (i2 = 0; i2 < Np[0]; i2++) {
for (i3 = 0; i3 < N; i3++)
xj[i1] += x[i2][i1][i3];
}
xj[i1] /= (Np[0] * N);
}
let xij = new Array();
for (i1 = 0; i1 < Np[0]; i1++) {
xij[i1] = new Array();
for (i2 = 0; i2 < Np[1]; i2++) {
xij[i1][i2] = 0.0;
for (i3 = 0; i3 < N; i3++)
xij[i1][i2] += x[i1][i2][i3];
xij[i1][i2] /= N;
}
}
xa = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
xa += x[i1][i2][i3];
}
}
xa /= (Np[0] * Np[1] * N);
SP = 0.0;
for (i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= (Np[1] * N);
SQ = 0.0;
for (i1 = 0; i1 < Np[1]; i1++)
SQ += (xj[i1] - xa) * (xj[i1] - xa);
SQ *= (Np[0] * N);
SI = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++)
SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
}
SI *= N;
SE = 0.0;
for (i1 = 0; i1 < Np[0]; i1++) {
for (i2 = 0; i2 < Np[1]; i2++) {
for (i3 = 0; i3 < N; i3++)
SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
}
}
VP = SP / (Np[0] - 1);
VQ = SQ / (Np[1] - 1);
VI = SI / ((Np[0] - 1) * (Np[1] - 1));
VE = SE / (Np[0] * Np[1] * (N - 1));
FP = VP / VE;
FQ = VQ / VE;
FI = VI / VE;
p = 0.01 * a;
dof2 = Np[0] * Np[1] * (N - 1);
str = "";
str = str + "全変動： 平方和 " + (SP+SQ+SI+SE) + " 自由度 " + (Np[0]*Np[1]*N-1) + "\n";
dof1 = Np[0] - 1;
str = str + name[0] + "の水準間： 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "％値 " + p_f(sw) + "\n";
dof1 = Np[1] - 1;
str = str + name[1] + "の水準間： 平方和 " + SQ + " 自由度 " + (Np[1]-1) + " 不偏分散 " + VQ + " F " + FQ + " " + a + "％値 " + p_f(sw) + "\n";
dof1 = (Np[0] - 1) * (Np[1] - 1);
str = str + "相互作用： 平方和 " + SI + " 自由度 " + ((Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + VI + " F " + FI + " " + a + "％値 " + p_f(sw) + "\n";
str = str + "水準内： 平方和 " + SE + " 自由度 " + (Np[0]*Np[1]*(N-1)) + " 不偏分散 " + VE + "\n";
document.getElementById("ans").value = str;
}
}
/****************************************/
/* ｆ分布の計算（P(X = ff), P(X < ff)） */
/* dd : P(X = ff) */
/* df1,df2 : 自由度 */
/* return : P(X < ff) */
/****************************************/
function f(ff, dd, df1, df2)
{
let pi = Math.PI;
let pp;
let u;
let x;
let ia;
let ib;
let i1;
if (ff < 1.0e-10)
ff = 1.0e-10;
x = ff * df1 / (ff * df1 + df2);
if(df1%2 == 0) {
if (df2%2 == 0) {
u = x * (1.0 - x);
pp = x;
ia = 2;
ib = 2;
}
else {
u = 0.5 * x * Math.sqrt(1.0-x);
pp = 1.0 - Math.sqrt(1.0-x);
ia = 2;
ib = 1;
}
}
else {
if (df2%2 == 0) {
u = 0.5 * Math.sqrt(x) * (1.0 - x);
pp = Math.sqrt(x);
ia = 1;
ib = 2;
}
else {
u = Math.sqrt(x*(1.0-x)) / pi;
pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
ia = 1;
ib = 1;
}
}
if (ia != df1) {
for (i1 = ia; i1 <= df1-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= x * (i1 + ib) / i1;
}
}
if (ib != df2) {
for (i1 = ib; i1 <= df2-2; i1 += 2) {
pp += 2.0 * u / i1;
u *= (1.0 - x) * (i1 + df1) / i1;
}
}
dd[0] = u / ff;
return pp;
}
/****************************************/
/* ｆ分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/****************************************/
function p_f(ind)
{
let a = 0.0;
let a1 = 0.0;
let b = 0.0;
let b1 = 0.0;
let df1 = 0.0;
let df2 = 0.0;
let e = 0.0;
let ff = 0.0;
let f0;
let x;
let y1;
let y2;
let yq = 0.0;
let sw = 0;
let MAX = 340;
/*
初期値計算の準備
while文は，大きな自由度によるガンマ関数の
オーバーフローを避けるため
*/
while (sw >= 0) {
df1 = 0.5 * (dof1 - sw);
df2 = 0.5 * dof2;
a = 2.0 / (9.0 * (dof1 - sw));
a1 = 1.0 - a;
b = 2.0 / (9.0 * dof2);
b1 = 1.0 - b;
yq = p_normal(ind);
e = b1 * b1 - b * yq * yq;
if (e > 0.8 || (dof1+dof2-sw) <= MAX)
sw = -1;
else {
sw += 1;
if ((dof1-sw) == 0)
sw = -2;
}
}
if (sw == -2)
ind[0] = -1;
/*
f0 : 初期値
*/
else {
if (e > 0.8) {
x = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e;
f0 = Math.pow(x, 3.0);
}
else {
y1 = Math.pow(dof2, df2-1.0);
y2 = Math.pow(dof1, df2);
x = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) * 2.0 * y1 / y2 / p;
f0 = Math.pow(x, 2.0/dof2);
}
/*
ニュートン法
*/
ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind);
}
return ff;
}
/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/* 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-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 < 2; i1++) {
if (sel.options[i1].selected) {
method = i1;
break;
}
}
// 一元配置法
if (method == 0) {
document.getElementById("name2_t").style.display = "none";
document.getElementById("N").value = "6";
document.getElementById("np_1").value = "3";
document.getElementById("data").value = "3.1 4.7 5.1\n4.1 5.6 3.7\n3.3 4.3 4.5\n3.9 5.9 6.0\n3.7 6.1 3.9\n2.4 4.2 5.4";
}
// 二元配置法
else {
document.getElementById("name2_t").style.display = "";
document.getElementById("N").value = "3";
document.getElementById("np_1").value = "5";
document.getElementById("data").value = "3 4 12 -4 -4\n8 -8 31 12 19\n7 -5 8 0 23\n8 -10 9 10 15\n-5 11 26 -1 13\n10 -6 13 -7 -6";
}
}
</SCRIPT>
</HEAD>
<BODY STYLE="font-size: 130%; background-color: #eeffee;">
<H2 STYLE="text-align:center"><B>分散分析</B></H2>
<DL>
<DT> テキストフィールドおよびテキストエリアには，例として，一元配置法及び一元配置法によって分散分析を行う場合に対する値が設定されています．他の問題を実行する場合は，それらを適切に修正してください．
</DL>
<DIV STYLE="text-align:center">
方法：<SELECT ID="method" onChange="set_m()" STYLE="font-size:100%">
<OPTION SELECTED>一元配置法
<OPTION>二元配置法
</SELECT>
データ数：<INPUT ID="N" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="6">
有意水準（％）：<INPUT ID="a" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="5"><BR><BR>
<SPAN ID="name1_t">因子１と水準数：<INPUT ID="name1" STYLE="font-size: 100%;" TYPE="text" SIZE="15" VALUE="因子１"> <INPUT ID="np_1" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="3"> </SPAN>
<SPAN ID="name2_t" STYLE="display: none">因子２と水準数：<INPUT ID="name2" STYLE="font-size: 100%;" TYPE="text" SIZE="15" VALUE="因子２"> <INPUT ID="np_2" STYLE="font-size: 100%;" TYPE="text" SIZE="3" VALUE="2"></SPAN><BR><BR>
データ：<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">
3.1 4.7 5.1
4.1 5.6 3.7
3.3 4.3 4.5
3.9 5.9 6.0
3.7 6.1 3.9
2.4 4.2 5.4
</TEXTAREA>
<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
<TEXTAREA ID="ans" COLS="80" ROWS="10" STYLE="font-size: 100%;"></TEXTAREA>
</DIV>
</BODY>
</HTML>
<?php
/**************************************/
/* 分散分析（一元配置法と二元配置法） */
/* coded by Y.Suganuma */
/**************************************/
$Np = array(1, 1);
$name = array(2);
$p = 0.0;
$dof1 = 1;
$dof2 = 1;
fscanf(STDIN, "%d %d %d", $method, $N, $a);
if ($method == 1 || $method == 2) {
for ($i1 = 0; $i1 < $method; $i1++)
fscanf(STDIN, "%s %d", $name[$i1], $Np[$i1]);
$x = array($Np[0]);
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
$x[$i1] = array($Np[1]);
for ($i2 = 0; $i2 < $Np[1]; $i2++)
$x[$i1][$i2] = array($N);
}
for ($i1 = 0; $i1 < $Np[1]; $i1++) {
for ($i2 = 0; $i2 < $N; $i2++) {
$str = trim(fgets(STDIN));
$x[0][$i1][$i2] = floatval(strtok($str, " "));
for ($i3 = 1; $i3 < $Np[0]; $i3++)
$x[$i3][$i1][$i2] = floatval(strtok(" "));
}
}
aov($method, $Np, $N, $name, $a, $x);
}
else
printf("一元配置法，または，二元配置法だけです！\n");
/**************************************/
/* 分散分析（一元配置法と二元配置法） */
/* method : 1 or 2 */
/* Np[0] : 因子1の水準数 */
/* [1] : 因子2の水準数 */
/* N : データ数 */
/* name[0] : 因子1の名前 */
/* [1] : 因子2の名前 */
/* a : a ％値 */
/* x : データ */
/* coded by Y.Suganuma */
/**************************************/
function aov($method, $Np, $N, $name, $a, $x)
{
global $p, $dof1, $dof2;
// 一元配置法
if ($method == 1) {
$xi = array($Np[0]);
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
$xi[$i1] = 0.0;
for ($i2 = 0; $i2 < $N; $i2++)
$xi[$i1] += $x[$i1][0][$i2];
$xi[$i1] /= $N;
}
$xa = 0.0;
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
for ($i2 = 0; $i2 < $N; $i2++)
$xa += $x[$i1][0][$i2];
}
$xa /= ($Np[0] * $N);
$SP = 0.0;
for ($i1 = 0; $i1 < $Np[0]; $i1++)
$SP += ($xi[$i1] - $xa) * ($xi[$i1] - $xa);
$SP *= $N;
$SE = 0.0;
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
for ($i2 = 0; $i2 < $N; $i2++)
$SE += ($x[$i1][0][$i2] - $xi[$i1]) * ($x[$i1][0][$i2] - $xi[$i1]);
}
$VP = $SP / ($Np[0] - 1);
$VE = $SE / ($Np[0] * ($N - 1));
$FP = $VP / $VE;
$p = 0.01 * $a;
$dof2 = $Np[0] * ($N - 1);
printf("全変動： 平方和 %.2f 自由度 %d\n", $SP+$SE, $Np[0]*$N-1);
$dof1 = $Np[0] - 1;
printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", $name[0], $SP, $Np[0]-1, $VP, $FP, $a, p_f($sw));
printf("水準内： 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", $SE, $Np[0]*($N-1), $VE);
}
// 二元配置法
else {
$xi = array($Np[0]);
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
$xi[$i1] = 0.0;
for ($i2 = 0; $i2 < $Np[1]; $i2++) {
for ($i3 = 0; $i3 < $N; $i3++)
$xi[$i1] += $x[$i1][$i2][$i3];
}
$xi[$i1] /= ($Np[1] * $N);
}
$xj = array($Np[1]);
for ($i1 = 0; $i1 < $Np[1]; $i1++) {
$xj[$i1] = 0.0;
for ($i2 = 0; $i2 < $Np[0]; $i2++) {
for ($i3 = 0; $i3 < $N; $i3++)
$xj[$i1] += $x[$i2][$i1][$i3];
}
$xj[$i1] /= ($Np[0] * $N);
}
$xij = array($Np[0]);
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
$xij[$i1] = array($Np[1]);
for ($i2 = 0; $i2 < $Np[1]; $i2++) {
$xij[$i1][$i2] = 0.0;
for ($i3 = 0; $i3 < $N; $i3++)
$xij[$i1][$i2] += $x[$i1][$i2][$i3];
$xij[$i1][$i2] /= $N;
}
}
$xa = 0.0;
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
for ($i2 = 0; $i2 < $Np[1]; $i2++) {
for ($i3 = 0; $i3 < $N; $i3++)
$xa += $x[$i1][$i2][$i3];
}
}
$xa /= ($Np[0] * $Np[1] * $N);
$SP = 0.0;
for ($i1 = 0; $i1 < $Np[0]; $i1++)
$SP += ($xi[$i1] - $xa) * ($xi[$i1] - $xa);
$SP *= ($Np[1] * $N);
$SQ = 0.0;
for ($i1 = 0; $i1 < $Np[1]; $i1++)
$SQ += ($xj[$i1] - $xa) * ($xj[$i1] - $xa);
$SQ *= ($Np[0] * $N);
$SI = 0.0;
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
for ($i2 = 0; $i2 < $Np[1]; $i2++)
$SI += ($xij[$i1][$i2] - $xi[$i1] - $xj[$i2] + $xa) * ($xij[$i1][$i2] - $xi[$i1] - $xj[$i2] + $xa);
}
$SI *= $N;
$SE = 0.0;
for ($i1 = 0; $i1 < $Np[0]; $i1++) {
for ($i2 = 0; $i2 < $Np[1]; $i2++) {
for ($i3 = 0; $i3 < $N; $i3++)
$SE += ($x[$i1][$i2][$i3] - $xij[$i1][$i2]) * ($x[$i1][$i2][$i3] - $xij[$i1][$i2]);
}
}
$VP = $SP / ($Np[0] - 1);
$VQ = $SQ / ($Np[1] - 1);
$VI = $SI / (($Np[0] - 1) * ($Np[1] - 1));
$VE = $SE / ($Np[0] * $Np[1] * ($N - 1));
$FP = $VP / $VE;
$FQ = $VQ / $VE;
$FI = $VI / $VE;
$p = 0.01 * $a;
$dof2 = $Np[0] * $Np[1] * ($N - 1);
printf("全変動： 平方和 %.2f 自由度 %d\n", $SP+$SQ+$SI+$SE, $Np[0]*$Np[1]*$N-1);
$dof1 = $Np[0] - 1;
printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", $name[0], $SP, $Np[0]-1, $VP, $FP, $a, p_f($sw));
$dof1 = $Np[1] - 1;
printf("%sの水準間： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", $name[1], $SQ, $Np[1]-1, $VQ, $FQ, $a, p_f($sw));
$dof1 = ($Np[0] - 1) * ($Np[1] - 1);
printf("相互作用： 平方和 %.2f 自由度 %d 不偏分散 %.4f F %.2f %d％値 %.2f\n", $SI, ($Np[0]-1)*($Np[1]-1), $VI, $FI, $a, p_f($sw));
printf("水準内： 平方和 %.2f 自由度 %d 不偏分散 %.4f\n", $SE, $Np[0]*$Np[1]*($N-1), $VE);
}
}
/*********************************************************/
/* 二分法による非線形方程式(f(x)=0)の解 */
/* f : f(x)を計算する関数名 */
/* x1,x2 : 初期値 */
/* eps1 : 終了条件１（｜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 = ff), P(X < ff)） */
/* dd : P(X = ff) */
/* df1,df2 : 自由度 */
/* return : P(X < ff) */
/****************************************/
function f($ff, &$dd, $df1, $df2)
{
$pi = 4.0 * atan(1.0);
if ($ff < 1.0e-10)
$ff = 1.0e-10;
$x = $ff * $df1 / ($ff * $df1 + $df2);
if($df1%2 == 0) {
if ($df2%2 == 0) {
$u = $x * (1.0 - $x);
$pp = $x;
$ia = 2;
$ib = 2;
}
else {
$u = 0.5 * $x * sqrt(1.0-$x);
$pp = 1.0 - sqrt(1.0-$x);
$ia = 2;
$ib = 1;
}
}
else {
if ($df2%2 == 0) {
$u = 0.5 * sqrt($x) * (1.0 - $x);
$pp = sqrt($x);
$ia = 1;
$ib = 2;
}
else {
$u = sqrt($x*(1.0-$x)) / $pi;
$pp = 1.0 - 2.0 * atan2(sqrt(1.0-$x), sqrt($x)) / $pi;
$ia = 1;
$ib = 1;
}
}
if ($ia != $df1) {
for ($i1 = $ia; $i1 <= $df1-2; $i1 += 2) {
$pp -= 2.0 * $u / $i1;
$u *= $x * ($i1 + $ib) / $i1;
}
}
if ($ib != $df2) {
for ($i1 = $ib; $i1 <= $df2-2; $i1 += 2) {
$pp += 2.0 * $u / $i1;
$u *= (1.0 - $x) * ($i1 + $df1) / $i1;
}
}
$dd = $u / $ff;
return $pp;
}
/****************************************/
/* ｆ分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/****************************************/
function p_f(&$ind)
{
global $p, $dof1, $dof2;
$ff = 0.0;
$sw = 0;
$MAX = 340;
/*
初期値計算の準備
while文は，大きな自由度によるガンマ関数の
オーバーフローを避けるため
*/
while ($sw >= 0) {
$df1 = 0.5 * ($dof1 - $sw);
$df2 = 0.5 * $dof2;
$a = 2.0 / (9.0 * ($dof1 - $sw));
$a1 = 1.0 - $a;
$b = 2.0 / (9.0 * $dof2);
$b1 = 1.0 - $b;
$yq = p_normal($ind);
$e = $b1 * $b1 - $b * $yq * $yq;
if ($e > 0.8 || ($dof1+$dof2-$sw) <= $MAX)
$sw = -1;
else {
$sw += 1;
if (($dof1-$sw) == 0)
$sw = -2;
}
}
if ($sw == -2)
$ind = -1;
else {
/*
f0 : 初期値
*/
if ($e > 0.8) {
$x = ($a1 * $b1 + $yq * sqrt($a1*$a1*$b+$a*$e)) / $e;
$f0 = pow($x, 3.0);
}
else {
$y1 = pow(floatval($dof2), $df2-1.0);
$y2 = pow(floatval($dof1), $df2);
$x = gamma($df1+$df2, $ind) / gamma($df1, $ind) / gamma($df2, $ind) *
2.0 * $y1 / $y2 / $p;
$f0 = pow($x, 2.0/$dof2);
}
/*
ニュートン法
*/
$ff = newton("f_f", "f_df", $f0, 1.0e-6, 1.0e-10, 100, $ind);
}
return $ff;
}
/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
function f_f($x)
{
global $p, $dof1, $dof2;
return f($x, $y, $dof1, $dof2) - 1.0 + $p;
}
/**************************/
/* P(X = x)（関数の微分） */
/**************************/
function f_df($x)
{
global $dof1, $dof2;
$z = f($x, $y, $dof1, $dof2);
return $y;
}
/*
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 // 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 // 因子の名前 水準の数
3.1 4.7 5.1 // 各水準に対する1番目のデータ
4.1 5.6 3.7 // 各水準に対する2番目のデータ
3.3 4.3 4.5 // 各水準に対する3番目のデータ
3.9 5.9 6.0 // 各水準に対する4番目のデータ
3.7 6.1 3.9 // 各水準に対する5番目のデータ
2.4 4.2 5.4 // 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 // 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 // 1番目の因子の名前 その水準の数
品種 2 // 2番目の因子の名前 その水準の数
3 4 12 -4 -4 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ
*/
?>
############################
# 分散分析
# 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 = ff), P(X < ff)）
# dd : P(X = ff)
# df1,df2 : 自由度
# return : P(X < ff)
#######################################
def f(ff, dd, df1, df2)
if ff < 1.0e-10
ff = 1.0e-10
end
x = ff * df1 / (ff * df1 + df2)
if df1%2 == 0
if df2%2 == 0
u = x * (1.0 - x)
pp = x
ia = 2
ib = 2
else
u = 0.5 * x * Math.sqrt(1.0-x)
pp = 1.0 - Math.sqrt(1.0-x)
ia = 2
ib = 1
end
else
if df2%2 == 0
u = 0.5 * Math.sqrt(x) * (1.0 - x)
pp = Math.sqrt(x)
ia = 1
ib = 2
else
u = Math.sqrt(x*(1.0-x)) / Math::PI
pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / Math::PI
ia = 1
ib = 1
end
end
if ia != df1
i1 = ia
while i1 < df1-1
pp -= 2.0 * u / i1
u *= x * (i1 + ib) / i1
i1 += 2
end
end
if ib != df2
i1 = ib
while i1 < df2-1
pp += 2.0 * u / i1
u *= (1.0 - x) * (i1 + df1) / i1
i1 += 2
end
end
dd[0] = u / ff
return pp
end
############################################
# 二分法による非線形方程式(f(x)=0)の解
# x1,x2 : 初期値
# eps1 : 終了条件１（｜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_f(ind)
f_snx = Proc.new { |sw, x|
y = Array.new(1)
z = f(x, y, $dof1, $dof2)
if sw == 0
z - 1.0 + $p
else
z = y[0]
end
}
max = 340
ff = 0.0
sw = 0
# 初期値計算の準備
# while文は，大きな自由度によるガンマ関数の
# オーバーフローを避けるため
while sw >= 0
df1 = 0.5 * ($dof1 - sw)
df2 = 0.5 * $dof2
a = 2.0 / (9.0 * ($dof1 - sw))
a1 = 1.0 - a
b = 2.0 / (9.0 * $dof2)
b1 = 1.0 - b
yq = p_normal(ind)
e = b1 * b1 - b * yq * yq
if e > 0.8 or $dof1+$dof2-sw <= max
sw = -1
else
sw += 1
if ($dof1-sw) == 0
sw = -2
end
end
end
if sw == -2
ind[0] = -1
else
# f0 : 初期値
if e > 0.8
x = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e
f0 = x ** 3.0
else
y1 = Float($dof2) ** (df2-1.0)
y2 = Float($dof1) ** df2
x = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) * 2.0 * y1 / y2 / $p
f0 = x ** (2.0/$dof2)
end
# ニュートン法
ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind, &f_snx)
end
return ff
end
#######################################
# 分散分析（一元配置法と二元配置法）
# method : 1 or 2
# np[0] : 因子1の水準数
# [1] : 因子2の水準数
# nn : データ数
# name[0] : 因子1の名前
# [1] : 因子2の名前
# a : a ％値
# x : データ
# coded by Y.Suganuma
#######################################
def aov(method, np, nn, name, a, x)
# 一元配置法
if method == 1
xi = Array.new(np[0])
for i1 in 0 ... np[0]
xi[i1] = 0.0
for i2 in 0 ... nn
xi[i1] += x[i1][0][i2]
end
xi[i1] /= nn
end
xa = 0.0
for i1 in 0 ... np[0]
for i2 in 0 ... nn
xa += x[i1][0][i2]
end
end
xa /= (np[0] * nn)
sp = 0.0
for i1 in 0 ... np[0]
sp += (xi[i1] - xa) * (xi[i1] - xa)
end
sp *= nn
se = 0.0
for i1 in 0 ... np[0]
for i2 in 0 ... nn
se += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1])
end
end
vp = sp / (np[0] - 1)
ve = se / (np[0] * (nn - 1))
fp = vp / ve
$p = 0.01 * a
$dof2 = np[0] * (nn - 1)
sw = Array.new(1)
print("全変動： 平方和 " + String(sp+se) + " 自由度 " + String(np[0]*nn-1) + "\n")
$dof1 = np[0] - 1
print(name[0])
print("の水準間： 平方和 " + String(sp) + " 自由度 " + String(np[0]-1) + " 不偏分散 " + String(vp) + " F " + String(fp) + " " + String(a) + "％値 " + String(p_f(sw)) + "\n")
print("水準内： 平方和 " + String(se) + " 自由度 " + String(np[0]*(nn-1)) + " 不偏分散 " + String(ve) + "\n")
# 二元配置法
else
xi = Array.new(np[0])
for i1 in 0 ... np[0]
xi[i1] = 0.0
for i2 in 0 ... np[1]
for i3 in 0 ... nn
xi[i1] += x[i1][i2][i3]
end
end
xi[i1] /= (np[1] * nn)
end
xj = Array.new(np[1])
for i1 in 0 ... np[1]
xj[i1] = 0.0
for i2 in 0 ... np[0]
for i3 in 0 ... nn
xj[i1] += x[i2][i1][i3]
end
end
xj[i1] /= (np[0] * nn)
end
xij = Array.new(np[0])
for i1 in 0 ... np[0]
xij[i1] = Array.new(np[1])
for i2 in 0 ... np[1]
xij[i1][i2] = 0.0
for i3 in 0 ... nn
xij[i1][i2] += x[i1][i2][i3]
end
xij[i1][i2] /= nn
end
end
xa = 0.0
for i1 in 0 ... np[0]
for i2 in 0 ... np[1]
for i3 in 0 ... nn
xa += x[i1][i2][i3]
end
end
end
xa /= (np[0] * np[1] * nn)
sp = 0.0
for i1 in 0 ... np[0]
sp += (xi[i1] - xa) * (xi[i1] - xa)
end
sp *= (np[1] * nn)
sq = 0.0
for i1 in 0 ... np[1]
sq += (xj[i1] - xa) * (xj[i1] - xa)
end
sq *= (np[0] * nn)
si = 0.0
for i1 in 0 ... np[0]
for i2 in 0 ... np[1]
si += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa)
end
end
si *= nn
se = 0.0
for i1 in 0 ... np[0]
for i2 in 0 ... np[1]
for i3 in 0 ... nn
se += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2])
end
end
end
vp = sp / (np[0] - 1)
vq = sq / (np[1] - 1)
vi = si / ((np[0] - 1) * (np[1] - 1))
ve = se / (np[0] * np[1] * (nn - 1))
fp = vp / ve
fq = vq / ve
fi = vi / ve
$p = 0.01 * a
$dof2 = np[0] * np[1] * (nn - 1)
sw = Array.new(1)
print("全変動： 平方和 " + String(sp+sq+si+se) + " 自由度 " + String(np[0]*np[1]*nn-1) + "\n")
$dof1 = np[0] - 1
print(name[0])
print("の水準間： 平方和 " + String(sp) + " 自由度 " + String(np[0]-1) + " 不偏分散 " + String(vp) + " F " + String(fp) + " " + String(a) + "％値 " + String(p_f(sw)) + "\n")
$dof1 = np[1] - 1
print(name[1])
print("の水準間： 平方和 " + String(sq) + " 自由度 " + String(np[1]-1) + " 不偏分散 " + String(vq) + " F " + String(fq) + " " + String(a) + "％値 " + String(p_f(sw)) + "\n")
$dof1 = (np[0] - 1) * (np[1] - 1)
print("相互作用： 平方和 " + String(si) + " 自由度 " + String((np[0]-1)*(np[1]-1)) + " 不偏分散 " + String(vi) + " F " + String(fi) + " " + String(a) + "％値 " + String(p_f(sw)) + "\n")
print("水準内： 平方和 " + String(se) + " 自由度 " + String(np[0]*np[1]*(nn-1)) + " 不偏分散 " + String(ve) + "\n")
end
end
$p = 0.0
$dof1 = 1
$dof2 = 2
name = ["", ""]
np = [1, 1]
ss = gets().split(" ")
method = Integer(ss[0]) # 因子の数
nn = Integer(ss[1]) # 各水準におけるデータ数
a = Float(ss[2]) # 有意水準（％）
if method == 1 or method == 2
for i1 in 0 ... method
ss = gets().split(" ")
name[i1] = ss[0]
np[i1] = Integer(ss[1])
end
x = Array.new(np[0])
for i1 in 0 ... np[0]
x[i1] = Array.new(np[1])
for i2 in 0 ... np[1]
x[i1][i2] = Array.new(nn)
end
end
for i1 in 0 ... np[1]
for i2 in 0 ... nn
ss = gets().split(" ")
for i3 in 0 ... np[0]
x[i3][i1][i2] = Float(ss[i3])
end
end
end
aov(method, np, nn, name, a, x)
else
print("一元配置法，または，二元配置法だけです！\n")
end
=begin
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 # 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 # 因子の名前 水準の数
3.1 4.7 5.1 # 各水準に対する1番目のデータ
4.1 5.6 3.7 # 各水準に対する2番目のデータ
3.3 4.3 4.5 # 各水準に対する3番目のデータ
3.9 5.9 6.0 # 各水準に対する4番目のデータ
3.7 6.1 3.9 # 各水準に対する5番目のデータ
2.4 4.2 5.4 # 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 # 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 # 1番目の因子の名前 その水準の数
品種 2 # 2番目の因子の名前 その水準の数
3 4 12 -4 -4 # 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 # 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 # 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 # 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 # 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 # 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ
=end
# -*- 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 = ff), P(X < ff)）
# dd : P(X = ff)
# df1,df2 : 自由度
# return : P(X < ff)
#######################################
def f(ff, dd, df1, df2) :
if ff < 1.0e-10 :
ff = 1.0e-10
x = ff * df1 / (ff * df1 + df2)
if df1%2 == 0 :
if df2%2 == 0 :
u = x * (1.0 - x)
pp = x
ia = 2
ib = 2
else :
u = 0.5 * x * sqrt(1.0-x)
pp = 1.0 - sqrt(1.0-x)
ia = 2
ib = 1
else :
if df2%2 == 0 :
u = 0.5 * sqrt(x) * (1.0 - x)
pp = sqrt(x)
ia = 1
ib = 2
else :
u = sqrt(x*(1.0-x)) / pi
pp = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi
ia = 1
ib = 1
if ia != df1 :
for i1 in range(ia, df1-1, 2) :
pp -= 2.0 * u / i1
u *= x * (i1 + ib) / i1
if ib != df2 :
for i1 in range(ib, df2-1, 2) :
pp += 2.0 * u / i1
u *= (1.0 - x) * (i1 + df1) / i1
dd[0] = u / ff
return pp
############################################
# 二分法による非線形方程式(f(x)=0)の解
# fn : f(x)を計算する関数名
# x1,x2 : 初期値
# eps1 : 終了条件１（｜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
##########################################
# 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)（関数値, f 分布）
#######################################
def f_f(x) :
y = np.empty(1, np.float)
return f(x, y, dof1, dof2) - 1.0 + p
#################################
# P(X = x)（関数の微分, f 分布）
#################################
def f_df(x) :
y = np.empty(1, np.float)
z = f(x, y, dof1, dof2)
return y[0]
######################################
# ｆ分布のｐ％値（P(X > u) = 0.01p）
# ind : >= 0 : normal（収束回数）
# = -1 : 収束しなかった
######################################
def p_f(ind) :
MAX = 340
ff = 0.0
sw = 0
# 初期値計算の準備
# while文は，大きな自由度によるガンマ関数の
# オーバーフローを避けるため
while sw >= 0 :
df1 = 0.5 * (dof1 - sw)
df2 = 0.5 * dof2
a = 2.0 / (9.0 * (dof1 - sw))
a1 = 1.0 - a
b = 2.0 / (9.0 * dof2)
b1 = 1.0 - b
yq = p_normal(ind)
e = b1 * b1 - b * yq * yq
if e > 0.8 or dof1+dof2-sw <= MAX :
sw = -1
else :
sw += 1
if (dof1-sw) == 0 :
sw = -2
if sw == -2 :
ind[0] = -1
else :
# f0 : 初期値
if e > 0.8 :
x = (a1 * b1 + yq * sqrt(a1*a1*b+a*e)) / e
f0 = x ** 3.0
else :
y1 = float(dof2) ** (df2-1.0)
y2 = float(dof1) ** df2
x = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) * 2.0 * y1 / y2 / p
f0 = x ** (2.0/dof2)
# ニュートン法
ff = newton(f_f, f_df, f0, 1.0e-6, 1.0e-10, 100, ind)
return ff
#######################################
# 分散分析（一元配置法と二元配置法）
# method : 1 or 2
# Np[0] : 因子1の水準数
# [1] : 因子2の水準数
# N : データ数
# name[0] : 因子1の名前
# [1] : 因子2の名前
# a : a ％値
# x : データ
# coded by Y.Suganuma
#######################################
def aov(method, Np, N, name, a, x) :
global p, dof1, dof2
# 一元配置法
if method == 1 :
xi = np.empty(Np[0], np.float)
for i1 in range(0, Np[0]) :
xi[i1] = 0.0
for i2 in range(0, N) :
xi[i1] += x[i1][0][i2]
xi[i1] /= N
xa = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, N) :
xa += x[i1][0][i2]
xa /= (Np[0] * N)
SP = 0.0
for i1 in range(0, Np[0]) :
SP += (xi[i1] - xa) * (xi[i1] - xa)
SP *= N
SE = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, N) :
SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1])
VP = SP / (Np[0] - 1)
VE = SE / (Np[0] * (N - 1))
FP = VP / VE
p = 0.01 * a
dof2 = Np[0] * (N - 1)
sw = np.empty(1, np.int)
print("全変動： 平方和 " + str(SP+SE) + " 自由度 " + str(Np[0]*N-1))
dof1 = Np[0] - 1
print(name[0] + "の水準間： 平方和 " + str(SP) + " 自由度 " + str(Np[0]-1) + " 不偏分散 " + str(VP) + " F " + str(FP) + " " + str(a) + "％値 " + str(p_f(sw)))
print("水準内： 平方和 " + str(SE) + " 自由度 " + str(Np[0]*(N-1)) + " 不偏分散 " + str(VE))
# 二元配置法
else :
xi = np.empty(Np[0], np.float)
for i1 in range(0, Np[0]) :
xi[i1] = 0.0
for i2 in range(0, Np[1]) :
for i3 in range(0, N) :
xi[i1] += x[i1][i2][i3]
xi[i1] /= (Np[1] * N)
xj = np.empty(Np[1], np.float)
for i1 in range(0, Np[1]) :
xj[i1] = 0.0
for i2 in range(0, Np[0]) :
for i3 in range(0, N) :
xj[i1] += x[i2][i1][i3]
xj[i1] /= (Np[0] * N)
xij = np.empty((Np[0], Np[1]), np.float)
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
xij[i1][i2] = 0.0
for i3 in range(0, N) :
xij[i1][i2] += x[i1][i2][i3]
xij[i1][i2] /= N
xa = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
for i3 in range(0, N) :
xa += x[i1][i2][i3]
xa /= (Np[0] * Np[1] * N)
SP = 0.0
for i1 in range(0, Np[0]) :
SP += (xi[i1] - xa) * (xi[i1] - xa)
SP *= (Np[1] * N)
SQ = 0.0
for i1 in range(0, Np[1]) :
SQ += (xj[i1] - xa) * (xj[i1] - xa)
SQ *= (Np[0] * N)
SI = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa)
SI *= N
SE = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
for i3 in range(0, N) :
SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2])
VP = SP / (Np[0] - 1)
VQ = SQ / (Np[1] - 1)
VI = SI / ((Np[0] - 1) * (Np[1] - 1))
VE = SE / (Np[0] * Np[1] * (N - 1))
FP = VP / VE
FQ = VQ / VE
FI = VI / VE
p = 0.01 * a
dof2 = Np[0] * Np[1] * (N - 1)
sw = np.empty(1, np.int)
print("全変動： 平方和 " + str(SP+SQ+SI+SE) + " 自由度 " + str(Np[0]*Np[1]*N-1))
dof1 = Np[0] - 1
print(name[0] + "の水準間： 平方和 " + str(SP) + " 自由度 " + str(Np[0]-1) + " 不偏分散 " + str(VP) + " F " + str(FP) + " " + str(a) + "％値 " + str(p_f(sw)))
dof1 = Np[1] - 1
print(name[1] + "の水準間： 平方和 " + str(SQ) + " 自由度 " + str(Np[1]-1) + " 不偏分散 " + str(VQ) + " F " + str(FQ) + " " + str(a) + "％値 " + str(p_f(sw)))
dof1 = (Np[0] - 1) * (Np[1] - 1)
print("相互作用： 平方和 " + str(SI) + " 自由度 " + str((Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + str(VI) + " F " + str(FI) + " " + str(a) + "％値 " + str(p_f(sw)))
print("水準内： 平方和 " + str(SE) + " 自由度 " + str(Np[0]*Np[1]*(N-1)) + " 不偏分散 " + str(VE))
p = 0.0
dof1 = 1
dof2 = 2
############################
# 分散分析
# coded by Y.Suganuma
############################
name = ["", ""]
Np = np.array([1, 1], np.int)
line = sys.stdin.readline()
ss = line.split()
method = int(ss[0]) # 因子の数
N = int(ss[1]) # 各水準におけるデータ数
a = float(ss[2]) # 有意水準（％）
if method == 1 or method == 2 :
for i1 in range(0, method) :
line = sys.stdin.readline()
ss = line.split()
name[i1] = ss[0]
Np[i1] = int(ss[1])
x = np.empty((Np[0], Np[1], N), np.float)
for i1 in range(0, Np[1]) :
for i2 in range(0, N) :
line = sys.stdin.readline()
ss = line.split()
for i3 in range(0, Np[0]) :
x[i3][i1][i2] = float(ss[i3])
aov(method, Np, N, name, a, x)
else :
print("一元配置法，または，二元配置法だけです！")
"""
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 # 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 # 因子の名前 水準の数
3.1 4.7 5.1 # 各水準に対する1番目のデータ
4.1 5.6 3.7 # 各水準に対する2番目のデータ
3.3 4.3 4.5 # 各水準に対する3番目のデータ
3.9 5.9 6.0 # 各水準に対する4番目のデータ
3.7 6.1 3.9 # 各水準に対する5番目のデータ
2.4 4.2 5.4 # 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 # 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 # 1番目の因子の名前 その水準の数
品種 2 # 2番目の因子の名前 その水準の数
3 4 12 -4 -4 # 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 # 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 # 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 # 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 # 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 # 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ
"""
/****************************/
/* 分散分析 */
/* coded by Y.Suganuma */
/****************************/
using System;
class Program
{
static void Main()
{
Test1 ts = new Test1();
}
}
class Test1
{
double p; // α％値を計算するとき時α/100を設定
int dof1, dof2; // 自由度
public Test1()
{
int[] Np = {1, 1};
char[] charSep = new char[] {' '};
String[] str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
int method = int.Parse(str[0]);
int N = int.Parse(str[1]);
int a = int.Parse(str[2]);
String[] name = new String [2];
if (method == 1 || method == 2) {
for (int i1 = 0; i1 < method; i1++) {
str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
name[i1] = str[0].Trim();
Np[i1] = int.Parse(str[1]);
}
double[][][]x = new double [Np[0]][][];
for (int i1 = 0; i1 < Np[0]; i1++) {
x[i1] = new double [Np[1]][];
for (int i2 = 0; i2 < Np[1]; i2++)
x[i1][i2] = new double [N];
}
for (int i1 = 0; i1 < Np[1]; i1++) {
for (int i2 = 0; i2 < N; i2++) {
str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
for (int i3 = 0; i3 < Np[0]; i3++) {
x[i3][i1][i2] = double.Parse(str[i3]);
}
}
}
aov(method, Np, N, name, a, x);
}
else
Console.WriteLine("一元配置法，または，二元配置法だけです！");
}
/**************************************/
/* 分散分析（一元配置法と二元配置法） */
/* method : 1 or 2 */
/* Np[0] : 因子1の水準数 */
/* [1] : 因子2の水準数 */
/* N : データ数 */
/* name[0] : 因子1の名前 */
/* [1] : 因子2の名前 */
/* a : a ％値 */
/* x : データ */
/* coded by Y.Suganuma */
/**************************************/
void aov(int method, int[] Np, int N, String[] name, int a, double[][][] x)
{
int sw = 0;
// 一元配置法
if (method == 1) {
double[] xi = new double [Np[0]];
for (int i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (int i2 = 0; i2 < N; i2++)
xi[i1] += x[i1][0][i2];
xi[i1] /= N;
}
double xa = 0.0;
for (int i1 = 0; i1 < Np[0]; i1++) {
for (int i2 = 0; i2 < N; i2++)
xa += x[i1][0][i2];
}
xa /= (Np[0] * N);
double SP = 0.0;
for (int i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= N;
double SE = 0.0;
for (int i1 = 0; i1 < Np[0]; i1++) {
for (int i2 = 0; i2 < N; i2++)
SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1]);
}
double VP = SP / (Np[0] - 1);
double VE = SE / (Np[0] * (N - 1));
double FP = VP / VE;
p = 0.01 * a;
dof2 = Np[0] * (N - 1);
Console.WriteLine("全変動： 平方和 " + (SP+SE) + " 自由度 " + (Np[0]*N-1));
dof1 = Np[0] - 1;
Console.WriteLine(name[0] + "の水準間： 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "％値 " + p_f(ref sw));
Console.WriteLine("水準内： 平方和 " + SE + " 自由度 " + (Np[0]*(N-1))+ " 不偏分散 " + VE);
}
// 二元配置法
else {
double[] xi = new double [Np[0]];
for (int i1 = 0; i1 < Np[0]; i1++) {
xi[i1] = 0.0;
for (int i2 = 0; i2 < Np[1]; i2++) {
for (int i3 = 0; i3 < N; i3++)
xi[i1] += x[i1][i2][i3];
}
xi[i1] /= (Np[1] * N);
}
double[] xj = new double [Np[1]];
for (int i1 = 0; i1 < Np[1]; i1++) {
xj[i1] = 0.0;
for (int i2 = 0; i2 < Np[0]; i2++) {
for (int i3 = 0; i3 < N; i3++)
xj[i1] += x[i2][i1][i3];
}
xj[i1] /= (Np[0] * N);
}
double[][] xij = new double [Np[0]][];
for (int i1 = 0; i1 < Np[0]; i1++) {
xij[i1] = new double [Np[1]];
for (int i2 = 0; i2 < Np[1]; i2++) {
xij[i1][i2] = 0.0;
for (int i3 = 0; i3 < N; i3++)
xij[i1][i2] += x[i1][i2][i3];
xij[i1][i2] /= N;
}
}
double xa = 0.0;
for (int i1 = 0; i1 < Np[0]; i1++) {
for (int i2 = 0; i2 < Np[1]; i2++) {
for (int i3 = 0; i3 < N; i3++)
xa += x[i1][i2][i3];
}
}
xa /= (Np[0] * Np[1] * N);
double SP = 0.0;
for (int i1 = 0; i1 < Np[0]; i1++)
SP += (xi[i1] - xa) * (xi[i1] - xa);
SP *= (Np[1] * N);
double SQ = 0.0;
for (int i1 = 0; i1 < Np[1]; i1++)
SQ += (xj[i1] - xa) * (xj[i1] - xa);
SQ *= (Np[0] * N);
double SI = 0.0;
for (int i1 = 0; i1 < Np[0]; i1++) {
for (int i2 = 0; i2 < Np[1]; i2++)
SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa);
}
SI *= N;
double SE = 0.0;
for (int i1 = 0; i1 < Np[0]; i1++) {
for (int i2 = 0; i2 < Np[1]; i2++) {
for (int i3 = 0; i3 < N; i3++)
SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2]);
}
}
double VP = SP / (Np[0] - 1);
double VQ = SQ / (Np[1] - 1);
double VI = SI / ((Np[0] - 1) * (Np[1] - 1));
double VE = SE / (Np[0] * Np[1] * (N - 1));
double FP = VP / VE;
double FQ = VQ / VE;
double FI = VI / VE;
p = 0.01 * a;
dof2 = Np[0] * Np[1] * (N - 1);
Console.WriteLine("全変動： 平方和 " + (SP+SQ+SI+SE) + " 自由度 " + (Np[0]*Np[1]*N-1));
dof1 = Np[0] - 1;
Console.WriteLine(name[0] + "の水準間： 平方和 " + SP + " 自由度 " + (Np[0]-1) + " 不偏分散 " + VP + " F " + FP + " " + a + "％値 " + p_f(ref sw));
dof1 = Np[1] - 1;
Console.WriteLine(name[1] + "の水準間： 平方和 " + SQ + " 自由度 " + (Np[1]-1) + " 不偏分散 " + VQ + " F " + FQ + " " + a + "％値 " + p_f(ref sw));
dof1 = (Np[0] - 1) * (Np[1] - 1);
Console.WriteLine("相互作用： 平方和 " + SI + " 自由度 " + ((Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + VI + " F " + FI + " " + a + "％値 " + p_f(ref sw));
Console.WriteLine("水準内： 平方和 " + SE + " 自由度 " + (Np[0]*Np[1]*(N-1)) + " 不偏分散 " + VE);
}
}
/****************************************/
/* ｆ分布の計算（P(X = ff), P(X < ff)） */
/* dd : P(X = ff) */
/* df1,df2 : 自由度 */
/* return : P(X < ff) */
/****************************************/
double f(double ff, ref double dd, int df1, int df2)
{
double pp, u;
int ia, ib;
if (ff < 1.0e-10)
ff = 1.0e-10;
double x = ff * df1 / (ff * df1 + df2);
if (df1%2 == 0) {
if (df2%2 == 0) {
u = x * (1.0 - x);
pp = x;
ia = 2;
ib = 2;
}
else {
u = 0.5 * x * Math.Sqrt(1.0-x);
pp = 1.0 - Math.Sqrt(1.0-x);
ia = 2;
ib = 1;
}
}
else {
if (df2%2 == 0) {
u = 0.5 * Math.Sqrt(x) * (1.0 - x);
pp = Math.Sqrt(x);
ia = 1;
ib = 2;
}
else {
u = Math.Sqrt(x*(1.0-x)) / Math.PI;
pp = 1.0 - 2.0 * Math.Atan2(Math.Sqrt(1.0-x), Math.Sqrt(x)) / Math.PI;
ia = 1;
ib = 1;
}
}
if (ia != df1) {
for (int i1 = ia; i1 <= df1-2; i1 += 2) {
pp -= 2.0 * u / i1;
u *= x * (i1 + ib) / i1;
}
}
if (ib != df2) {
for (int i1 = ib; i1 <= df2-2; i1 += 2) {
pp += 2.0 * u / i1;
u *= (1.0 - x) * (i1 + df1) / i1;
}
}
dd = u / ff;
return pp;
}
/****************************************/
/* ｆ分布のｐ％値（P(X > u) = 0.01p） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/****************************************/
double p_f(ref int ind)
{
int MAX = 340;
double a = 0.0, a1 = 0.0, b = 0.0, b1 = 0.0, df1 = 0.0, df2 = 0.0, e = 0.0, ff = 0.0, yq = 0.0;
int sw = 0;
/*
初期値計算の準備
while文は，大きな自由度によるガンマ関数の
オーバーフローを避けるため
*/
while (sw >= 0) {
df1 = 0.5 * (dof1 - sw);
df2 = 0.5 * dof2;
a = 2.0 / (9.0 * (dof1 - sw));
a1 = 1.0 - a;
b = 2.0 / (9.0 * dof2);
b1 = 1.0 - b;
yq = p_normal(ref ind);
e = b1 * b1 - b * yq * yq;
if (e > 0.8 || (dof1+dof2-sw) <= MAX)
sw = -1;
else {
sw += 1;
if ((dof1-sw) == 0)
sw = -2;
}
}
if (sw == -2)
ind = -1;
else {
/*
f0 : 初期値
*/
double f0 = 0.0;
if (e > 0.8) {
double x = (a1 * b1 + yq * Math.Sqrt(a1*a1*b+a*e)) / e;
f0 = Math.Pow(x, 3.0);
}
else {
double y1 = Math.Pow((double)dof2, df2-1.0);
double y2 = Math.Pow((double)dof1, df2);
double x = gamma(df1+df2, ref ind) / gamma(df1, ref ind) / gamma(df2, ref ind) *
2.0 * y1 / y2 / p;
f0 = Math.Pow(x, 2.0/dof2);
}
/*
ニュートン法
*/
ff = newton(f0, 1.0e-6, 1.0e-10, 100, ref ind, f_f, f_df);
}
return ff;
}
/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/* ier : =0 : normal */
/* =-1 : x=-n (n=0,1,2,･･･) */
/* return : 結果 */
/****************************************/
static double gamma(double x, ref int ier)
{
double g;
ier = 0;
if (x > 5.0) {
double v = 1.0 / x;
double s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
0.00078403922172) * v - 0.000229472093621) * v -
0.00268132716049) * v + 0.00347222222222) * v +
0.0833333333333) * v + 1.0;
g = 2.506628274631001 * Math.Exp(-x) * Math.Pow(x,x-0.5) * s;
}
else {
double err = 1.0e-20;
double w = x;
double t = 1.0;
if (x < 1.5) {
if (x < err) {
int k = (int)x;
double y = (double)k - x;
if (Math.Abs(y) < err || Math.Abs(1.0-y) < err)
ier = -1;
}
if (ier == 0) {
while (w < 1.5) {
t /= w;
w += 1.0;
}
}
}
else {
if (w > 2.5) {
while (w > 2.5) {
w -= 1.0;
t *= w;
}
}
}
w -= 2.0;
g = (((((((0.0021385778 * w - 0.0034961289) * w +
0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
0.999999926;
g *= t;
}
return g;
}
/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
double f_f(double x)
{
double y = 0.0;
return f(x, ref y, dof1, dof2) - 1.0 + p;
}
/**************************/
/* P(X = x)（関数の微分） */
/**************************/
double f_df(double x)
{
double y = 0.0;
double z = f(x, ref y, dof1, dof2);
return y;
}
/*****************************************************/
/* Newton法による非線形方程式(f(x)=0)の解 */
/* x1 : 初期値 */
/* eps1 : 終了条件１（｜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;
}
}
/*
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 // 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 // 因子の名前 水準の数
3.1 4.7 5.1 // 各水準に対する1番目のデータ
4.1 5.6 3.7 // 各水準に対する2番目のデータ
3.3 4.3 4.5 // 各水準に対する3番目のデータ
3.9 5.9 6.0 // 各水準に対する4番目のデータ
3.7 6.1 3.9 // 各水準に対する5番目のデータ
2.4 4.2 5.4 // 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 // 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 // 1番目の因子の名前 その水準の数
品種 2 // 2番目の因子の名前 その水準の数
3 4 12 -4 -4 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ
*/
'**************************'
' 分散分析 '
' coded by Y.Suganuma '
'**************************'
Imports System.Text.RegularExpressions
Module Test
Dim p As Double ' α％値を計算するとき時α/100を設定
Dim dof1 As Integer ' 自由度
Dim dof2 As Integer ' 自由度
Sub Main()
Dim Np() As Integer = {1, 1}
Dim MS As Regex = New Regex("\s+")
Dim str() As String = MS.Split(Console.ReadLine().Trim())
Dim method As Integer = Integer.Parse(str(0))
Dim N As Integer = Integer.Parse(str(1))
Dim a As Integer = Integer.Parse(str(2))
Dim name(2) As String
If method = 1 or method = 2
For i1 As Integer = 0 To method-1
str = MS.Split(Console.ReadLine().Trim())
name(i1) = str(0)
Np(i1) = Integer.Parse(str(1))
Next
Dim x(Np(0),Np(1),N) As Double
For i1 As Integer = 0 To Np(1)-1
For i2 As Integer = 0 To N-1
str = MS.Split(Console.ReadLine().Trim())
For i3 As Integer = 0 To Np(0)-1
x(i3,i1,i2) = Double.Parse(str(i3))
Next
Next
Next
aov(method, Np, N, name, a, x)
Else
Console.WriteLine("一元配置法，または，二元配置法だけです！")
End If
End Sub
'************************************'
' 分散分析（一元配置法と二元配置法） '
' method : 1 or 2 '
' Np(0) : 因子1の水準数 '
' (1) : 因子2の水準数 '
' N : データ数 '
' name(0) : 因子1の名前 '
' (1) : 因子2の名前 '
' a : a ％値 '
' x : データ '
' coded by Y.Suganuma '
'************************************'
Sub aov(method As Integer, Np() As Integer, N As Integer, name() As String, a As Integer, x(,,) As Double)
' P(X > x) - 1 + p（関数値）（ラムダ式）
Dim f_f = Function(v1) As Double
Dim v2 As Double = 0.0
Return f(v1, v2, dof1, dof2) - 1.0 + p
End Function
' P(X = x)（関数の微分）（ラムダ式）
Dim f_df = Function(v1) As Double
Dim v2 As Double = 0.0
Dim v3 As Double = f(v1, v2, dof1, dof2)
Return v2
End Function
' 1.0 - p - P(X>x)（関数値）（ラムダ式）
Dim normal_f = Function(v1) As Double
Dim v2 As Double = 0.0
Return 1.0 - p - normal(v1, v2)
End Function
Dim sw As Integer = 0
' 一元配置法
If method = 1
Dim xi(Np(0)) As Double
For i1 As Integer = 0 To Np(0)-1
xi(i1) = 0.0
For i2 As Integer = 0 To N-1
xi(i1) += x(i1,0,i2)
Next
xi(i1) /= N
Next
Dim xa As Double = 0.0
For i1 As Integer = 0 To Np(0)-1
For i2 As Integer = 0 To N-1
xa += x(i1,0,i2)
Next
Next
xa /= (Np(0) * N)
Dim SP As Double = 0.0
For i1 As Integer = 0 To Np(0)-1
SP += (xi(i1) - xa) * (xi(i1) - xa)
Next
SP *= N
Dim SE As Double = 0.0
For i1 As Integer = 0 To Np(0)-1
For i2 As Integer = 0 To N-1
SE += (x(i1,0,i2) - xi(i1)) * (x(i1,0,i2) - xi(i1))
Next
Next
Dim VP As Double = SP / (Np(0) - 1)
Dim VE As Double = SE / (Np(0) * (N - 1))
Dim FP As Double = VP / VE
p = 0.01 * a
dof2 = Np(0) * (N - 1)
Console.WriteLine("全変動： 平方和 " & (SP+SE) & " 自由度 " & (Np(0)*N-1))
dof1 = Np(0) - 1
Console.WriteLine(name(0) & "の水準間： 平方和 " & SP & " 自由度 " & (Np(0)-1) & " 不偏分散 " & VP & " F " & FP & " " & a & "％値 " & p_f(sw, f_f, f_df, normal_f))
Console.WriteLine("水準内： 平方和 " & SE & " 自由度 " & (Np(0)*(N-1)) & " 不偏分散 " & VE)
' 二元配置法
Else
Dim xi(Np(0)) As Double
For i1 As Integer = 0 To Np(0)-1
xi(i1) = 0.0
For i2 As Integer = 0 To Np(1)-1
For i3 As Integer = 0 To N-1
xi(i1) += x(i1,i2,i3)
Next
Next
xi(i1) /= (Np(1) * N)
Next
Dim xj(Np(1)) As Double
For i1 As Integer = 0 To Np(1)-1
xj(i1) = 0.0
For i2 As Integer = 0 To Np(0)-1
For i3 As Integer = 0 To N-1
xj(i1) += x(i2,i1,i3)
Next
Next
xj(i1) /= (Np(0) * N)
Next
Dim xij(Np(0),Np(1)) As Double
For i1 As Integer = 0 To Np(0)-1
For i2 As Integer = 0 To Np(1)-1
xij(i1,i2) = 0.0
For i3 As Integer = 0 To N-1
xij(i1,i2) += x(i1,i2,i3)
Next
xij(i1,i2) /= N
Next
Next
Dim xa As Double = 0.0
For i1 As Integer = 0 To Np(0)-1
For i2 As Integer = 0 To Np(1)-1
For i3 As Integer = 0 To N-1
xa += x(i1,i2,i3)
Next
Next
Next
xa /= (Np(0) * Np(1) * N)
Dim SP As Double = 0.0
For i1 As Integer = 0 To Np(0)-1
SP += (xi(i1) - xa) * (xi(i1) - xa)
Next
SP *= (Np(1) * N)
Dim SQ As Double = 0.0
For i1 As Integer = 0 To Np(1)-1
SQ += (xj(i1) - xa) * (xj(i1) - xa)
Next
SQ *= (Np(0) * N)
Dim SI As Double = 0.0
For i1 As Integer = 0 To Np(0)-1
For i2 As Integer = 0 To Np(1)-1
SI += (xij(i1,i2) - xi(i1) - xj(i2) + xa) * (xij(i1,i2) - xi(i1) - xj(i2) + xa)
Next
Next
SI *= N
Dim SE As Double = 0.0
For i1 As Integer = 0 To Np(0)-1
For i2 As Integer = 0 To Np(1)-1
For i3 As Integer = 0 To N-1
SE += (x(i1,i2,i3) - xij(i1,i2)) * (x(i1,i2,i3) - xij(i1,i2))
Next
Next
Next
Dim VP As Double = SP / (Np(0) - 1)
Dim VQ As Double = SQ / (Np(1) - 1)
Dim VI As Double = SI / ((Np(0) - 1) * (Np(1) - 1))
Dim VE As Double = SE / (Np(0) * Np(1) * (N - 1))
Dim FP As Double = VP / VE
Dim FQ As Double = VQ / VE
Dim FI As Double = VI / VE
p = 0.01 * a
dof2 = Np(0) * Np(1) * (N - 1)
Console.WriteLine("全変動： 平方和 " & (SP+SQ+SI+SE) & " 自由度 " & (Np(0)*Np(1)*N-1))
dof1 = Np(0) - 1
Console.WriteLine(name(0) & "の水準間： 平方和 " & SP & " 自由度 " & (Np(0)-1) & " 不偏分散 " & VP & " F " & FP & " " & a & "％値 " & p_f(sw, f_f, f_df, normal_f))
dof1 = Np(1) - 1
Console.WriteLine(name(1) & "の水準間： 平方和 " & SQ & " 自由度 " & (Np(1)-1) & " 不偏分散 " & VQ & " F " & FQ & " " & a & "％値 " & p_f(sw, f_f, f_df, normal_f))
dof1 = (Np(0) - 1) * (Np(1) - 1)
Console.WriteLine("相互作用： 平方和 " & SI & " 自由度 " & ((Np(0)-1)*(Np(1)-1)) & " 不偏分散 " & VI & " F " & FI & " " & a & "％値 " & p_f(sw, f_f, f_df, normal_f))
Console.WriteLine("水準内： 平方和 " & SE & " 自由度 " & (Np(0)*Np(1)*(N-1)) & " 不偏分散 " & VE)
End If
End Sub
'**************************************'
' ｆ分布の計算（P(X = ff), P(X < ff)） '
' dd : P(X = ff) '
' df1,df2 : 自由度 '
' return : P(X < ff) '
'**************************************'
Function f(ff As Double, ByRef dd As Double, df1 As Integer, df2 As Integer)
Dim pp As Double
Dim u As Double
Dim ia As Integer
Dim ib As Integer
If ff < 1.0e-10
ff = 1.0e-10
End If
Dim x As Double = ff * df1 / (ff * df1 + df2)
If (df1 Mod 2) = 0
If (df2 Mod 2) = 0
u = x * (1.0 - x)
pp = x
ia = 2
ib = 2
Else
u = 0.5 * x * Math.Sqrt(1.0-x)
pp = 1.0 - Math.Sqrt(1.0-x)
ia = 2
ib = 1
End If
Else
If (df2 Mod 2) = 0
u = 0.5 * Math.Sqrt(x) * (1.0 - x)
pp = Math.Sqrt(x)
ia = 1
ib = 2
Else
u = Math.Sqrt(x*(1.0-x)) / Math.PI
pp = 1.0 - 2.0 * Math.Atan2(Math.Sqrt(1.0-x), Math.Sqrt(x)) / Math.PI
ia = 1
ib = 1
End If
End If
If ia <> df1
For i1 As Integer = ia To df1-2 Step 2
pp -= 2.0 * u / i1
u *= x * (i1 + ib) / i1
Next
End If
If ib <> df2
For i1 As Integer = ib To df2-2 Step 2
pp += 2.0 * u / i1
u *= (1.0 - x) * (i1 + df1) / i1
Next
End If
dd = u / ff
Return pp
End Function
'**************************************'
' ｆ分布のｐ％値（P(X > u) = 0.01p） '
' ind : >= 0 : normal（収束回数） '
' = -1 : 収束しなかった '
'**************************************'
Function p_f(ByRef ind As Integer, f_f As Func(Of Double, Double),
f_df As Func(Of Double, Double), normal_f As Func(Of Double, Double))
Dim MAX As Integer = 340
Dim a As Double = 0.0
Dim a1 As Double = 0.0
Dim b As Double = 0.0
Dim b1 As Double = 0.0
Dim df1 As Double = 0.0
Dim df2 As Double = 0.0
Dim e As Double = 0.0
Dim ff As Double = 0.0
Dim yq As Double = 0.0
Dim sw As Integer = 0
'
' 初期値計算の準備
' while文は，大きな自由度によるガンマ関数の
' オーバーフローを避けるため
'
Do While sw >= 0
df1 = 0.5 * (dof1 - sw)
df2 = 0.5 * dof2
a = 2.0 / (9.0 * (dof1 - sw))
a1 = 1.0 - a
b = 2.0 / (9.0 * dof2)
b1 = 1.0 - b
yq = p_normal(ind, normal_f)
e = b1 * b1 - b * yq * yq
If e > 0.8 or (dof1+dof2-sw) <= MAX
sw = -1
Else
sw += 1
If (dof1-sw) = 0
sw = -2
End If
End If
Loop
If sw = -2
ind = -1
Else
'
' f0 : 初期値
'
Dim f0 As Double = 0.0
If e > 0.8
Dim x As Double = (a1 * b1 + yq * Math.Sqrt(a1*a1*b+a*e)) / e
f0 = Math.Pow(x, 3.0)
Else
Dim y1 As Double = Math.Pow(dof2, df2-1.0)
Dim y2 As Double = Math.Pow(dof1, df2)
Dim x As Double = gamma(df1+df2, ind) / gamma(df1, ind) /
gamma(df2, ind) * 2.0 * y1 / y2 / p
f0 = Math.Pow(x, 2.0/dof2)
End If
'
' ニュートン法
'
ff = newton(f0, 1.0e-6, 1.0e-10, 100, ind, f_f, f_df)
End If
Return ff
End Function
'**************************************'
' Γ（ｘ）の計算（ガンマ関数，近似式） '
' 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
/*
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 // 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 // 因子の名前 水準の数
3.1 4.7 5.1 // 各水準に対する1番目のデータ
4.1 5.6 3.7 // 各水準に対する2番目のデータ
3.3 4.3 4.5 // 各水準に対する3番目のデータ
3.9 5.9 6.0 // 各水準に対する4番目のデータ
3.7 6.1 3.9 // 各水準に対する5番目のデータ
2.4 4.2 5.4 // 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 // 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 // 1番目の因子の名前 その水準の数
品種 2 // 2番目の因子の名前 その水準の数
3 4 12 -4 -4 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ
*/
| 情報学部 | 菅沼ホーム | 目次 | 索引 |