二項分布

/****************************/
/* 二項分布の計算           */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>

/**************/
/* 組合せ nCr */
/**************/
double comb(int n, int r)
{
	double c = 1.0, x = 1.0, y = 1.0;
	int i1;

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

double binomial(int x, int n, double p, double *pr)
{
	double f = 0.0, q;
	int i1;

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

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

	return 0;
}