乱数プログラム 使用例３

/*	test3.c		正規分布乱数の検定(Random number of Normal Distribution)	*/

#include <stdio.h>
#include <stdlib.h>
#include "random.h"

#define		N		50			/* 分割数 - 2 */

double qnorm(double u);
double pnorm(double qn);
double qchi(double x2, int n);
double pchi(double qc, int n);

main(int argc, char *argv[])
{
	double av, st, ave, du, m, m1, q, q0, u, sigma, x, x2, x2a, x2b, x2c, xm, w;
	double p[N + 2];
	int i, i1, i2, j, k, loop, n, nmax = N + 1, n10;
	int b[N + 2], f[N + 2];

	n = (argc < 2)? 10000: atoi(argv[1]);
	if (n < 10)	n = 10000;
	n10 = n / 10;

	av = (double)N / 2. + 1.;
	st = (double)N / 10.;
	du = 1. / st;

	for(i = N / 2 + 1, u = du, q0 = 0.5; i <= N; i++)
	{
		q = qnorm(u);
		p[i] = q - q0;
		q0 = q;
		u += du;
	}
	p[i] = 1. - q;
	for(i = N / 2 + 1; i <= N + 1; i++)	p[N + 1 - i] = p[i]; 
	printf("個数 = %10d                             棄却水準 χo^2\n", n);

	printf("関数\t平均値      標準偏差       χ^2      [  1%%      0.1%%    0.01%% ]\n");

	for(loop = 0; loop < _M; loop++)
	{
		fprintf(stderr, "*** %s  start   ", _pname[loop]);
 
		init_gaussd(_init[loop], _func[loop], 0);

		ave = sigma = m = 0.0;
		for(i = 0; i <= N + 1; i++)	b[i] = 0.0;
		while(m < n)
		{
			x = gaussd(0., 1.);
			i = (int)(x * st + av);
			if(i < 0)	i = 0;
			else if(i > nmax)	i = nmax;
			b[i]++;
			m1 = m++;
			ave += (xm = (x -= ave) / m);
			sigma += m1 * x * xm;
			if((int)m % n10 == 0)	putc('.', stderr);
		}

		fprintf(stderr, "  test start");

/*	検定	*/
		for(i = 0; i <= N + 1; i++)	f[i] = (int)((double)n * p[i]);
		i1 = 0;
		i2 = N + 1;
		for(i = N / 2; i >= 0; i--)
		{
			if(b[i] < 5)
			{
				i1 = i + 1;
				for(j = k = 0; i >= 0; i--)
				{
					j += b[i];
					k += f[i];
					b[i] = f[i] = 0;
				}
				b[i1] += j;
				f[i1] += k;
				break;
			}
		}

		for(i = N / 2 + 1; i <= N + 1; i++)
		{
			if(b[i] < 5)
			{
				i2 = i - 1;
				for(j = k = 0; i <= N + 1; i++)
				{
					j += b[i];
					k += f[i];
					b[i] = f[i] = 0;
				}
				b[i2] += j;
				f[i2] += k;
				break;
			}
		}
		x2 = 0.;
		for(i = i1; i <= i2; i++)
		{
			w = (double)(b[i] - f[i]);
			x2 += (w * w / (double)f[i]);
		}

		putc('\n', stderr);
		x2a = pchi(0.01, i2 - i1);
		x2b = pchi(0.001, i2 - i1);
		x2c = pchi(0.0001, i2 - i1);
		printf("%s\t%11.8f %11.8f  (%8.2f)  [ %6.2f  %6.2f  %6.2f ]\n",
			_pname[loop], ave, sqrt(sigma / m1), x2, x2a, x2b, x2c);
	}
	fprintf(stderr, "***** end of test3 *****\n");

	return 1;
}

