<?php
/****************************/
/* ｔ分布の計算 */
/* coded by Y.Suganuma */
/****************************/
/*********************************************************/
/* 二分法による非線形方程式(f(x)=0)の解 */
/* f : f(x)を計算する関数名 */
/* x1,x2 : 初期値 */
/* eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜f(x(k))｜＜eps2） */
/* max : 最大試行回数 */
/* ind : 実際の試行回数 */
/* （負の時は解を得ることができなかった） */
/* return : 解 */
/*********************************************************/
function bisection($f, $x1, $x2, $eps1, $eps2, $max, &$ind)
{
$x0 = 0.0;
$f1 = $f($x1);
$f2 = $f($x2);
if ($f1*$f2 > 0.0)
$ind = -1;
else {
$ind = 0;
if ($f1*$f2 == 0.0)
$x0 = ($f1 == 0.0) ? $x1 : $x2;
else {
$sw = 0;
while ($sw == 0 && $ind >= 0) {
$sw = 1;
$ind += 1;
$x0 = 0.5 * ($x1 + $x2);
$f0 = $f($x0);
if (abs($f0) > $eps2) {
if ($ind <= $max) {
if (abs($x1-$x2) > $eps1 && abs($x1-$x2) > $eps1*abs($x2)) {
$sw = 0;
if ($f0*$f1 < 0.0) {
$x2 = $x0;
$f2 = $f0;
}
else {
$x1 = $x0;
$f1 = $f0;
}
}
}
else
$ind = -1;
}
}
}
}
return $x0;
}
/*************************************************/
/* 標準正規分布N(0,1)の計算（P(X = x), P(X < x)）*/
/* w : P(X = x) */
/* return : P(X < x) */
/*************************************************/
function normal($x, &$w)
{
$pi = 4.0 * atan(1.0);
/*
確率密度関数（定義式）
*/
$w = exp(-0.5 * $x * $x) / sqrt(2.0 * $pi);
/*
確率分布関数（近似式を使用）
*/
$y = 0.70710678118654 * abs($x);
$z = 1.0 + $y * (0.0705230784 + $y * (0.0422820123 +
$y * (0.0092705272 + $y * (0.0001520143 + $y * (0.0002765672 +
$y * 0.0000430638)))));
$P = 1.0 - pow($z, -16.0);
if ($x < 0.0)
$P = 0.5 - 0.5 * $P;
else
$P = 0.5 + 0.5 * $P;
return $P;
}
/******************************************************************/
/* 標準正規分布N(0,1)のｐ％値（P(X > u) = 0.01p）（二分法を使用） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/******************************************************************/
function p_normal(&$ind)
{
$u = bisection("normal_f", -7.0, 7.0, 1.0e-6, 1.0e-10, 100, $sw);
$ind = $sw;
return $u;
}
/******************************/
/* 1.0 - p - P(X>x)（関数値） */
/******************************/
function normal_f($x)
{
global $p;
return 1.0 - $p - normal($x, $y);
}
/*****************************************************/
/* Newton法による非線形方程式(f(x)=0)の解 */
/* f : f(x)を計算する関数名 */
/* df : f(x)の微分を計算する関数名 */
/* x0 : 初期値 */
/* eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜f(x(k))｜＜eps2） */
/* max : 最大試行回数 */
/* ind : 実際の試行回数 */
/* （負の時は解を得ることができなかった） */
/* return : 解 */
/*****************************************************/
function newton($f, $df, $x0, $eps1, $eps2, $max, &$ind)
{
$x1 = $x0;
$x = $x1;
$ind = 0;
$sw = 0;
while ($sw == 0 && $ind >= 0) {
$sw = 1;
$ind += 1;
$g = $f($x1);
if (abs($g) > $eps2) {
if ($ind <= $max) {
$dg = $df($x1);
if (abs($dg) > $eps2) {
$x = $x1 - $g / $dg;
if (abs($x-$x1) > $eps1 && abs($x-$x1) > $eps1*abs($x)) {
$x1 = $x;
$sw = 0;
}
}
else
$ind = -1;
}
else
$ind = -1;
}
}
return $x;
}
/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/* ier : =0 : normal */
/* =-1 : x=-n (n=0,1,2,･･･) */
/* return : 結果 */
/****************************************/
function gamma($x, &$ier)
{
$ier = 0;
if ($x > 5.0) {
$v = 1.0 / $x;
$s = ((((((-0.000592166437354 * $v + 0.0000697281375837) * $v +
0.00078403922172) * $v - 0.000229472093621) * $v -
0.00268132716049) * $v + 0.00347222222222) * $v +
0.0833333333333) * $v + 1.0;
$g = 2.506628274631001 * exp(-$x) * pow($x,$x-0.5) * $s;
}
else {
$err = 1.0e-20;
$w = $x;
$t = 1.0;
if ($x < 1.5) {
if ($x < $err) {
$k = intval($x);
$y = floatval($k) - $x;
if (abs($y) < $err || abs(1.0-$y) < $err)
$ier = -1;
}
if ($ier == 0) {
while ($w < 1.5) {
$t /= $w;
$w += 1.0;
}
}
}
else {
if ($w > 2.5) {
while ($w > 2.5) {
$w -= 1.0;
$t *= $w;
}
}
}
$w -= 2.0;
$g = (((((((0.0021385778 * $w - 0.0034961289) * $w +
0.0122995771) * $w - 0.00012513767) * $w + 0.0740648982) * $w +
0.0815652323) * $w + 0.411849671) * $w + 0.422784604) * $w +
0.999999926;
$g *= $t;
}
return $g;
}
/**************************************************/
/* ｔ分布の計算（P(X = tt), P(X < tt)） */
/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
/* dd : P(X = tt) */
/* df : 自由度 */
/* return : P(X < tt) */
/**************************************************/
function t($tt, &$dd, $df)
{
$pi = 4.0 * atan(1.0);
$sign = ($tt < 0.0) ? -1.0 : 1.0;
if (abs($tt) < 1.0e-10)
$tt = $sign * 1.0e-10;
$t2 = $tt * $tt;
$x = $t2 / ($t2 + $df);
if($df%2 != 0) {
$u = sqrt($x*(1.0-$x)) / $pi;
$p = 1.0 - 2.0 * atan2(sqrt(1.0-$x), sqrt($x)) / $pi;
$ia = 1;
}
else {
$u = sqrt($x) * (1.0 - $x) / 2.0;
$p = sqrt($x);
$ia = 2;
}
if ($ia != $df) {
for ($i1 = $ia; $i1 <= $df-2; $i1 += 2) {
$p += 2.0 * $u / $i1;
$u *= (1.0 + $i1) / $i1 * (1.0 - $x);
}
}
$dd = $u / abs($tt);
$pp = 0.5 + 0.5 * $sign * $p;
return $pp;
}
/**************************************************/
/* ｔ分布のｐ％値（P(X > u) = 0.01p） */
/* （自由度が∞の時の値はN(0,1)を利用して下さい） */
/* ind : >= 0 : normal（収束回数） */
/* = -1 : 収束しなかった */
/**************************************************/
function p_t(&$ind)
{
global $p, $dof;
$pi = 4.0 * atan(1.0);
$tt = 0.0;
$pis = sqrt($pi);
$df = floatval($dof);
$df2 = 0.5 * $df;
// 自由度＝１
if ($dof == 1)
$tt = tan($pi*(0.5-$p));
else {
// 自由度＝２
if ($dof == 2) {
$c = ($p > 0.5) ? -1.0 : 1.0;
$p2 = (1.0 - 2.0 * $p);
$p2 *= $p2;
$tt = $c * sqrt(2.0 * $p2 / (1.0 - $p2));
}
// 自由度＞２
else {
$yq = p_normal($ind); // 初期値計算のため
if ($ind >= 0) {
$x = 1.0 - 1.0 / (4.0 * $df);
$e = $x * $x - $yq * $yq / (2.0 * $df);
if ($e > 0.5)
$t0 = $yq / sqrt($e);
else {
$x = sqrt($df) / ($pis * $p * $df * gamma($df2, $ind) / gamma($df2+0.5, $ind));
$t0 = pow($x, 1.0/$df);
}
// ニュートン法
$tt = newton("t_f", "t_df", $t0, 1.0e-6, 1.0e-10, 100, $ind);
}
}
}
return $tt;
}
/********************************/
/* 1.0 - p - P(X > x)（関数値） */
/********************************/
function t_f($x)
{
global $p, $dof;
return t($x, $y, $dof) - 1.0 + $p;
}
/**************************/
/* P(X = x)（関数の微分） */
/**************************/
function t_df($x)
{
global $dof;
$z = t($x, $y, $dof);
return $y;
}
/********/
/* main */
/********/
printf("自由度は？ ");
fscanf(STDIN, "%d", $dof);
printf("目的とする結果は？ \n");
printf(" =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n");
printf(" =1 : ｐ％値（ P(X > u) = 0.01p となるuの値） ");
fscanf(STDIN, "%d", $sw);
if ($sw == 0) {
printf("グラフ出力？（=1: yes, =0: no） ");
fscanf(STDIN, "%d", $sw);
// 密度関数と分布関数の値
if ($sw == 0) {
printf(" データは？ ");
fscanf(STDIN, "%lf", $x);
$y = t($x, $z, $dof);
printf("P(X = %f) = %f, P( X < %f) = %f (自由度 = %d)\n", $x, $z, $x, $y, $dof);
}
// グラフ出力
else {
printf(" 密度関数のファイル名は？ ");
fscanf(STDIN, "%s", $file1);
printf(" 分布関数のファイル名は？ ");
fscanf(STDIN, "%s", $file2);
$out1 = fopen($file1, "wb");
$out2 = fopen($file2, "wb");
printf(" データの下限は？ ");
fscanf(STDIN, "%lf", $x1);
printf(" データの上限は？ ");
fscanf(STDIN, "%lf", $x2);
printf(" 刻み幅は？ ");
fscanf(STDIN, "%lf", $h);
for ($x = $x1; $x < $x2+0.5*$h; $x += $h) {
$y = t($x, $z, $dof);
fwrite($out1, $x." ".$z."\n");
fwrite($out2, $x." ".$y."\n");
}
}
}
// ％値
else {
printf("％の値は？ ");
fscanf(STDIN, "%lf", $x);
$p = 0.01 * $x;
if ($p < 1.0e-7)
printf("%f％値 = ∞ (自由度 = %d)\n", $x, $dof);
else if ((1.0-$p) < 1.0e-7)
printf("%f％値 = -∞ (自由度 = %d)\n", $x, $dof);
else {
$y = p_t($sw);
printf("%f％値 = %f sw %d (自由度 = %d)\n", $x, $y, $sw, $dof);
}
}
?>