127 lines
3.1 KiB
C
127 lines
3.1 KiB
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"
|
||
|
||
|
||
/* 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
|
||
*/
|
||
|
||
fprintf(stderr, "\n\n# Mode comparison\n");
|
||
|
||
// print the results
|
||
double mode_e = numeric_mode(min, max);
|
||
uncert mode_o = bootstrap_mode(r, sample, samples, 100);
|
||
fprintf(stderr, "\n## Results\n");
|
||
fprintf(stderr, "expected mode: %.7f\n", mode_e);
|
||
fprintf(stderr, "observed mode: %.4f±%.4f\n", mode_o.n, mode_o.s);
|
||
|
||
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
|
||
*
|
||
* Find the bins x₋ and x₊.
|
||
*/
|
||
|
||
|
||
/* 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);
|
||
|
||
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);
|
||
|
||
|
||
// clean up and exit
|
||
gsl_histogram_free(hist);
|
||
gsl_rng_free(r);
|
||
free(sample);
|
||
return EXIT_SUCCESS;
|
||
}
|