ポアソン分布

/****************************/
/* ポアソン分布の計算       */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>

/**************/
/* 階乗の計算 */
/**************/
double fact(int n)
{
	double x = 1.0;
	int i1;

	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)                   */
/********************************************/
#include <math.h>

double Poisson(int x, double ram, double *pr)
{
	double f = 0.0;
	int i1;

	*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 */
/********/
int main()
{
	double ram, pr, f;
	int x, up, sw;
	char file1[100], file2[100];
	FILE *out1, *out2;

	printf("λ は？ ");
	scanf("%lf", &ram);
	printf("グラフ出力？（=1: yes,  =0: no） ");
	scanf("%d", &sw);
					// 密度関数と分布関数の値
	if (sw == 0) {
		printf("   データは？ ");
		scanf("%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("   密度関数のファイル名は？ ");
		scanf("%s", file1);
		printf("   分布関数のファイル名は？ ");
		scanf("%s", file2);
		printf("   データの上限は？ ");
		scanf("%d", &up);
		out1 = fopen(file1,"w");
		out2 = fopen(file2,"w");
		for (x = 0; x <= up; x++) {
			f = Poisson(x, ram, &pr);
			fprintf(out1, "%d %f\n", x, pr);
			fprintf(out2, "%d %f\n", x, f);
		}
	}

	return 0;
}