|
Packit |
67cb25 |
/* histogram/stat2d.c
|
|
Packit |
67cb25 |
* Copyright (C) 2002 Achim Gaedke
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This library is free software; you can redistribute it and/or
|
|
Packit |
67cb25 |
* modify it under the terms of the GNU General Public License as
|
|
Packit |
67cb25 |
* published by the Free Software Foundation; either version 3 of the
|
|
Packit |
67cb25 |
* License, or (at your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful,
|
|
Packit |
67cb25 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License along
|
|
Packit |
67cb25 |
* with this library; if not, write to the Free Software Foundation, Inc.,
|
|
Packit |
67cb25 |
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/***************************************************************
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* File histogram/stat2d.c:
|
|
Packit |
67cb25 |
* Routine to return statistical values of the content of a 2D hisogram.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Contains the routines:
|
|
Packit |
67cb25 |
* gsl_histogram2d_sum sum up all bin values
|
|
Packit |
67cb25 |
* gsl_histogram2d_xmean determine mean of x values
|
|
Packit |
67cb25 |
* gsl_histogram2d_ymean determine mean of y values
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Author: Achim Gaedke Achim.Gaedke@zpr.uni-koeln.de
|
|
Packit |
67cb25 |
* Jan. 2002
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
***************************************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_histogram2d.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
sum up all bins of histogram2d
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram2d_sum (const gsl_histogram2d * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = h->nx * h->ny;
|
|
Packit |
67cb25 |
double sum = 0;
|
|
Packit |
67cb25 |
size_t i = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (i < n)
|
|
Packit |
67cb25 |
sum += h->bin[i++];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return sum;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram2d_xmean (const gsl_histogram2d * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t nx = h->nx;
|
|
Packit |
67cb25 |
const size_t ny = h->ny;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the bin-weighted arithmetic mean M of a histogram using the
|
|
Packit |
67cb25 |
recurrence relation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
|
|
Packit |
67cb25 |
W(n) = W(n-1) + w(n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long double wmean = 0;
|
|
Packit |
67cb25 |
long double W = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < nx; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0;
|
|
Packit |
67cb25 |
double wi = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < ny; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wij = h->bin[i * ny + j];
|
|
Packit |
67cb25 |
if (wij > 0)
|
|
Packit |
67cb25 |
wi += wij;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (wi > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
W += wi;
|
|
Packit |
67cb25 |
wmean += (xi - wmean) * (wi / W);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return wmean;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram2d_ymean (const gsl_histogram2d * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t nx = h->nx;
|
|
Packit |
67cb25 |
const size_t ny = h->ny;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the bin-weighted arithmetic mean M of a histogram using the
|
|
Packit |
67cb25 |
recurrence relation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
|
|
Packit |
67cb25 |
W(n) = W(n-1) + w(n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long double wmean = 0;
|
|
Packit |
67cb25 |
long double W = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < ny; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0;
|
|
Packit |
67cb25 |
double wj = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < nx; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wij = h->bin[i * ny + j];
|
|
Packit |
67cb25 |
if (wij > 0)
|
|
Packit |
67cb25 |
wj += wij;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (wj > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
W += wj;
|
|
Packit |
67cb25 |
wmean += (yj - wmean) * (wj / W);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return wmean;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram2d_xsigma (const gsl_histogram2d * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double xmean = gsl_histogram2d_xmean (h);
|
|
Packit |
67cb25 |
const size_t nx = h->nx;
|
|
Packit |
67cb25 |
const size_t ny = h->ny;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the bin-weighted arithmetic mean M of a histogram using the
|
|
Packit |
67cb25 |
recurrence relation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
|
|
Packit |
67cb25 |
W(n) = W(n-1) + w(n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long double wvariance = 0;
|
|
Packit |
67cb25 |
long double W = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < nx; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = (h->xrange[i + 1] + h->xrange[i]) / 2 - xmean;
|
|
Packit |
67cb25 |
double wi = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < ny; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wij = h->bin[i * ny + j];
|
|
Packit |
67cb25 |
if (wij > 0)
|
|
Packit |
67cb25 |
wi += wij;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (wi > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
W += wi;
|
|
Packit |
67cb25 |
wvariance += ((xi * xi) - wvariance) * (wi / W);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xsigma = sqrt (wvariance);
|
|
Packit |
67cb25 |
return xsigma;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram2d_ysigma (const gsl_histogram2d * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ymean = gsl_histogram2d_ymean (h);
|
|
Packit |
67cb25 |
const size_t nx = h->nx;
|
|
Packit |
67cb25 |
const size_t ny = h->ny;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the bin-weighted arithmetic mean M of a histogram using the
|
|
Packit |
67cb25 |
recurrence relation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
|
|
Packit |
67cb25 |
W(n) = W(n-1) + w(n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long double wvariance = 0;
|
|
Packit |
67cb25 |
long double W = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < ny; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean;
|
|
Packit |
67cb25 |
double wj = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < nx; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wij = h->bin[i * ny + j];
|
|
Packit |
67cb25 |
if (wij > 0)
|
|
Packit |
67cb25 |
wj += wij;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (wj > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
W += wj;
|
|
Packit |
67cb25 |
wvariance += ((yj * yj) - wvariance) * (wj / W);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ysigma = sqrt (wvariance);
|
|
Packit |
67cb25 |
return ysigma;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram2d_cov (const gsl_histogram2d * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double xmean = gsl_histogram2d_xmean (h);
|
|
Packit |
67cb25 |
const double ymean = gsl_histogram2d_ymean (h);
|
|
Packit |
67cb25 |
const size_t nx = h->nx;
|
|
Packit |
67cb25 |
const size_t ny = h->ny;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the bin-weighted arithmetic mean M of a histogram using the
|
|
Packit |
67cb25 |
recurrence relation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
|
|
Packit |
67cb25 |
W(n) = W(n-1) + w(n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long double wcovariance = 0;
|
|
Packit |
67cb25 |
long double W = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < ny; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (i = 0; i < nx; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0 - xmean;
|
|
Packit |
67cb25 |
double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean;
|
|
Packit |
67cb25 |
double wij = h->bin[i * ny + j];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (wij > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
W += wij;
|
|
Packit |
67cb25 |
wcovariance += ((xi * yj) - wcovariance) * (wij / W);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return wcovariance;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|