Giu/analistica
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

ex-1: add the possibility of sampling points from Moyal PDF

Browse Source
This commit is contained in:
Giù Marcer 2020-06-08 18:03:04 +02:00 committed by rnhmjoj
parent 87d506ea50
commit 297d9dacf6
1 changed files with 28 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
ex-1/main.c
View File

@ -1,30 +1,35 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <math.h>
#include <gsl/gsl_randist.h> #include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics_double.h> #include <gsl/gsl_statistics_double.h>
#include "landau.h" #include "landau.h"
#include "moyal.h"
#include "tests.h" #include "tests.h"
#include "bootstrap.h" #include "bootstrap.h"
/* Here we generate random numbers following /* Here we generate random numbers following
* the Landau distribution and run a series of * the Landau or the Moyal distribution and run a
* test to check if they really belong to such a * series of test to check if they seem to belong
* distribution. * to the Landau distribution.
*/ */
int main(int argc, char** argv) { int main(int argc, char** argv) {
size_t samples = 50000; size_t samples = 50000;
char* distr = "lan";
double m_params [2] = {-0.22278298, 1.1191486};
/* Process CLI arguments */ /* Process CLI arguments */
for (size_t i = 1; i < argc; i++) { for (size_t i = 1; i < argc; i++) {
if (!strcmp(argv[i], "-n")) samples = atol(argv[++i]); if (!strcmp(argv[i], "-n")) samples = atol(argv[++i]);
else if (!strcmp(argv[i], "-m")) distr = argv[++i];
else { else {
fprintf(stderr, "Usage: %s -[hn]\n", argv[0]); fprintf(stderr, "Usage: %s -[hnmp]\n", argv[0]);
fprintf(stderr, "\t-h\tShow this message.\n"); fprintf(stderr, " -h\t\tShow this message.\n");
fprintf(stderr, "\t-n N\tThe size of sample to generate. (default: 50000)\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; return EXIT_FAILURE;
} }
} }
@ -33,24 +38,30 @@ int main(int argc, char** argv) {
gsl_rng_env_setup(); gsl_rng_env_setup();
gsl_rng *r = gsl_rng_alloc(gsl_rng_default); gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
// prepare histogram // prepare data storage
double* sample = calloc(samples, sizeof(double)); double* sample = calloc(samples, sizeof(double));
double min = -10; double min = -10;
double max = 10; double max = 10;
/* Sample generation /* Sample generation
*
* Sample points from the Landau distribution
* using the GSL Landau generator.
*/ */
fprintf(stderr, "# Sampling\n"); fprintf(stderr, "# Sampling\n");
fprintf(stderr, "generating %ld points... ", samples); fprintf(stderr, "generating %ld points... ", samples);
double x; double x;
for(size_t i=0; i<samples; i++) { /* 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")){
x = gsl_ran_landau(r); x = gsl_ran_landau(r);
sample[i] = x; sample[i] = x;
} }
if (!strcmp(distr, "moy")){
x = gsl_rng_uniform(r);
sample[i] = moyal_qdf(x, m_params);
}
}
fprintf(stderr, "done\n"); fprintf(stderr, "done\n");
// sort the sample: needed for HSM and ks tests // sort the sample: needed for HSM and ks tests
@ -65,11 +76,12 @@ int main(int argc, char** argv) {
fprintf(stderr, "\n\n# Kolmogorov-Smirnov test\n"); fprintf(stderr, "\n\n# Kolmogorov-Smirnov test\n");
double D = 0; double D = 0;
double d; double d;
for(size_t i=0; i<samples; i++) { for(size_t i = 0; i < samples; i++) {
d = fabs(landau_cdf(sample[i], NULL) - ((double)i+1)/samples); d = fabs(landau_cdf(sample[i], NULL) - ((double)i+1)/samples);
if (d > D) if (d > D)
D = d; D = d;
} }
double beta = kolmogorov_cdf(D, samples); double beta = kolmogorov_cdf(D, samples);
// print the results // print the results