htswanalysis/MACS/lib/gsl/gsl-1.11/histogram/test1d_resample.c

   1 /* histogram/test1d_resample.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23 #include <gsl/gsl_histogram.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_ieee_utils.h>
  26
  27 #include "urand.c"
  28
  29 void
  30 test1d_resample (void)
  31 {
  32   size_t i;
  33   int status = 0;
  34
  35   gsl_histogram *h;
  36
  37   gsl_ieee_env_setup ();
  38
  39   h = gsl_histogram_calloc_uniform (10, 0.0, 1.0);
  40
  41   gsl_histogram_increment (h, 0.1);
  42   gsl_histogram_increment (h, 0.2);
  43   gsl_histogram_increment (h, 0.2);
  44   gsl_histogram_increment (h, 0.3);
  45
  46   {
  47     gsl_histogram_pdf *p = gsl_histogram_pdf_alloc (10);
  48
  49     gsl_histogram *hh = gsl_histogram_calloc_uniform (100, 0.0, 1.0);
  50
  51     gsl_histogram_pdf_init (p, h);
  52
  53     for (i = 0; i < 100000; i++)
  54       {
  55         double u = urand();
  56         double x = gsl_histogram_pdf_sample (p, u);
  57         gsl_histogram_increment (hh, x);
  58       }
  59
  60     for (i = 0; i < 100; i++)
  61       {
  62         double y = gsl_histogram_get (hh, i) / 2500;
  63         double x, xmax;
  64         size_t k;
  65         double ya;
  66
  67         gsl_histogram_get_range (hh, i, &x, &xmax);
  68
  69         gsl_histogram_find (h, x, &k);
  70         ya = gsl_histogram_get (h, k);
  71
  72         if (ya == 0)
  73           {
  74             if (y != 0)
  75               {
  76                 printf ("%d: %g vs %g\n", (int) i, y, ya);
  77                 status = 1;
  78               }
  79           }
  80         else
  81           {
  82             double err = 1 / sqrt (gsl_histogram_get (hh, i));
  83             double sigma = fabs ((y - ya) / (ya * err));
  84             if (sigma > 3)
  85               {
  86                 status = 1;
  87                 printf ("%g vs %g err=%g sigma=%g\n", y, ya, err, sigma);
  88               }
  89           }
  90       }
  91
  92     gsl_histogram_pdf_free (p) ;
  93     gsl_histogram_free (hh);
  94
  95     gsl_test (status, "gsl_histogram_pdf_sample within statistical errors");
  96   }
  97
  98   gsl_histogram_free (h);
  99 }