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

   1 /* histogram/stat2d.c
   2  * Copyright (C) 2002  Achim Gaedke
   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  *
  22  * File histogram/stat2d.c:
  23  * Routine to return statistical values of the content of a 2D hisogram.
  24  *
  25  * Contains the routines:
  26  * gsl_histogram2d_sum sum up all bin values
  27  * gsl_histogram2d_xmean determine mean of x values
  28  * gsl_histogram2d_ymean determine mean of y values
  29  *
  30  * Author: Achim Gaedke Achim.Gaedke@zpr.uni-koeln.de
  31  * Jan. 2002
  32  *
  33  ***************************************************************/
  34
  35 #include <config.h>
  36 #include <math.h>
  37 #include <gsl/gsl_errno.h>
  38 #include <gsl/gsl_histogram2d.h>
  39
  40 /*
  41   sum up all bins of histogram2d
  42  */
  43
  44 double
  45 gsl_histogram2d_sum (const gsl_histogram2d * h)
  46 {
  47   const size_t n = h->nx * h->ny;
  48   double sum = 0;
  49   size_t i = 0;
  50
  51   while (i < n)
  52     sum += h->bin[i++];
  53
  54   return sum;
  55 }
  56
  57 double
  58 gsl_histogram2d_xmean (const gsl_histogram2d * h)
  59 {
  60   const size_t nx = h->nx;
  61   const size_t ny = h->ny;
  62   size_t i;
  63   size_t j;
  64
  65   /* Compute the bin-weighted arithmetic mean M of a histogram using the
  66      recurrence relation
  67
  68      M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
  69      W(n) = W(n-1) + w(n)
  70
  71    */
  72
  73   long double wmean = 0;
  74   long double W = 0;
  75
  76   for (i = 0; i < nx; i++)
  77     {
  78       double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0;
  79       double wi = 0;
  80
  81       for (j = 0; j < ny; j++)
  82         {
  83           double wij = h->bin[i * ny + j];
  84           if (wij > 0)
  85             wi += wij;
  86         }
  87       if (wi > 0)
  88         {
  89           W += wi;
  90           wmean += (xi - wmean) * (wi / W);
  91         }
  92     }
  93
  94   return wmean;
  95 }
  96
  97 double
  98 gsl_histogram2d_ymean (const gsl_histogram2d * h)
  99 {
 100   const size_t nx = h->nx;
 101   const size_t ny = h->ny;
 102   size_t i;
 103   size_t j;
 104
 105   /* Compute the bin-weighted arithmetic mean M of a histogram using the
 106      recurrence relation
 107
 108      M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
 109      W(n) = W(n-1) + w(n)
 110
 111    */
 112
 113   long double wmean = 0;
 114   long double W = 0;
 115
 116   for (j = 0; j < ny; j++)
 117     {
 118       double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0;
 119       double wj = 0;
 120
 121       for (i = 0; i < nx; i++)
 122         {
 123           double wij = h->bin[i * ny + j];
 124           if (wij > 0)
 125             wj += wij;
 126         }
 127
 128       if (wj > 0)
 129         {
 130           W += wj;
 131           wmean += (yj - wmean) * (wj / W);
 132         }
 133     }
 134
 135   return wmean;
 136 }
 137
 138 double
 139 gsl_histogram2d_xsigma (const gsl_histogram2d * h)
 140 {
 141   const double xmean = gsl_histogram2d_xmean (h);
 142   const size_t nx = h->nx;
 143   const size_t ny = h->ny;
 144   size_t i;
 145   size_t j;
 146
 147   /* Compute the bin-weighted arithmetic mean M of a histogram using the
 148      recurrence relation
 149
 150      M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
 151      W(n) = W(n-1) + w(n)
 152
 153    */
 154
 155   long double wvariance = 0;
 156   long double W = 0;
 157
 158   for (i = 0; i < nx; i++)
 159     {
 160       double xi = (h->xrange[i + 1] + h->xrange[i]) / 2 - xmean;
 161       double wi = 0;
 162
 163       for (j = 0; j < ny; j++)
 164         {
 165           double wij = h->bin[i * ny + j];
 166           if (wij > 0)
 167             wi += wij;
 168         }
 169
 170       if (wi > 0)
 171         {
 172           W += wi;
 173           wvariance += ((xi * xi) - wvariance) * (wi / W);
 174         }
 175     }
 176
 177   {
 178     double xsigma = sqrt (wvariance);
 179     return xsigma;
 180   }
 181 }
 182
 183 double
 184 gsl_histogram2d_ysigma (const gsl_histogram2d * h)
 185 {
 186   const double ymean = gsl_histogram2d_ymean (h);
 187   const size_t nx = h->nx;
 188   const size_t ny = h->ny;
 189   size_t i;
 190   size_t j;
 191
 192   /* Compute the bin-weighted arithmetic mean M of a histogram using the
 193      recurrence relation
 194
 195      M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
 196      W(n) = W(n-1) + w(n)
 197
 198    */
 199
 200   long double wvariance = 0;
 201   long double W = 0;
 202
 203   for (j = 0; j < ny; j++)
 204     {
 205       double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean;
 206       double wj = 0;
 207
 208       for (i = 0; i < nx; i++)
 209         {
 210           double wij = h->bin[i * ny + j];
 211           if (wij > 0)
 212             wj += wij;
 213         }
 214       if (wj > 0)
 215         {
 216           W += wj;
 217           wvariance += ((yj * yj) - wvariance) * (wj / W);
 218         }
 219     }
 220
 221   {
 222     double ysigma = sqrt (wvariance);
 223     return ysigma;
 224   }
 225 }
 226
 227 double
 228 gsl_histogram2d_cov (const gsl_histogram2d * h)
 229 {
 230   const double xmean = gsl_histogram2d_xmean (h);
 231   const double ymean = gsl_histogram2d_ymean (h);
 232   const size_t nx = h->nx;
 233   const size_t ny = h->ny;
 234   size_t i;
 235   size_t j;
 236
 237   /* Compute the bin-weighted arithmetic mean M of a histogram using the
 238      recurrence relation
 239
 240      M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
 241      W(n) = W(n-1) + w(n)
 242
 243    */
 244
 245   long double wcovariance = 0;
 246   long double W = 0;
 247
 248   for (j = 0; j < ny; j++)
 249     {
 250       for (i = 0; i < nx; i++)
 251         {
 252           double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0 - xmean;
 253           double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean;
 254           double wij = h->bin[i * ny + j];
 255
 256           if (wij > 0)
 257             {
 258               W += wij;
 259               wcovariance += ((xi * yj) - wcovariance) * (wij / W);
 260             }
 261         }
 262     }
 263
 264   return wcovariance;
 265
 266 }