96 lines
2.3 KiB
C
96 lines
2.3 KiB
C
/* This file contains functions to
|
|
* compute the Landau distribution
|
|
* PDF, CDF and QDF functions.
|
|
*/
|
|
|
|
#include <gsl/gsl_roots.h>
|
|
#include <gsl/gsl_randist.h>
|
|
#include <gsl/gsl_integration.h>
|
|
|
|
|
|
/* This is a wrapper needed by `landau_cdf` because
|
|
* the numerical integration expects a function
|
|
* with parameters.
|
|
*/
|
|
double landau_pdf(double x, void* params) {
|
|
return gsl_ran_landau_pdf(x);
|
|
}
|
|
|
|
|
|
/* The cumulative function of the Landau distribution
|
|
* calculated by numerical integration.
|
|
*/
|
|
double landau_cdf(double x, void* params) {
|
|
// create the integrand
|
|
gsl_function pdf;
|
|
pdf.function = &landau_pdf;
|
|
pdf.params = NULL;
|
|
|
|
// set up the integration
|
|
double res, err;
|
|
size_t iter = 500;
|
|
gsl_integration_workspace* w =
|
|
gsl_integration_workspace_alloc(iter);
|
|
|
|
// clip values too small
|
|
if (x < -200) x = -200;
|
|
|
|
// We integrate the pdf in [x, +∞) instead
|
|
// of (-∞, x] to avoid a singularity and
|
|
// then return 1 - P.
|
|
gsl_integration_qagiu(
|
|
&pdf, // integrand
|
|
x, // lower bound
|
|
1e-10, 1e-6, // abs and rel error
|
|
iter, // max iteration
|
|
w, // "workspace"
|
|
&res, &err); // result, abs error
|
|
|
|
// free the memory
|
|
gsl_integration_workspace_free(w);
|
|
|
|
return 1 - res;
|
|
}
|
|
|
|
|
|
/* Another wrapper: this time it is needed by
|
|
* `landau_qdf` to solve the equation:
|
|
* `landau_cdf(x) = p0` for x
|
|
*/
|
|
double landau_cdf_root(double x, void* params) {
|
|
double p0 = *(double*)params;
|
|
return p0 - landau_cdf(x, NULL);
|
|
}
|
|
|
|
|
|
/* The quantile function (inverse CDF) of the Landau
|
|
* distribution calculated by numerical root method.
|
|
*/
|
|
double landau_qdf(double p0) {
|
|
// create function
|
|
gsl_function cdf;
|
|
cdf.function = &landau_cdf_root;
|
|
cdf.params = &p0;
|
|
|
|
// use the Brent method
|
|
gsl_root_fsolver* s =
|
|
gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
|
|
|
|
// search interval
|
|
double low = -1000, // lower value
|
|
upp = 100000; // upper value
|
|
gsl_root_fsolver_set(s, &cdf, low, upp);
|
|
|
|
// iterative search
|
|
size_t iter = 1000; // max iteration
|
|
int stat = GSL_CONTINUE;
|
|
for (size_t i=0; stat==GSL_CONTINUE && i<iter; i++) {
|
|
stat = gsl_root_fsolver_iterate(s);
|
|
low = gsl_root_fsolver_x_lower(s);
|
|
upp = gsl_root_fsolver_x_upper(s);
|
|
stat = gsl_root_test_interval(low, upp, 0, 0.01);
|
|
}
|
|
|
|
return gsl_root_fsolver_root(s);
|
|
}
|