LXR iproute2/netem/stats.c

   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