正規分布

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

double normal(double, double *);
double p_normal(int *);
double bisection(double(*)(double), double, double, double, double, int, int *);
double normal_f(double);

double p;     // α％値を計算するとき時α/100を設定
              // 必ず，この位置で定義しておく必要がある

int main()
{
	double h, x, xx, x1, x2, y, z, mean, sd;
	int sw;
	char file1[100], file2[100];
	FILE *out1, *out2;

	printf("平均値は？ ");
	scanf("%lf",&mean);
	printf("標準偏差は？ ");
	scanf("%lf",&sd);
	printf("目的とする結果は？ \n");
	printf("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n");
	printf("     =1 : ｐ％値（ P(X > u) = 0.01p となるuの値） ");
	scanf("%d", &sw);

	if (sw == 0) {

		printf("グラフ出力？（=1: yes,  =0: no） ");
		scanf("%d", &sw);
					// 密度関数と分布関数の値
		if (sw == 0) {
			printf("   データは？ ");
			scanf("%lf", &x);
			xx = (x - mean) / sd;
			y  = normal(xx, &z);
			printf("P(X = %f) = %f,  P(X < %f) = %f\n", x, z, x, y);
		}
					// グラフ出力
		else {
			printf("   密度関数のファイル名は？ ");
			scanf("%s", file1);
			printf("   分布関数のファイル名は？ ");
			scanf("%s", file2);
			out1 = fopen(file1,"w");
			out2 = fopen(file2,"w");
			printf("   データの下限は？ ");
			scanf("%lf", &x1);
			printf("   データの上限は？ ");
			scanf("%lf", &x2);
			printf("   刻み幅は？ ");
			scanf("%lf", &h);
			for (x = x1; x < x2+0.5*h; x += h) {
				xx = (x - mean) / sd;
				y  = normal(xx, &z);
				fprintf(out1, "%f %f\n", x, z);
				fprintf(out2, "%f %f\n", x, y);
			}
		}
	}
					// ％値
	else {
		printf("％の値は？ ");
		scanf("%lf", &x);
		p = 0.01 * x;
		if (p < 1.0e-7)
			printf("%f％値 = ∞\n", x);
		else if ((1.0-p)< 1.0e-7)
			printf("%f％値 = -∞\n", x);
		else {
			y = sd * p_normal(&sw) + mean;
			printf("%f％値 = %f  sw %d\n", x, y, sw);
		}
	}

	return 0;
}

/*************************************************/
/* 標準正規分布N(0,1)の計算（P(X = x), P(X < x)）*/
/*      w : P(X = x)                             */
/*      return : P(X < x)                        */
/*************************************************/
#include <math.h>

double normal(double x, double *w)
{
	double pi = 4.0 * atan(1.0);
	double y, z, P;
/*
     確率密度関数（定義式）
*/
	*w = exp(-0.5 * x * x) / sqrt(2.0*pi);
/*
     確率分布関数（近似式を使用）
*/
	y = 0.70710678118654 * fabs(x);
	z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
        y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
        y * 0.0000430638)))));
	P = 1.0 - pow(z,-16.0);

	if (x < 0.0)
		P = 0.5 - 0.5 * P;
	else
		P = 0.5 + 0.5 * P;

	return P;
}

/******************************************************************/
/* 標準正規分布N(0,1)のｐ％値（P(X > u) = 0.01p）（二分法を使用） */
/*      ind : >= 0 : normal（収束回数）                           */
/*            = -1 : 収束しなかった                               */
/******************************************************************/
double p_normal(int *ind)
{
	double u;
	int sw;

	u    = bisection(normal_f, -7.0, 7.0, 1.0e-6, 1.0e-10, 100, &sw);
	*ind = sw;

	return u;
}

/******************************/
/* 1.0 - p - P(X>x)（関数値） */
/******************************/
double normal_f(double x)
{
	double y;

	return 1.0 - p - normal(x, &y);
}

/*********************************************************/
/* 二分法による非線形方程式(f(x)=0)の解                  */
/*      fn : f(x)を計算する関数名                        */
/*      x1,x2 : 初期値                                   */
/*      eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
/*      eps2 : 終了条件２（｜f(x(k))｜＜eps2）           */
/*      max : 最大試行回数                               */
/*      ind : 実際の試行回数                             */
/*            （負の時は解を得ることができなかった）     */
/*      return : 解                                      */
/*********************************************************/
#include <math.h>

double bisection(double(*f)(double), double x1, double x2,
                 double eps1, double eps2, int max, int *ind)
{
	double f0, f1, f2, x0 = 0.0;
	int sw;

	f1 = (*f)(x1);
	f2 = (*f)(x2);

	if (f1*f2 > 0.0)
		*ind = -1;

	else {
		*ind = 0;
		if (f1*f2 == 0.0)
			x0 = (f1 == 0.0) ? x1 : x2;
		else {
			sw = 0;
			while (sw == 0 && *ind >= 0) {
				sw    = 1;
				*ind += 1;
				x0    = 0.5 * (x1 + x2);
				f0    = (*f)(x0);

				if (fabs(f0) > eps2) {
					if (*ind <= max) {
						if (fabs(x1-x2) > eps1 && fabs(x1-x2) > eps1*fabs(x2)) {
							sw = 0;
							if (f0*f1 < 0.0) {
								x2 = x0;
								f2 = f0;
							}
							else {
								x1 = x0;
								f1 = f0;
							}
						}
					}
					else
						*ind = -1;
				}
			}
		}
	}

	return x0;
}