|
［ 簡単な説明 ］ 一様分布乱数の平均、分散、及び一様さのχ二乗検定です。 ０〜１の一様分布乱数ですから、平均値は 0.5、分散は 0.833333...(=1/12) が理想値です。 １７種類の擬似一様乱数関数を評価しています。 注１ ： χ二乗検定のため、科学技術計算ライブラリの chi( )、qnorm( )、pnorm( )、qchi( )、pchi( )をコピーして使用しています。 注２ ： プログラム中の _init 及び _func は、ヘッダファイル "random.h" 内で、定義されている配列変数で、一様乱数の初期化関数名 及び 生成関数名が、配列で定義されています。 |
（先頭に戻る）
/* test2.c 一様分布乱数テスト */
#include <stdio.h>
#include <stdlib.h>
#include "random.h"
#define ulong unsigned long
void error(char *s, int r);
double chi(int n);
int qnorm(double u,double *qn);
int pnorm(double qn, double *u);
int qchi(double x2, int nu, double *qc);
int pchi(double qc, int nu, double *x2);
#define N1 20 /* 分割数 */
#define N2 50 /* 分割数 */
#define N3 200 /* 分割数 */
static int flag[_M] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
main(int argc, char *argv[])
{
double ave, dn1, dn2, dn3, k1, k2, k3, kk1, kk2, kk3, m, m1, n;
double sigma1, sigma2, sigma3, x, xm, x2;
double Ave, Sigma;
double (*func)(void);
ulong b1[N1], b2[N2], b3[N3];
int i, i1, i2, i3, k;
n = (argc < 2)? 10000.: atof(argv[1]);
fprintf(stderr, "n = %.0f\n", n);
dn1 = (double)N1;
dn2 = (double)N2;
dn3 = (double)N3;
fprintf(stderr, "Now, Calcurate start!\n");
printf("個数 = %.0f １次元分割数 = %d , %d , %d\n", n, N1, N2, N3);
k1 = chi(N1 - 1);
k2 = chi(N2 - 1);
k3 = chi(N3 - 1);
printf(" χ^2\n");
printf(" 関数 平均値 分散 分割%3d %3d %3d\n", N1, N2, N3);
for(k = 0; k < _M; k++)
{
if(flag[k])
{
for(i = 0; i < N1; i++) b1[i] = 0;
for(i = 0; i < N2; i++) b2[i] = 0;
for(i = 0; i < N3; i++) b3[i] = 0;
fprintf(stderr, "***** %s() ", _pname[k]);
Ave = Sigma = m = 0.0;
if(_init[k] != NULL) _init[k](1L);
func = _func[k];
do
{
x = func();
i1 = x * dn1;
i2 = x * dn2;
i3 = x * dn3;
b1[i1]++;
b2[i2]++;
b3[i3]++;
m1 = m++;
x -= Ave;
Ave += (xm = x / m);
Sigma += m1 * x * xm;
} while(m < n);
putc('\n', stderr);
sigma1 = sigma2 = sigma3 = 0.0;
ave = n / dn1;
for(i = 0; i < N1; i++)
{
x = b1[i] - ave;
sigma1 += x * x;
}
ave = n / dn2;
for(i = 0; i < N2; i++)
{
x = b2[i] - ave;
sigma2 += x * x;
}
ave = n / dn3;
for(i = 0; i < N3; i++)
{
x = b3[i] - ave;
sigma3 += x * x;
}
kk1 = sigma1 / n * dn1;
kk2 = sigma2 / n * dn2;
kk3 = sigma3 / n * dn3;
printf("%8s %8.6f %8.6f %6.1f %6.1f %6.1f %s %s %s\n",
_pname[k], Ave, Sigma / (n - 1.), kk1, kk2, kk3,
(kk1 < k1)? " ": "×",
(kk2 < k2)? " ": "×",
(kk3 < k3)? " ": "×");
}
}
fprintf(stderr, "Calcurate End!\n\n");
return 1;
}
void error(char *s, int r)
{
fprintf(stderr, "Error : %s\n", s);
exit(r);
}
double chi(int n)
{
static double p[8] = {0.01, 0.1, 0.5, 1., 2., 5., 10., 20.};
double c;
int j;
printf(" n α =");
for(j = 0; j < 8; j++) printf("%5.2f%% ", p[j]);
putchar('\n');
printf("%4d ", n);
for(j = 0; j < 8; j++)
{
pchi(p[j] / 100., n, &c);
printf("%8.2f", c);
}
putchar('\n');
return c;
}
int qnorm(double u,double *qn)
{
int i;
static double a[9] = { 1.24818987e-4, -1.075204047e-3,
5.198775019e-3, -1.9198292004e-2,
5.9054035642e-2, -1.51968751364e-1,
3.19152932694e-1, -5.319230073e-1,
7.97884560593e-1 },
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,
-1.0557625006e-2, 1.1630447319e-2,
-9.279453341e-3, 5.353579108e-3,
-2.141268741e-3, 5.35310549e-4,
9.99936657524e-1 };
double w, y, z;
if(u == 0.0) z = 0.0;
else
{
y = fabs(u) * 0.5;
if(y >= 3.0) z = 1.0;
else
{
if(y < 1.0)
{
w = y * y;
z = a[0];
for(i = 1; i < 9; i++) z = z * w + a[i];
z *= (y * 2.0);
}
else
{
y -= 2.0;
z = b[0];
for(i = 1; i < 15; i++) z = z * y + b[i];
}
}
}
if(u < 0.0) *qn = (1.0 - z) * 0.5;
else *qn = (1.0 + z) * 0.5;
return 0;
}
/* normal distribution */
int pnorm(double qn, double *u)
{
int i, j;
static double b[11]= { 0.1570796288e1, 0.3706987906e-1,
-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, w2, w3;
if(qn <= 0.0 || qn >= 1.0) return 999;
if(qn == 0.5) *u = 0.0;
else
{
w1 = qn;
if(qn > 0.5) w1 = 1.0 - w1;
w2 = 4.0 * w1 * (1.0 - w1);
w3 = - log(w2);
w1 = b[0];
for(i = 1; i < 11; i++)
{
w2 = i;
w1 += (b[i] * pow(w3, w2));
}
w3 = sqrt(w1 * w3);
if(qn < 0.5) w3 = - w3;
*u = w3;
}
return 0;
}
/* chi-square distribution */
int qchi(double x2, int nu, double *qc)
{
double f, w, pw, x, qw, a, a1, a2;
if(nu < 1) return 999;
if(x2 <= 0.0) *qc = 1.0;
else if(x2 > 400.0) *qc = 0.0;
else if(nu > 10)
{
f = nu;
w = 2.0 / (9.0 * f);
pw = pow(x2 / f, 1.0 / 3.0);
x = (pw - 1.0 + w) / sqrt(w);
qnorm(x, &qw);
*qc = 1.0 - qw;
}
else
{
w = exp(- x2 / 2.0);
x = sqrt(x2);
if((nu % 2) == 0)
{
a1 = 2.0;
*qc = w;
if(nu == 2) return 0;
}
else
{
qnorm(x, &qw);
*qc = 2.0 * (1.0 - qw);
if(nu == 1) return 0;
a1 = 1.0;
w = 0.797884560750774 * w / x;
}
a2 = nu - 2;
for(a = a1; a <= a2; a += 2.0)
{
w *= (x2 / a);
*qc += w;
}
}
return 0;
}
/* chi-square distribution */
int pchi(double qc, int nu, double *x2)
{
int i, j;
double f, u, w, w1, w2, c, c1, c2, c3, gam, x, y, qch, wx;
if(qc <= 0.0 || qc >= 1.0 || nu < 1) return 999;
if(nu == 1)
{
w1 = qc * 0.5;
pnorm(w1, &u);
*x2 = u * u;
}
else if(nu == 2) *x2 = - 2.0 * log(qc);
else
{
f = nu;
pnorm(qc, &u);
u = - u;
if(nu > 10)
{
w1 = sqrt(2.0 * f);
w2 = u * u;
c1 = f + w1 * u + 2.0 * (w2 - 1.0) / 3.0;
c2 = (w2 - 7.0) * u / (9.0 * w1);
c3 = ((3.0 * w2 + 7.0) * w2 - 16.0) * 2.0 / (405.0 * f);
wx = c1 + c2 - c3;
if(wx < 0.0) wx = 0.0;
}
else
{
w = 2.0 / (9.0 * f);
w = 1.0 - w + u * sqrt(w);
wx = f * w * w * w;
if(wx < 0.0) wx = 0.0;
if((nu % 2) == 0) gam = 1.0;
else gam = 1.772453850905516;
j = (nu + 1) / 2 - 1;
w = f / 2.0;
for(i = 1; i <= j; i++)
{
c = i;
gam *= (w - c);
}
x = wx / 2.0;
y = w - 1.0;
c1 = pow(x, y);
c2 = exp(-x) / 2.0;
c3 = c1 * c2;
qchi(wx, nu, &qch);
wx += ((qch - qc) * gam / c3);
}
*x2 = wx;
}
return 0;
}
|
（先頭に戻る）
個数 = 10000 １次元分割数 = 20 , 50 , 200
n α = 0.01% 0.10% 0.50% 1.00% 2.00% 5.00% 10.00% 20.00%
19 50.77 43.81 38.58 36.19 33.69 30.14 27.20 23.90
49 94.59 85.35 78.23 74.92 71.41 66.34 62.04 57.08
199 281.87 266.39 254.13 248.33 242.08 232.91 224.96 215.57
χ^2
関数 平均値 分散 分割 20 50 200 時間[sec]
urnd0 0.504872 0.083440 20.5 46.8 196.5 1.516
urnd1 0.495858 0.084288 23.1 66.4 209.4 1.500 ×
urnd2 0.499196 0.083227 20.2 38.9 183.6 1.516
urnd3 0.503982 0.083137 13.3 46.8 192.9 1.500
urnd4 0.498005 0.084786 24.0 57.9 195.8 1.531 × ×
urnd5 0.497543 0.083469 18.5 36.1 206.9 1.563
urnd6 0.496068 0.082186 17.2 41.8 186.7 1.500
mrnd 0.496310 0.084270 18.5 42.3 210.6 1.531
rrnd 0.498968 0.083952 22.7 41.6 179.0 2.094
ornd0 0.499616 0.083963 25.0 70.5 238.6 1.563 × × ×
ornd1 0.495985 0.082711 24.2 43.9 198.6 1.578 ×
ornd2 0.500468 0.083656 10.0 35.7 203.2 1.563
brnd 0.499566 0.083800 24.7 69.3 237.9 1.672 × × ×
ernd 0.501506 0.082683 21.5 32.1 180.5 1.844
krnd 0.501401 0.084116 27.8 69.0 230.0 1.563 × × ×
drnd 0.505637 0.082814 15.1 44.5 240.2 1.953 ×
qrnd 0.501766 0.083250 14.8 50.5 192.9 1.516
個数 = 100000 １次元分割数 = 20 , 50 , 200
n α = 0.01% 0.10% 0.50% 1.00% 2.00% 5.00% 10.00% 20.00%
19 50.77 43.81 38.58 36.19 33.69 30.14 27.20 23.90
49 94.59 85.35 78.23 74.92 71.41 66.34 62.04 57.08
199 281.87 266.39 254.13 248.33 242.08 232.91 224.96 215.57
χ^2
関数 平均値 分散 分割 20 50 200 時間[sec]
urnd0 0.500248 0.083326 14.9 53.0 186.1 15.672
urnd1 0.500668 0.083258 13.1 41.9 196.1 15.719
urnd2 0.499322 0.083636 18.5 53.3 207.0 15.734
urnd3 0.500502 0.083408 11.8 43.9 196.0 15.813
urnd4 0.501603 0.083232 30.0 82.9 220.2 15.766 × × ×
urnd5 0.499909 0.083222 15.7 35.2 175.3 15.766
urnd6 0.500882 0.083323 16.7 57.1 226.4 15.656 × ×
mrnd 0.500605 0.083535 26.4 45.2 184.7 16.078 ×
rrnd 0.499971 0.083636 13.2 36.7 167.9 21.594
ornd0 0.499398 0.083541 17.3 54.5 210.9 16.469
ornd1 0.498780 0.083248 10.2 27.8 168.7 16.609
ornd2 0.499891 0.083125 17.5 45.3 215.5 16.391
brnd 0.499425 0.083529 17.3 54.3 210.3 17.438
ernd 0.499511 0.083275 23.1 50.4 212.4 19.078
krnd 0.500453 0.083304 9.7 32.8 164.8 16.078
drnd 0.501743 0.083057 21.0 52.6 214.4 20.547
qrnd 0.500728 0.083220 9.5 45.7 179.7 15.703
個数 = 1000000 １次元分割数 = 20 , 50 , 200
n α = 0.01% 0.10% 0.50% 1.00% 2.00% 5.00% 10.00% 20.00%
19 50.77 43.81 38.58 36.19 33.69 30.14 27.20 23.90
49 94.59 85.35 78.23 74.92 71.41 66.34 62.04 57.08
199 281.87 266.39 254.13 248.33 242.08 232.91 224.96 215.57
χ^2
関数 平均値 分散 分割 20 50 200 時間[sec]
urnd0 0.500195 0.083217 11.8 60.7 225.4 157.266 × ×
urnd1 0.499969 0.083301 24.5 51.9 214.2 157.391 ×
urnd2 0.499932 0.083501 19.5 48.8 197.9 158.297
urnd3 0.500188 0.083333 14.3 37.7 153.4 158.625
urnd4 0.500206 0.083253 16.0 52.4 223.1 158.469 ×
urnd5 0.499896 0.083317 8.1 52.0 213.3 158.047
urnd6 0.500172 0.083322 25.9 62.6 238.9 157.453 × × ×
mrnd 0.499962 0.083316 21.8 59.7 203.4 161.391 ×
rrnd 0.500631 0.083335 25.5 61.2 201.3 216.969 × ×
ornd0 0.500256 0.083361 29.0 56.5 204.5 165.781 ×
ornd1 0.499637 0.083363 21.5 54.8 229.1 167.125 ×
ornd2 0.499792 0.083280 9.7 26.5 182.2 164.844
brnd 0.500254 0.083360 29.0 56.3 204.0 175.281 ×
ernd 0.500158 0.083377 8.2 30.6 168.3 192.109
krnd 0.500007 0.083374 10.5 42.9 230.0 161.484 ×
drnd 0.500598 0.083415 25.3 45.6 186.7 206.406 ×
qrnd 0.500116 0.083388 23.4 55.5 234.1 157.891 ×
個数 = 1000000 １次元分割数 = 20 , 50 , 200
n α = 0.01% 0.10% 0.50% 1.00% 2.00% 5.00% 10.00% 20.00%
19 50.77 43.81 38.58 36.19 33.69 30.14 27.20 23.90
49 94.59 85.35 78.23 74.92 71.41 66.34 62.04 57.08
199 281.87 266.39 254.13 248.33 242.08 232.91 224.96 215.57
χ^2
関数 平均値 分散 分割 20 50 200 時間[sec]
urnd0 0.500195 0.083217 11.8 60.7 225.4 148.891 × ×
urnd1 0.499969 0.083301 24.5 51.9 214.2 148.813 ×
urnd2 0.499932 0.083501 19.5 48.8 197.9 149.797
urnd3 0.500188 0.083333 14.3 37.7 153.4 149.609
urnd4 0.500206 0.083253 16.0 52.4 223.1 149.203 ×
urnd5 0.499896 0.083317 8.1 52.0 213.3 149.031
urnd6 0.500172 0.083322 25.9 62.6 238.9 149.578 × × ×
mrnd 0.499962 0.083316 21.8 59.7 203.4 152.203 ×
rrnd 0.499827 0.083391 21.9 56.3 223.7 238.766 ×
ornd0 0.500256 0.083361 29.0 56.5 204.5 159.828 ×
ornd1 0.499637 0.083363 21.5 54.8 229.1 160.094 ×
ornd2 0.499792 0.083280 9.7 26.5 182.2 158.813
brnd 0.500254 0.083360 29.0 56.3 204.0 171.094 ×
ernd 0.500158 0.083377 8.2 30.6 168.3 184.969
krnd 0.500007 0.083374 10.5 42.9 230.0 156.125 ×
drnd 0.500598 0.083415 25.3 45.6 186.7 213.766 ×
qrnd 0.500116 0.083388 23.4 55.5 234.1 150.406 ×
|
戻る（メインメニュー）
戻る（乱数メニュー）
注意事項