From 6ad77013b6a0089b0884059037e4b6e62d485d6e Mon Sep 17 00:00:00 2001 From: rnhmjoj Date: Mon, 6 Apr 2020 18:43:00 +0000 Subject: [PATCH] ex-1: generalise numeric_fwhm to take any pdf --- ex-1/main.c | 27 ++++++++++-- ex-1/tests.c | 113 +++++++++++++++++++++++++++++++-------------------- ex-1/tests.h | 13 +++--- 3 files changed, 101 insertions(+), 52 deletions(-) diff --git a/ex-1/main.c b/ex-1/main.c index 3e204ea..f73f1d5 100644 --- a/ex-1/main.c +++ b/ex-1/main.c @@ -72,9 +72,10 @@ int main(int argc, char** argv) { /* Mode comparison * - * Find the bin with the maximum number of events + * 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 @@ -102,8 +103,28 @@ int main(int argc, char** argv) { /* FWHM comparison * - * Find the bins x₋ and x₊. + * 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); + //uncert fwhm_o = bootstrap_fwhm(r, sample, samples, 100); + + // 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 = 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 diff --git a/ex-1/tests.c b/ex-1/tests.c index d2cf017..95d2035 100644 --- a/ex-1/tests.c +++ b/ex-1/tests.c @@ -6,6 +6,7 @@ #include #include #include +#include #include #include "landau.h" @@ -119,67 +120,91 @@ double numeric_mode(double min, double max, } -/* This is the function to be minimized in `numeric_FWHM`. +/* A structure containing the half-maximum + * of a PDF and a gsl_function of the PDF + * itself, used by `numeric_fwhm` to compute + * the FWHM. */ -double abs_landau_pdf(double x, void* params_) { - double* params = ((double *) params_); - return fabs(gsl_ran_landau_pdf(x) - params[0]); +struct fwhm_params { + double halfmax; + gsl_function *pdf; +}; + + +/* This is the implicit equation that is solved + * in `numeric_fwhm` by a numerical root-finding + * method. This function takes a point `x`, the + * parameters defined in `fwhm_params` and returns + * the value of: + * + * f(x) - f(max)/2 + * + * where f is the PDF of interest. + */ +double fwhm_equation(double x, void* params) { + struct fwhm_params p = *((struct fwhm_params*) params); + return p.pdf->function(x, p.pdf->params) - p.halfmax; } -/* Numerically computes the FWHM of Landau - * distribution by maximising the derivative. +/* Numerically computes the FWHM of a PDF + * distribution using the definition. * The `min,max` parameters are the initial search - * interval for the optimisation. `mode` can be - * computer with `numeric_mode(min, max)`. + * interval for the root search of equation + * `fwhm_equation`. Two searches are performed in + * [min, mode] and [mode, max] that will give the + * two solutions x₋, x₊. The FWHM is then x₊-x₋. */ -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; +double numeric_fwhm(double min, double max, + gsl_function *pdf) { + /* Create the gls_function structure + * for the equation to be solved. + */ + double mode = numeric_mode(min, max, pdf); + struct fwhm_params p = { + pdf->function(mode, pdf->params)/2, + pdf + }; + gsl_function equation; + equation.function = &fwhm_equation; + equation.params = &p; + + const gsl_root_fsolver_type * T = gsl_root_fsolver_brent; + gsl_root_fsolver * s = gsl_root_fsolver_alloc(T); // initialize minimization for x₋ - double guess = mode - 1; - double fmin, fmax; - int iter = 0; + double x, fmin, fmax; 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; + gsl_root_fsolver_set(s, &equation, min, mode); + status = GSL_CONTINUE; + for (int iter = 0; status == GSL_CONTINUE && iter < max_iter; iter++) + { status = gsl_root_fsolver_iterate(s); + x = gsl_root_fsolver_root(s); + fmin = gsl_root_fsolver_x_lower(s); + fmax = gsl_root_fsolver_x_upper(s); + status = gsl_min_test_interval(fmin, fmax, 0, prec); + } + double x_low = x; // initialize minimization for x₊ - guess = mode + 1; - gsl_min_fminimizer_set(s, &pdf, guess, mode, max); + gsl_root_fsolver_set(s, &equation, 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); + status = GSL_CONTINUE; + for (int iter = 0; status == GSL_CONTINUE && iter < max_iter; iter++) + { status = gsl_root_fsolver_iterate(s); + x = gsl_root_fsolver_root(s); + fmin = gsl_root_fsolver_x_lower(s); + fmax = gsl_root_fsolver_x_upper(s); + status = gsl_min_test_interval(fmin, fmax, 0, prec); + } + double x_upp = x; - } while (status == GSL_CONTINUE && iter < max_iter); - double x_upp = guess; - - // Free memory - gsl_min_fminimizer_free(s); + // free memory + gsl_root_fsolver_free(s); return x_upp - x_low; } diff --git a/ex-1/tests.h b/ex-1/tests.h index f8f1f6b..9c41cbc 100644 --- a/ex-1/tests.h +++ b/ex-1/tests.h @@ -22,10 +22,13 @@ double numeric_mode(double min, double max, gsl_function *pdf); -/* Numerically computes the FWHM of Landau - * distribution by maximising the derivative. +/* Numerically computes the FWHM of a PDF + * distribution using the definition. * The `min,max` parameters are the initial search - * interval for the optimisation. `mode` can be - * computer with `numeric_mode(min, max)`. + * interval for the root search of equation + * `fwhm_equation`. Two searches are performed in + * [min, mode] and [mode, max] that will give the + * two solutions x₋, x₊. The FWHM is then x₊-x₋. */ -double numeric_fwhm(double min, double max, double mode); +double numeric_fwhm(double min, double max, + gsl_function *pdf);