ポアソン分布

<?php

/****************************/
/* ポアソン分布の計算       */
/*      coded by Y.Suganuma */
/****************************/

/**************/
/* 階乗の計算 */
/**************/
function fact($n)
{
	$x = 1.0;

	for ($i1 = 2; $i1 <= $n; $i1++)
		$x *= $i1;

	return $x;
}
/********************************************/
/* ポアソン分布の計算（P(X = x), P(X < x)） */
/*      x : データ                          */
/*      ram : パラメータ                    */
/*      pr : P(X = x)                       */
/*      return : P(X < x)                   */
/********************************************/
function Poisson($x, $ram, &$pr)
{
	$f = 0.0;

	$pr = pow($ram, $x) * exp(-$ram) / fact($x);
	$f  = $pr;
	for ($i1 = 0; $i1 < $x; $i1++)
		$f += pow($ram, $i1) * exp(-$ram) / fact($i1);

	return $f;
}

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

	printf("λ は？ ");
	fscanf(STDIN, "%lf", $ram);
	printf("グラフ出力？（=1: yes,  =0: no） ");
	fscanf(STDIN, "%d", $sw);
					// 密度関数と分布関数の値
	if ($sw == 0) {
		printf("   データは？ ");
		fscanf(STDIN, "%d", $x);
		$f = Poisson($x, $ram, $pr);
		printf("P(X = %d) = %f,  P( X < %d) = %f (λ = %f)\n", $x, $pr, $x, $f, $ram);
	}
					// グラフ出力
	else {
		printf("   密度関数のファイル名は？ ");
		fscanf(STDIN, "%s", $file1);
		printf("   分布関数のファイル名は？ ");
		fscanf(STDIN, "%s", $file2);
		printf("   データの上限は？ ");
		fscanf(STDIN, "%d", $up);
		$out1 = fopen($file1, "wb");
		$out2 = fopen($file2, "wb");
		for ($x = 0; $x <= $up; $x++) {
			$f = Poisson($x, $ram, $pr);
			fwrite($out1, $x." ".$pr."\n");
			fwrite($out2, $x." ".$f."\n");
		}
	}

?>