1 /*
2 * Experimental data distribution table generator
3 * Taken from the uncopyrighted NISTnet code (public domain).
4 *
5 * Rread in a series of "random" data values, either
6 * experimentally or generated from some probability distribution.
7 * From this, report statistics.
8 */
9
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <math.h>
13 #include <malloc.h>
14 #include <sys/types.h>
15 #include <sys/stat.h>
16
17 void
stats(FILE * fp)18 stats(FILE *fp)
19 {
20 struct stat info;
21 double *x;
22 int limit;
23 int n=0, i;
24 double mu=0.0, sigma=0.0, sumsquare=0.0, sum=0.0, top=0.0, rho=0.0;
25 double sigma2=0.0;
26
27 fstat(fileno(fp), &info);
28 if (info.st_size > 0) {
29 limit = 2*info.st_size/sizeof(double); /* @@ approximate */
30 } else {
31 limit = 10000;
32 }
33 x = (double *)malloc(limit*sizeof(double));
34
35 for (i=0; i<limit; ++i){
36 fscanf(fp, "%lf", &x[i]);
37 if (feof(fp))
38 break;
39 sumsquare += x[i]*x[i];
40 sum += x[i];
41 ++n;
42 }
43 mu = sum/(double)n;
44 sigma = sqrt((sumsquare - (double)n*mu*mu)/(double)(n-1));
45
46 for (i=1; i < n; ++i){
47 top += ((double)x[i]-mu)*((double)x[i-1]-mu);
48 sigma2 += ((double)x[i-1] - mu)*((double)x[i-1] - mu);
49
50 }
51 rho = top/sigma2;
52
53 printf("mu = %12.6f\n", mu);
54 printf("sigma = %12.6f\n", sigma);
55 printf("rho = %12.6f\n", rho);
56 /*printf("sigma2 = %10.4f\n", sqrt(sigma2/(double)(n-1)));*/
57 /*printf("correlation rho = %10.6f\n", top/((double)(n-1)*sigma*sigma));*/
58 }
59
60
61 int
main(int argc,char ** argv)62 main(int argc, char **argv)
63 {
64 FILE *fp;
65
66 if (argc > 1) {
67 fp = fopen(argv[1], "r");
68 if (!fp) {
69 perror(argv[1]);
70 exit(1);
71 }
72 } else {
73 fp = stdin;
74 }
75 stats(fp);
76 return 0;
77 }
78