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

   1 /* gsl_histogram_stat.c
   2  * Copyright (C) 2000  Simone Piccardi
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License as
   6  * published by the Free Software Foundation; either version 3 of the
   7  * License, or (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19 /***************************************************************
  20  *
  21  * File gsl_histogram_stat.c:
  22  * Routines for statisticalcomputations on histograms.
  23  * Need GSL library and header.
  24  * Contains the routines:
  25  * gsl_histogram_mean    compute histogram mean
  26  * gsl_histogram_sigma   compute histogram sigma
  27  *
  28  * Author: S. Piccardi
  29  * Jan. 2000
  30  *
  31  ***************************************************************/
  32 #include <config.h>
  33 #include <stdlib.h>
  34 #include <math.h>
  35 #include <gsl/gsl_errno.h>
  36 #include <gsl/gsl_histogram.h>
  37
  38 /* FIXME: We skip negative values in the histogram h->bin[i] < 0,
  39    since those correspond to negative weights (BJG) */
  40
  41 double
  42 gsl_histogram_mean (const gsl_histogram * h)
  43 {
  44   const size_t n = h->n;
  45   size_t i;
  46
  47   /* Compute the bin-weighted arithmetic mean M of a histogram using the
  48      recurrence relation
  49
  50      M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
  51      W(n) = W(n-1) + w(n)
  52
  53    */
  54
  55   long double wmean = 0;
  56   long double W = 0;
  57
  58   for (i = 0; i < n; i++)
  59     {
  60       double xi = (h->range[i + 1] + h->range[i]) / 2;
  61       double wi = h->bin[i];
  62
  63       if (wi > 0)
  64         {
  65           W += wi;
  66           wmean += (xi - wmean) * (wi / W);
  67         }
  68     }
  69
  70   return wmean;
  71 }
  72
  73 double
  74 gsl_histogram_sigma (const gsl_histogram * h)
  75 {
  76   const size_t n = h->n;
  77   size_t i;
  78
  79   long double wvariance = 0 ;
  80   long double wmean = 0;
  81   long double W = 0;
  82
  83   /* FIXME: should use a single pass formula here, as given in
  84      N.J.Higham 'Accuracy and Stability of Numerical Methods', p.12 */
  85
  86   /* Compute the mean */
  87
  88   for (i = 0; i < n; i++)
  89     {
  90       double xi = (h->range[i + 1] + h->range[i]) / 2;
  91       double wi = h->bin[i];
  92
  93       if (wi > 0)
  94         {
  95           W += wi;
  96           wmean += (xi - wmean) * (wi / W);
  97         }
  98     }
  99
 100   /* Compute the variance */
 101
 102   W = 0.0;
 103
 104   for (i = 0; i < n; i++)
 105     {
 106       double xi = ((h->range[i + 1]) + (h->range[i])) / 2;
 107       double wi = h->bin[i];
 108
 109       if (wi > 0) {
 110         const long double delta = (xi - wmean);
 111         W += wi ;
 112         wvariance += (delta * delta - wvariance) * (wi / W);
 113       }
 114     }
 115
 116   {
 117     double sigma = sqrt (wvariance) ;
 118     return sigma;
 119   }
 120 }
 121
 122
 123 /*
 124   sum up all bins of histogram
 125  */
 126
 127 double
 128 gsl_histogram_sum(const gsl_histogram * h)
 129 {
 130   double sum=0;
 131   size_t i=0;
 132   size_t n;
 133   n=h->n;
 134
 135   while(i < n)
 136     sum += h->bin[i++];
 137
 138   return sum;
 139 }
 140