二項分布

<?php

/****************************/
/* 二項分布の計算           */
/*      coded by Y.Suganuma */
/****************************/

/**************/
/* 組合せ nCr */
/**************/
function comb($n, $r)
{
	$c = 1.0;
	$x = 1.0;
	$y = 1.0;

	if ($r > 0) {
		if ($r > $n-$r)
			$r = $n - $r;
		for ($i1 = $n; $i1 >= $n-$r+1; $i1--)
			$x *= $i1;
		for ($i1 = 2; $i1 <= $r; $i1++)
			$y *= $i1;
		$c = $x / $y;
	}

	return $c;
}

/****************************************/
/* 二項分布の計算（P(X = x), P(X < x)） */
/*      x : データ                      */
/*      n,p : パラメータ                */
/*      pr : P(X = x)                   */
/*      return : P(X < x)               */
/****************************************/
function binomial($x, $n, $p, &$pr)
{
	$f  = 0.0;
	$q  = 1.0 - $p;
	$pr = comb($n, $x) * pow($p, $x) * pow($q, $n-$x);
	$f  = $pr;
	for ($i1 = 0; $i1 < $x; $i1++)
		$f += comb($n, $i1) * pow($p, $i1) * pow($q, $n-$i1);

	return $f;
}

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

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

?>