analistica/ex-1/main.c

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

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


/* Here we generate random numbers following
 * the Landau or the Moyal distribution and run a
 * series of test to check if they seem to belong
 * to the Landau distribution.
 */
int main(int argc, char** argv) {
  size_t samples = 50000;
  char* distr = "lan";
  double m_params[2] = {-0.22278298, 1.1191486};

  /* Process CLI arguments */
  for (size_t i = 1; i < (size_t)argc; i++) {
    if      (!strcmp(argv[i], "-n")) samples = atol(argv[++i]);
    else if (!strcmp(argv[i], "-m")) distr = argv[++i];
    else {
      fprintf(stderr, "Usage: %s -[hnmp]\n", argv[0]);
      fprintf(stderr, "  -h\t\tShow this message.\n");
      fprintf(stderr, "  -n N\t\tThe size of sample to generate. (default: 50000)\n");
      fprintf(stderr, "  -m MODE\tUse Landau 'lan' or Moyal 'moy' distribution. "
          "(default: 'lan')\n");
      return EXIT_FAILURE;
    }
  }

  // initialize an RNG
  gsl_rng_env_setup();
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  // prepare data storage
  double* sample = calloc(samples, sizeof(double));
  double min = -10;
  double max =  10;


  /* Sample generation
   */
  fprintf(stderr, "# Sampling\n");
  fprintf(stderr, "generating %ld points... ", samples);

  /* Sample points from the Landau distribution
   * using the GSL Landau generator or from the
   * Moyal distribution using inverse sampling.
   */
  for(size_t i=0; i < samples; i++){
    if (!strcmp(distr, "lan"))
      sample[i]  = gsl_ran_landau(r);
    if (!strcmp(distr, "moy"))
      sample[i] = moyal_qdf(gsl_rng_uniform(r), m_params);
  }
  fprintf(stderr, "done\n");

  // sort the sample: needed for HSM and KS tests
  qsort(sample, samples, sizeof(double), &cmp_double);


  /* Kolmogorov-Smirnov test
   *
   * Compute the D statistic and its
   * associated probability.
   */
  fprintf(stderr, "\n\n# Kolmogorov-Smirnov test\n");
  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;
  }
  double beta = kolmogorov_cdf(D, samples);

  // print the results
  fprintf(stderr, "\n## Results\n");
  fprintf(stderr, "D=%g\n", D);
  fprintf(stderr, "p=%.3f\n", 1 - beta);


  /* Mode comparison
   *
   * Compute the half-sample mode by bootstrapping
   * and compare the result with the value found by
   * numerical maximisation of the PDF.
   */
  fprintf(stderr, "\n\n# Mode comparison\n");

  /* A structure used by the optimisation
   * routines in numeric_mode and others
   * functions below.
   */
  gsl_function pdf;
  pdf.function = &landau_pdf;
  pdf.params = NULL;

  // number of bootstrap samples
  const size_t boots = 100;

  double mode_e = numeric_mode(min, max, &pdf, 1);
  uncert mode_o = bootstrap_mode(r, sample, samples, boots);

  // print the results
  fprintf(stderr, "\n## Results\n");
  fprintf(stderr, "expected mode: %.8f\n", mode_e);
  fprintf(stderr, "observed mode: %.4f±%.4f\n", mode_o.n, mode_o.s);

  // t-test
  double t = fabs(mode_e - mode_o.n)/mode_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);


  /* FWHM comparison
   *
   * Estimate the FWHM of the sample by constructing
   * an empirical PDF via a KDE method and applying
   * the definition on it (numerical solution of
   * `f(x) = max/2` ⇒ x₁-x₀). This is again bootstrapped
   * to estimate the standard errors and compared against
   * the numerical value of FWHM from the true PDF.
   */
  fprintf(stderr, "\n\n# FWHM comparison\n");

  double fwhm_e = numeric_fwhm(min, max, &pdf, 1);
  uncert fwhm_o = bootstrap_fwhm(r, min, max, sample, samples, boots);

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

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


  /* 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, boots);

  // 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);

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


  /* Infinite moments test
   *
   * Apply the Trapani test for infinite moment
   * to the mean and variance.
   *
   * Use r=n^0.75 points in both cases and rescale
   * the variance by the ψ=2 moment. The mean is not
   * rescaled.
   */
  fprintf(stderr, "\n\n# Infinite moments test\n");
  double p_mean = trapani(r, 0.75, 0, 1, sample, samples);
  double p_var  = trapani(r, 0.75, 1, 2, sample, samples);

  // print the results
  fprintf(stderr, "\n## Results\n");
  fprintf(stderr, "mean: p=%.3f\n", p_mean);
  fprintf(stderr, "variance: p=%.3f\n", p_var);

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