maximum_randN.c File Reference

Functions
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                               {
  if (argc < 3) {
    printf("Usage: 1.This 2.ndata 3.InFile\n");
    exit(EXIT_FAILURE);
  }
  Integer_t i, j, k, n, count, itemp;
  Integer_t const ndata = atoi(argv[1]);
  Integer_t const MAX = 1000;
 
  Integer_t seed=0;
  Integer_t nu = 7;
  double max, maxData, dtemp;
  double rand, rand1, rand2;
  double SampleMean, SampleStddev;
  double CDF[MAX]; for(i = 0;i < MAX; ++i) CDF[i] = 0;
 
  Integer_t subMax;
  Integer_t *histogram, maxCount;
  double mean, stddev, mode;
  double binwidth, binratio=0.01;
 
  Integer_t const BUFSIZE = 1000;
  char str[BUFSIZE];
  char *filename = str;
 
  FILE *Infp;
  strcpy(str, argv[2]);
  Infp = fopen(filename, "r");
  if (Infp == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  skip_header(Infp);
  count=0;
  while( fscanf(Infp, "%lld %lld %le", &i, &j, &dtemp) != EOF ) {
    count+=1;
    CDF[i] = dtemp;
    printf("(count=%lld) %lld %lld %le\n", count, i, j, CDF[i]);
    fflush(stdout);
  }
 
 
  FILE *Outfp;
  sprintf(str, "Dist_maximum_raw.txt");
  Outfp = fopen(filename, "w");
  if (Outfp == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  fprintf(Outfp, "#(ndata=%lld, seed=%lld) Distribution of maximum of |x| \n\n", ndata, seed);
 
  init_genrand(seed);
  maxData = NAN;
  mean    = 0;
  stddev  = 0;
  for(i = 0;i < ndata; ++i) {
    //******************** Choose the number of random variables
    rand1 = genrand_real3();
    for(j = 0;j < MAX; ++j) if(CDF[j] > rand1) break;
    n = (Integer_t)( (j-1)/10.0 );
    // n = 1;
    //******************** (END) Choose the number of random variables
    max = NAN;
    for(j = 0;j < n; ++j) {
      SampleMean = 0;
      SampleStddev = 0;
      for(k = 0;k <= nu; ++k) {
        rand1 = genrand_real3();
        rand2 = genrand_real3();
        rand  = -sqrt(-2.0*log(rand1)) *cos(2*M_PI*rand2);
        SampleMean += fabs(rand);
        SampleStddev += rand*rand;
      }
      SampleMean /= (double)(nu+1);
      SampleStddev /= (double)(nu+1);
      SampleStddev = SampleStddev -SampleMean*SampleMean;
      SampleStddev = sqrt( SampleStddev/nu );
 
      rand = fabs( SampleMean/SampleStddev );
      max = fmax(max, rand);
      // max = rand;
    }
    // subMax  = (Integer_t)(max/binwidth) +1;
    // if(subMax > MAX) {
    //   i -=1;
    //   continue;
    // }
    maxData = fmax(maxData, max);
    mean   += max;
    stddev += max*max;
    fprintf(Outfp, "%lld %le\n", i, max);
  }
  mean    /= (double)ndata;
  stddev  /= (double)ndata;
  stddev   = sqrt( stddev -mean*mean );
  binwidth = binratio*stddev;
  subMax   = (Integer_t)(maxData/binwidth) +1;
  histogram = alloc_ivector(subMax); for(i = 0;i < subMax; ++i) histogram[i] = 0;
 
  sprintf(str, "Dist_maximum_raw.txt");
  Infp = fopen(filename, "r");
  if (Infp == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  skip_header(Infp);
 
  count=0;
  while ( fscanf(Infp, "%lld %le", &i, &max) != EOF ) {
    itemp = (Integer_t)(max/binwidth);
    if(itemp >= subMax) continue;
    if(itemp < 0) continue;
    histogram[itemp] += 1;
    count += 1;
  }
 
  maxCount = 0;
  for(i = 0;i < subMax; ++i) {
    if(maxCount < histogram[i]) {
      mode = binwidth*(i+0.5);
      maxCount  = histogram[i];
    }
  }
 
  sprintf(str, "Dist_maximum.txt");
  Outfp = fopen(filename, "w");
  if (Outfp == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  fprintf(Outfp, "# Input: %s\n", argv[2]);
  fprintf(Outfp, "# mean=%le, stddev=%le, mode=%le\n", mean, stddev, mode);
  fprintf(Outfp, "# 1.Maximum value 2.Probability density\n\n");
  for(i = 0;i < subMax; ++i) {
    fprintf(Outfp, "%le %le\n", binwidth*(i+0.5), histogram[i]/(binwidth*count));
  }
 
  return 0;
}