/* To see where a sample mean ranks among samples. Normal distribution */ #include #include #include #define n 25 #define ss 100000000 typedef double R; static R v(int j){return j<2 ? j+1 : 3.1415926535897933*v(j-2)*2/j;} // Each invocation of rnd returns a new random normal deviate. R rnd(){static int f = 0; R static const scale = 0x1.p-30; static R m; if(f) {f = 0; return m;} while (1){R x = scale*random()-1, y = scale*random()-1, s = x*x +y*y; if(s<1) {R r = sqrt(- 2*log(s)/s); f = 1; m = x*r; return y*r;}}} int main(){int j = 10000; if(RAND_MAX != 0x7fffffff) {printf("futz\n"); exit(0);} printf("%d samples for the %d-simplex\n", ss, n-1); if(0) while(j--) printf("%20.14f\n", rnd()); if(0){R s1=0, s2=0, s3=0; while (j--){R s = rnd(); s1 += s; s2 += s*s; s3 += s*s*s;} printf("%20.14f %20.14f %20.14f\n", s1, s2, s3);} {int j = ss; int h[n+1]; {int k=n+1; while(k--) h[k]=0;} while(j--){R a[n]; R s = 0; {int k=n; while(k--) {a[k] = rnd(); s += a[k];}} s /= n; {int k=n, m=0; while(k--) if(s