F 分布

/****************************/
/* F分布の計算             */
/*      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 f(double, double *, int, int);
double p_f(int *);
double gamma(double, int *);
double newton(double(*)(double), double(*)(double), double, double,
              double, int, int *);
double f_f(double);
double f_df(double);

double p;         // α%値を計算するとき時α/100を設定
int dof1, dof2;   // 自由度
                  // 上の3つの変数は,必ず,この位置で定義しておく必要がある

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

	printf("自由度1は? ");
	scanf("%d", &dof1);
	printf("自由度2は? ");
	scanf("%d", &dof2);
	printf("目的とする結果は? \n");
	printf("     =0 : 確率の計算( P(X = x) 及び P(X < x) の値)\n");
	printf("     =1 : p%値( 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);
			y = f(x, &z, dof1, dof2);
			printf("P(X = %f) = %f,  P(X < %f) = %f (自由度 = %d, %d)\n",
                   x, z, x, y, dof1, dof2);
		}
					// グラフ出力
		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) {
				y = f(x, &z, dof1, dof2);
				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%値 = ∞(自由度 = %d, %d)\n", x, dof1, dof2);
		else if ((1.0-p) < 1.0e-7)
			printf("%f%値 = 0(自由度 = %d, %d)\n", x, dof1, dof2);
		else {
			y = p_f(&sw);
			printf("%f%値 = %f  sw %d (自由度 = %d, %d)\n",
                   x, y, sw, dof1, dof2);
		}
	}

	return 0;
}

/****************************************/
/* f分布の計算(P(X = ff), P(X < ff)) */
/*      dd : P(X = ff)                  */
/*      df1,df2 : 自由度                */
/*      return : P(X < ff)              */
/****************************************/
#include <math.h>

double f(double ff, double *dd, int df1, int df2)
{
	double pi = 4.0 * atan(1.0);
	double pp, u, x;
	int ia, ib, i1;

	if (ff < 1.0e-10)
		ff = 1.0e-10;

	x  = ff * df1 / (ff * df1 + df2);

	if(df1%2 == 0) {
		if (df2%2 == 0) {
			u  = x * (1.0 - x);
			pp = x;
			ia = 2;
			ib = 2;
		}
		else {
			u  = 0.5 * x * sqrt(1.0-x);
			pp = 1.0 - sqrt(1.0-x);
			ia = 2;
			ib = 1;
		}
	}

	else {
		if (df2%2 == 0) {
			u  = 0.5 * sqrt(x) * (1.0 - x);
			pp = sqrt(x);
			ia = 1;
			ib = 2;
		}
		else {
			u  = sqrt(x*(1.0-x)) / pi;
			pp = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi;
			ia = 1;
			ib = 1;
		}
	}

	if (ia != df1) {
		for (i1 = ia; i1 <= df1-2; i1 += 2) {
			pp -= 2.0 * u / i1;
			u  *= x * (i1 + ib) / i1;
		}
	}

	if (ib != df2) {
		for (i1 = ib; i1 <= df2-2; i1 += 2) {
			pp += 2.0 * u / i1;
			u  *= (1.0 - x) * (i1 + df1) / i1;
		}
	}

	*dd = u / ff;

	return pp;
}

/****************************************/
/* f分布のp%値(P(X > u) = 0.01p)   */
/*      ind : >= 0 : normal(収束回数) */
/*            = -1 : 収束しなかった     */
/****************************************/
#include <math.h>

#define MAX 340

double p_f(int *ind)
{
	double a, a1, b, b1, df1, df2, e, ff = 0.0, f0, x, y1, y2, yq;
	int sw = 0;
/*
     初期値計算の準備
          while文は,大きな自由度によるガンマ関数の
          オーバーフローを避けるため
*/
	while (sw >= 0) {

		df1 = 0.5 * (dof1 - sw);
		df2 = 0.5 * dof2;
		a   = 2.0 / (9.0 * (dof1 - sw));
		a1  = 1.0 - a;
		b   = 2.0 / (9.0 * dof2);
		b1  = 1.0 - b;

		yq  = p_normal(ind);

		e   = b1 * b1 - b * yq * yq;

		if (e > 0.8 || (dof1+dof2-sw) <= MAX)
			sw = -1;
		else {
			sw += 1;
			if ((dof1-sw) == 0)
				sw = -2;
		}
	}

	if (sw == -2)
		*ind = -1;

	else {
/*
     f0 : 初期値
*/
		if (e > 0.8) {
			x  = (a1 * b1 + yq * sqrt(a1*a1*b+a*e)) / e;
			f0 = pow(x, 3.0);
		}
		else {
			y1 = pow((double)dof2, df2-1.0);
			y2 = pow((double)dof1, df2);
			x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
                 2.0 * y1 / y2 / p;
			f0 = pow(x, 2.0/dof2);
		}
/*
     ニュートン法
*/
		ff = newton(f_f, f_df, f0, 1.0e-6, 1.0e-10, 100, ind);
	}

	return ff;
}

/****************************************/
/* Γ(x)の計算(ガンマ関数,近似式) */
/*      ier : =0 : normal               */
/*            =-1 : x=-n (n=0,1,2,・・・)  */
/*      return : 結果                   */
/****************************************/
#include <math.h>

double gamma(double x, int *ier)
{
	double err, g, s, t, v, w, y;
	long k;

	*ier = 0;

	if (x > 5.0) {
		v = 1.0 / x;
		s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
            0.00078403922172) * v - 0.000229472093621) * v -
            0.00268132716049) * v + 0.00347222222222) * v +
            0.0833333333333) * v + 1.0;
		g = 2.506628274631001 * exp(-x) * pow(x,x-0.5) * s;
	}

	else {

		err = 1.0e-20;
		w   = x;
		t   = 1.0;

		if (x < 1.5) {

			if (x < err) {
				k = (long)x;
				y = (double)k - x;
				if (fabs(y) < err || fabs(1.0-y) < err)
					*ier = -1;
			}

			if (*ier == 0) {
				while (w < 1.5) {
					t /= w;
					w += 1.0;
				}
			}
		}

		else {
			if (w > 2.5) {
				while (w > 2.5) {
					w -= 1.0;
					t *= w;
				}
			}
		}

		w -= 2.0;
		g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
             0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
             0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
             0.999999926;
		g *= t;
	}

	return g;
}

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

	return f(x, &y, dof1, dof2) - 1.0 + p;
}

/**************************/
/* P(X = x)(関数の微分) */
/**************************/
double f_df(double x)
{
	double y, z;

	z = f(x, &y, dof1, dof2);

	return y;
}

/*****************************************************/
/* Newton法による非線形方程式(f(x)=0)の解            */
/*      fn : f(x)を計算する関数名                    */
/*      dfn : f(x)の微分を計算する関数名             */
/*      x0 : 初期値                                  */
/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
/*      eps2 : 終了条件2(|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;
}

/*************************************************/
/* 標準正規分布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%値(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 : 終了条件1(|x(k+1)-x(k)|<eps1)       */
/*      eps2 : 終了条件2(|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;
}