正規分布

<?php

/****************************/
/* 正規分布の計算           */
/*      coded by Y.Suganuma */
/****************************/

	printf("平均値は？ ");
	fscanf(STDIN, "%lf", $mean);
	printf("標準偏差は？ ");
	fscanf(STDIN, "%lf", $sd);
	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);
			$xx = ($x - $mean) / $sd;
			$y  = normal($xx, $z);
			printf("P(X = %f) = %f,  P(X < %f) = %f\n", $x, $z, $x, $y);
		}
					// グラフ出力
		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) {
				$xx = ($x - $mean) / $sd;
				$y  = normal($xx, $z);
				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％値 = ∞\n", $x);
		else if ((1.0-$p)< 1.0e-7)
			printf("%f％値 = -∞\n", $x);
		else {
			$y = $sd * p_normal($sw) + $mean;
			printf("%f％値 = %f  sw %d\n", $x, $y, $sw);
		}
	}

/*************************************************/
/* 標準正規分布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);
}

/*********************************************************/
/* 二分法による非線形方程式(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;
}

?>