/* 中心極限定理(central limit theorem) 確率変数の和の漸近分布が元の確率変数によらず、正規分布であることを示した。 これのヒストグラムを作成し、正規分布がどうかを確認する。 */ #include#include #include #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define NTAB 32 #define NDIV (1+(IM-1)/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) unsigned long int idum =1024; double mean1(int n); /*平均を求める*/ double mean2(int n); float ran1(long *); /*乱数*/ main() { int i,n,numrnd; double a; numrnd=200; n=30; for(i=1;i<=numrnd;i++) { a=mean2(n); printf("%f\n", a); /*printf("%d %f\n",i, a);*/ } return; } double mean1(int n) { double a,b,c; int i; a=0; for(i=1;i<=n;i++){ b=ran1(&idum); a=a+b; } c=a/n; return c; } double mean2(int n) { double a,b; a=mean1(n); b=sqrt(12.0*n)*(a-0.5); return b; } float ran1(long *idum) { int j; long k; static long iy=0; static long iv[NTAB]; float temp; if(*idum <= 0 || !iy){ if(-(*idum) < 1) *idum =1; else *idum = -(*idum); for (j=NTAB +7;j>=0; j--){ k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; if (*idum < 0) *idum += IM; if( j < NTAB) iv[j]=*idum; } iy=iv[0]; } k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; if (*idum < 0) *idum += IM; j=iy/NDIV; iy=iv[j]; iv[j]= *idum; if((temp=AM*iy) > RNMX) return RNMX; else return temp; }
ヒストグラムを作る
/* ヒストグラムの作成 10-1.txt n=5 10-2.txt n=30 */ #include[出力結果]#include #include #define Num 200 main() { int i,rank,rank2,histo[201],histo2[201]; double a[Num+1],val; FILE *fp; if((fp=fopen("10-1.txt","r") )==NULL) { printf("ファイルが見つかりません. -- 10-1.txt\n"); exit(EXIT_FAILURE); } for (i=0;i<=200;i++){ a[i]=0.0; histo[i]=0; histo2[i]=0; } for (i=0;i =0) histo[rank]++; else { rank2=-rank; histo2[rank2]++; } } for(i=0;i<40;i++){ /*printf("%d %f\n", i, a[i]);*/ /*printf("%d %d\n", i, histo[i]);*/ printf("%d %d\n", -i, histo2[i]); } fclose (fp); return; }
n=5のときのヒストグラム
n=30のときのヒストグラム
比較すると、n=30の方が正規分布に近いことが分かる。
Back to C Language