double qnorm(double u)
{
	static double a[9] = {	 1.24818987e-4, -1.075204047e-3, 5.198775019e-3,
							-0.019198292004, 0.059054035642,-0.151968751364,
							 0.319152932694,-0.5319230073,   0.797884560593};
	static double b[15] = {	-4.5255659e-5,   1.5252929e-4,  -1.9538132e-5,
							-6.76904986e-4,  1.390604284e-3,-7.9462082e-4,
							-2.034254874e-3, 6.549791214e-3,-0.010557625006,
							 0.011630447319,-9.279453341e-3, 5.353579108e-3,
							-2.141268741e-3, 5.35310549e-4,  0.999936657524};
	double w, y, z;
	int i;

	if(u == 0.)	return 0.5;
	y = u / 2.;
	if(y < -6.)	return 0.;
	if(y > 6.)		return 1.;
	if(y < 0.)		y = - y;
	if(y < 1.)
	{
		w = y * y;
		z = a[0];
		for(i = 1; i < 9; i++)		z = z * w + a[i];
		z *= (y * 2.);
	}
	else
	{
		y -= 2.;
		z = b[0];
		for(i = 1; i < 15; i++)	z = z * y + b[i];
	}

	if(u < 0.)	return (1. - z) / 2.;
	return (1. + z) / 2.;
}

double pnorm(double qn)
{
	static double b[11] = {	 1.570796288,     0.03706987906,  -0.8364353589e-3,
							-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
							-0.104527497e-5,  0.8360937017e-7,-0.3231081277e-8,
							 0.3657763036e-10,0.6936233982e-12};
	double w1, w3;
	int i;

	if(qn < 0. || 1. < qn)
	{
		fprintf(stderr, "Error : qn <= 0 or qn >= 1  in pnorm()!\n");
		return 0.;
	}
	if(qn == 0.5)	return 0.;

	w1 = qn;
	if(qn > 0.5)	w1 = 1. - w1;
	w3 = -log(4. * w1 * (1. - w1));
	w1 = b[0];
	for(i = 1; i < 11; i++)	w1 += (b[i] * pow(w3, (double)i));
	if(qn > 0.5)	return sqrt(w1 * w3);
	return -sqrt(w1 * w3);
}
 
double qchi(double x2, int n)
{
	double w, pw, x, qc;
	int i, i1;

	if(n < 1)
	{
		fprintf(stderr,"Error : 自由度 < 1 in qchi()!\n");
		return 0.;
	}
	if(x2 <= 0.)	return 1.;
	if(x2 > 400.)	return 0.;
	if(n > 10)
	{
		w = 2. / (9. * (double)n);
		pw = pow(x2 / (double)n, 1. / 3.);
		return 1. - qnorm((pw - 1. + w) / sqrt(w));
	}
	w = exp(-x2 / 2.);
	if(n == 2)	return w;
	x = sqrt(x2);
	if(n == 1)	return 2. * (1. - qnorm(x));
	if((n % 2) == 0)
	{
		i1 = 2;
		qc = w;
	}
	else
	{
 		i1 = 1;
		qc = 2. * (1. - qnorm(x));
		w *= (0.797884560750774 / x);
	}
	for(i = i1; i <= n - 2; i += 2)
	{
		w *= (x2 / (double)i);
		qc += w;
	}
	return qc;
}

double pchi(double qc, int n)
{
	double c1, c2, gam, x, w, wx;
	int i, j;

	if(qc <= 0. || qc >= 1. || n < 1)
	{
		fprintf(stderr,"Error : Illigal parameter in pchi()!\n");
		return 0.;
	}
	if(n == 1)
	{
		w = pnorm(qc / 2.);
		return w * w;
	}
	if(n == 2)	return -2. * log(qc);

	x = -pnorm(qc);
	if(n > 10)
	{
		w = x * x;
		wx = sqrt(2. * (double)n);
		c1 = (w - 7.) * x / 9. / wx;
		c2 = ((3. * w + 7.) * w - 16.) * 2. / 405. / (double)n;
		wx = (double)n + wx * x + 2. * (w - 1.) / 3. + c1 - c2;
		if(wx < 0.)	return 0.;
		return wx;
	}

	w = 2. / 9. / (double)n;
	w = 1. - w + x * sqrt(w);
	wx = (double)n * w * w * w;
	if((n % 2) == 0)	gam = 1.;
	else				gam = 1.772453850905516;
	j = (n + 1) / 2 - 1;
	w = (double)n / 2.;
	for(i = 1; i <= j; i++)	gam *= (w - (double)i);
	x = wx / 2.;
	c1 = pow(x, w - 1.);
	c2 = exp(-x) / 2.;
	return wx + (qchi(wx, n) - qc) * gam / c1 / c2;
}