analistica/ex-1/tests.c

/* This file contains functions to perform
 * statistical tests on the points sampled
 * from a Landau distribution.
 */

#include <stdio.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_sum.h>

#include "landau.h"


/* Kolmogorov distribution CDF
 * for sample size n and statistic D
 */
double kolmogorov_cdf(double D, int n) {
  double x = sqrt(n) * D;

  // trick to reduce estimate error
  x += 1/(6 * sqrt(n)) + (x - 1)/(4 * n);

  // calculate the first n_terms of the series
  // Σ_k=1 exp(-(2k - 1)²π²/8x²)
  size_t n_terms = 30;
  double *terms = calloc(n_terms, sizeof(double));
  for (size_t k=0; k<n_terms; k++) {
    terms[k] = exp(-pow((2*(double)(k + 1) - 1)*M_PI/x, 2) / 8);
  }

  // do a transform to accelerate the convergence
  double sum, abserr;
  gsl_sum_levin_utrunc_workspace* s = gsl_sum_levin_utrunc_alloc(n_terms);
  gsl_sum_levin_utrunc_accel(terms, n_terms, s, &sum, &abserr);

  fprintf(stderr, "\n## Kolmogorov CDF\n");
  fprintf(stderr, "accel sum: %f\n", sum);
  fprintf(stderr, "plain sum: %f\n", s->sum_plain);
  fprintf(stderr, "err: %f\n",   abserr);

  gsl_sum_levin_utrunc_free(s);
  free(terms);

  return sqrt(2*M_PI)/x * sum;
}


/* This is a wrapper needed by `numeric_mode` because
 * the minimization expects a function to be minimized and not
 * maximized.
 */
double neg_landau_pdf(double x, void* params) {
  return (-1) * gsl_ran_landau_pdf(x);
}


/* Numerically computes the mode of a Landau
 * distribution by maximising the derivative.
 * The min,max parameters are the initial search
 * interval for the optimisation.
 */
double numeric_mode(double min, double max) {
  // create function
  gsl_function pdf;
  pdf.function = &neg_landau_pdf;
  pdf.params   = NULL;

  // initialize minimization
  double guess = 0;
  int    iter = 0;
  int    max_iter = 100;
  double prec = 1e-7;
  int    status;
  const gsl_min_fminimizer_type * T = gsl_min_fminimizer_goldensection;
  gsl_min_fminimizer * s = gsl_min_fminimizer_alloc(T);
  gsl_min_fminimizer_set(s, &pdf, guess, min, max);

  // minimization
  do {
      iter++;
      status = gsl_min_fminimizer_iterate(s);
      guess = gsl_min_fminimizer_x_minimum(s);
      min = gsl_min_fminimizer_x_lower(s);
      max = gsl_min_fminimizer_x_upper(s);
      status = gsl_min_test_interval(min, max, prec, prec);

  } while (status == GSL_CONTINUE && iter < max_iter);

  fprintf(stderr, "mode error: %.3g\n", max - min);

  // Free memory
  gsl_min_fminimizer_free(s);
  return guess;
}


/* This is the function to be minimized in `numeric_FWHM`.
 */
double abs_landau_pdf(double x, void* params_) {
  double* params = ((double *) params_);
  return fabs(gsl_ran_landau_pdf(x) - params[0]);
}


/* Numerically computes the FWHM of Landau
 * distribution by maximising the derivative.
 * The `min,max` parameters are the initial search
 * interval for the optimisation. `mode` can be
 * computer with `numeric_mode(min, max)`.
 */
double numeric_fwhm(double min, double max, double mode) {
  // create function
  gsl_function pdf;
  pdf.function = &abs_landau_pdf;
  double params [1];
  params[0]= gsl_ran_landau_pdf(mode)/2;
  pdf.params = params;

  // initialize minimization for x₋
  double guess = mode - 1;
  double fmin, fmax;
  int    iter = 0;
  int    max_iter = 100;
  double prec = 1e-7;
  int    status;
  const gsl_min_fminimizer_type * T = gsl_min_fminimizer_goldensection;
  gsl_min_fminimizer * s = gsl_min_fminimizer_alloc(T);
  gsl_min_fminimizer_set(s, &pdf, guess, min, mode);

  // minimization
  do {
      iter++;
      status = gsl_min_fminimizer_iterate(s);
      guess = gsl_min_fminimizer_x_minimum(s);
      fmin = gsl_min_fminimizer_x_lower(s);
      fmax = gsl_min_fminimizer_x_upper(s);
      status = gsl_min_test_interval(fmin, fmax, prec, prec);
  } while (status == GSL_CONTINUE && iter < max_iter);
  double x_low = guess;

  // initialize minimization for x₊
  guess = mode + 1;
  gsl_min_fminimizer_set(s, &pdf, guess, mode, max);

  // minimization
  do {
      iter++;
      status = gsl_min_fminimizer_iterate(s);
      guess = gsl_min_fminimizer_x_minimum(s);
      fmin = gsl_min_fminimizer_x_lower(s);
      fmax = gsl_min_fminimizer_x_upper(s);
      status = gsl_min_test_interval(fmin, fmax, prec, prec);

  } while (status == GSL_CONTINUE && iter < max_iter);
  double x_upp = guess;

  // Free memory
  gsl_min_fminimizer_free(s);
  return x_upp - x_low;
}