]>
Commit | Line | Data |
---|---|---|
ebfd0f31 SH |
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 | |
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 | |
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 | } |