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
17void
18stats(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
61int
62main(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