Stephen Hemminger | ebfd0f3 | 2006-06-15 16:22:09 -0700 | [diff] [blame] | 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 | } |