analistica/ex-1/main.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_statistics_double.h>

#include "landau.h"
#include "tests.h"
#include "bootstrap.h"


/* Function that compare doubles for sorting:
 * x > y  ⇒ 1
 * x == y ⇒ 0
 * x < y  ⇒ -1
 */
int cmp_double (const void *xp, const void *yp) {
  double x = *(double*)xp,
         y = *(double*)yp;
  return x > y ? 1 : -1;
}


/* Here we generate random numbers in a uniform 
 * range and by using the quantile we map them
 * to a Landau distribution. Then we generate an
 * histogram to check the correctness.
 */
int main(int argc, char** argv) {
  // initialize an RNG
  gsl_rng_env_setup();
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  // prepare histogram
  size_t samples = 100000;
  double* sample = calloc(samples, sizeof(double));
  size_t bins = 40;
  double min = -20;
  double max = 20;
  gsl_histogram* hist = gsl_histogram_alloc(bins);
  gsl_histogram_set_ranges_uniform(hist, min, max);


  /* Sample generation
   *
   * Sample points from the Landau
   * distribution and fill the histogram.
   */
  fprintf(stderr, "# Sampling\n");
  fprintf(stderr, "generating %ld points... ", samples);
  double x;
  for(size_t i=0; i<samples; i++) {
    x = gsl_ran_landau(r);
    sample[i] = x;
    gsl_histogram_increment(hist, x);
  }
  fprintf(stderr, "done\n");

  // sort the sample
  qsort(sample, samples, sizeof(double), &cmp_double);


  /* Kolmogorov-Smirnov test
   *
   * Compute the D statistic and its
   * associated probability.
   */
  double D = 0;
  double d;
  for(size_t i=0; i<samples; i++) {
    d = fabs(landau_cdf(sample[i], NULL) - ((double)i+1)/samples);
    if (d > D)
      D = d;
  }
  fprintf(stderr, "\n\n# Kolmogorov-Smirnov test\n");
  double beta = kolmogorov_cdf(D, samples);

  // print the results
  fprintf(stderr, "\n## Results\n");
  fprintf(stderr, "D: %f\n", D);
  fprintf(stderr, "α: %g\n", 1 - beta);


  /* Mode comparison
   *
   * Find the bin with the maximum number of events
   */
  double mode_o, maxbin = 0;
  double f_mode_o = 0;
  double m1, m2   = 0;
  for(size_t i=0; i<bins; i++) {
    m1 = hist->bin[i];
    if (m1 > m2){
      m2 = m1;
      maxbin = (double)i;
      f_mode_o = hist->bin[i];
    }
  }
  fprintf(stderr, "\n\n# Mode comparison\n");
  fprintf(stderr, "step: %.2f\n ", (max - min)/bins);
  f_mode_o = f_mode_o/samples;
  mode_o = min + (maxbin + 0.5)*(max - min)/bins;

  // print the results
  double mode_e = numeric_mode(min, max);
  fprintf(stderr, "\n## Results\n");
  fprintf(stderr, "expected mode: %.7f\n", mode_e);
  fprintf(stderr, "observed mode: %.3f\n", mode_o);


  /* FWHM comparison
   *
   * Find the bins x₋ and x₊.
   */
  double half = f_mode_o*samples/2;
  m2 = samples;
  double x_low = 0;
  double x_upp = 0;
  double diff;
  for(size_t i=0; i<maxbin; i++) {
    m1 = hist->bin[i];
    diff = fabs(m1 - half);
    if (diff < m2){
      m2 = diff;
      x_low = (double)i;
    }
  }
  m2 = samples;
  for(size_t i=maxbin; i<bins; i++) {
    m1 = hist->bin[i];
    diff = fabs(m1 - half);
    if (diff < m2){
      m2 = diff;
      x_upp = (double)i;
    }
  }
  x_low = min + (x_low + 0.5)*(max - min)/bins;
  x_upp = min + (x_upp + 0.5)*(max - min)/bins;
  double fwhm_o = x_upp - x_low;

  // print the results
  fprintf(stderr, "\n\n# FWHM comparison\n");
  double fwhm_e = numeric_fwhm(min, max, mode_e);
  
  fprintf(stderr, "\n## Results\n");
  fprintf(stderr, "expected FWHM: %.7f\n", fwhm_e);
  fprintf(stderr, "observed FWHM: %.3f\n", fwhm_o);


  /* Median comparison
   *
   * Compute the median of the sample by bootstrapping
   * it and comparing it with the QDF(1/2).
   */
  fprintf(stderr, "\n\n# Median comparison\n");
  double med_e = landau_qdf(0.5);
  uncert med_o = bootstrap_median(r, sample, samples, 100);

  // print the results
  fprintf(stderr, "\n## Results\n");
  fprintf(stderr, "expected median: %.7f\n", med_e);
  fprintf(stderr, "observed median: %.4f±%.4f\n", med_o.n, med_o.s);

  double t = (med_e - med_o.n)/med_o.s;
  double p = 1 - erf(t/sqrt(2));
  fprintf(stderr, "\n## t-test\n");
  fprintf(stderr, "t=%.3f\n", t);
  fprintf(stderr, "p=%.3f\n", p);


  // clean up and exit
  gsl_histogram_free(hist);
  gsl_rng_free(r);
  free(sample);
  return EXIT_SUCCESS;
}