/****************************/
/* 分散分析 */
/* coded by Y.Suganuma */
/****************************/
import java.io.*;
import java.util.*;
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番目のデータ