<?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番目のデータ
*/
?>