指数分布

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

double p;     // α％値を計算するとき時α/100を設定
double ram;     // 母数

/****************************************/
/* 指数分布の計算（P(X = x), P(X < x)） */
/*      x : データ                      */
/*      ram : 母数                      */
/*      pr : P(X = x)                   */
/*      return : P(X < x)               */
/****************************************/
#include <math.h>

double exponential(double x, double ram, double *pr)
{
	double f = 0.0, y;

	if (x < 0)
		*pr = 0.0;
	else {
		y   = exp(-ram * x);
		*pr = ram * y;
		f   = 1.0 - y;
	}

	return f;
}

/******************************/
/* P(X > x) - 1 + p（関数値） */
/******************************/
double exponential_f(double x)
{
	return p - exp(-ram * x);
}

/**************************/
/* P(X = x)（関数の微分） */
/**************************/
double exponential_df(double x)
{
	return ram * exp(-ram * x);
}

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

double newton(double(*f)(double), double(*df)(double), double x0,
              double eps1, double eps2, int max, int *ind)
{
	double g, dg, x, x1;
	int sw;

	x1   = x0;
	x    = x1;
	*ind = 0;
	sw   = 0;

	while (sw == 0 && *ind >= 0) {

		sw    = 1;
		*ind += 1;
		g     = (*f)(x1);

		if (fabs(g) > eps2) {
			if (*ind <= max) {
				dg = (*df)(x1);
				if (fabs(dg) > eps2) {
					x = x1 - g / dg;
					if (fabs(x-x1) > eps1 && fabs(x-x1) > eps1*fabs(x)) {
						x1 = x;
						sw = 0;
					}
				}
				else
					*ind = -1;
			}
			else
				*ind = -1;
		}
	}

	return x;
}

/********/
/* main */
/********/
int main()
{
	double x, pr, f, up, h;
	int sw;
	char file1[100], file2[100];
	FILE *out1, *out2;

	printf("母数は？ ");
	scanf("%lf", &ram);
	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);
			f = exponential(x, ram, &pr);
			printf("P(X = %f) = %f,  P( X < %f) = %f (母数 = %f)\n", x, pr, x, f, ram);
		}
					// グラフ出力
		else {
			printf("   密度関数のファイル名は？ ");
			scanf("%s", file1);
			printf("   分布関数のファイル名は？ ");
			scanf("%s", file2);
			out1 = fopen(file1,"w");
			out2 = fopen(file2,"w");
			printf("   データの上限は？ ");
			scanf("%lf", &up);
			printf("   刻み幅は？ ");
			scanf("%lf", &h);
			for (x = 0; x < up+0.5*h; x += h) {
				f = exponential(x, ram, &pr);
				fprintf(out1, "%f %f\n", x, pr);
				fprintf(out2, "%f %f\n", x, f);
			}
		}
	}
					// ％値
	else {
		printf("％の値は？ ");
		scanf("%lf", &x);
		p = 0.01 * x;
		if (p < 1.0e-7)
			printf("%f％値 = ∞ (母数 = %f)\n", x, ram);
		else {
			f = newton(exponential_f, exponential_df, 0.0, 1.0e-6, 1.0e-10, 100, &sw);
			printf("%f％値 = %f  sw %d (母数 = %f)\n", x, f, sw, ram);
		}
	}

	return 0;
}