指数分布

<?php

/****************************/
/* 指数分布の計算           */
/*      coded by Y.Suganuma */
/****************************/

/****************************************/
/* 指数分布の計算（P(X = x), P(X < x)） */
/*      x : データ                      */
/*      ram : 母数                      */
/*      pr : P(X = x)                   */
/*      return : P(X < x)               */
/****************************************/

function exponential($x, $ram, &$pr)
{
	$f = 0.0;

	if ($x < 0)
		$pr = 0.0;
	else {
		$y  = exp(-$ram * $x);
		$pr = $ram * $y;
		$f  = 1.0 - $y;
	}

	return $f;
}

/******************************/
/* P(X > x) - 1 + p（関数値） */
/******************************/
function exponential_f($x)
{
	global $p, $ram;
	return $p - exp(-$ram * $x);
}

/**************************/
/* P(X = x)（関数の微分） */
/**************************/
function exponential_df($x)
{
	global $ram;
	return $ram * exp(-$ram * $x);
}

/*****************************************************/
/* 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;
}

/********/
/* main */
/********/

	printf("母数は？ ");
	fscanf(STDIN, "%lf", $ram);
	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);
			$f = exponential($x, $ram, $pr);
			printf("P(X = %f) = %f,  P( X < %f) = %f (母数 = %f)\n", $x, $pr, $x, $f, $ram);
		}
					// グラフ出力
		else {
			printf("   密度関数のファイル名は？ ");
			fscanf(STDIN, "%s", $file1);
			printf("   分布関数のファイル名は？ ");
			fscanf(STDIN, "%s", $file2);
			$out1 = fopen($file1, "wb");
			$out2 = fopen($file2, "wb");
			printf("   データの上限は？ ");
			fscanf(STDIN, "%lf", $up);
			printf("   刻み幅は？ ");
			fscanf(STDIN, "%lf", $h);
			for ($x = 0; $x < $up+0.5*$h; $x += $h) {
				$f = exponential($x, $ram, $pr);
				fwrite($out1, $x." ".$pr."\n");
				fwrite($out2, $x." ".$f."\n");
			}
		}
	}
					// ％値
	else {
		printf("％の値は？ ");
		fscanf(STDIN, "%lf", $x);
		$p = 0.01 * $x;
		if ($p < 1.0e-7)
			printf("%f％値 = ∞ (母数 = %f)\n", $x, $ram);
		else {
			$f = newton("exponential_f", "exponential_df", 0.0, 1.0e-6, 1.0e-10, 100, $sw);
			printf("%f％値 = %f  sw %d (母数 = %f)\n", $x, $f, $sw, $ram);
		}
	}

?>