|
Packit |
67cb25 |
/* gsl_histogram_stat.c
|
|
Packit |
67cb25 |
* Copyright (C) 2000 Simone Piccardi
|
|
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 |
* File gsl_histogram_stat.c:
|
|
Packit |
67cb25 |
* Routines for statisticalcomputations on histograms.
|
|
Packit |
67cb25 |
* Need GSL library and header.
|
|
Packit |
67cb25 |
* Contains the routines:
|
|
Packit |
67cb25 |
* gsl_histogram_mean compute histogram mean
|
|
Packit |
67cb25 |
* gsl_histogram_sigma compute histogram sigma
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Author: S. Piccardi
|
|
Packit |
67cb25 |
* Jan. 2000
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
***************************************************************/
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_histogram.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* FIXME: We skip negative values in the histogram h->bin[i] < 0,
|
|
Packit |
67cb25 |
since those correspond to negative weights (BJG) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram_mean (const gsl_histogram * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = h->n;
|
|
Packit |
67cb25 |
size_t i;
|
|
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 < n; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = (h->range[i + 1] + h->range[i]) / 2;
|
|
Packit |
67cb25 |
double wi = h->bin[i];
|
|
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_histogram_sigma (const gsl_histogram * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = h->n;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long double wvariance = 0 ;
|
|
Packit |
67cb25 |
long double wmean = 0;
|
|
Packit |
67cb25 |
long double W = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Use a two-pass algorithm for stability. Could also use a single
|
|
Packit |
67cb25 |
pass formula, as given in N.J.Higham 'Accuracy and Stability of
|
|
Packit |
67cb25 |
Numerical Methods', p.12 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the mean */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = (h->range[i + 1] + h->range[i]) / 2;
|
|
Packit |
67cb25 |
double wi = h->bin[i];
|
|
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 |
/* Compute the variance */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
W = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = ((h->range[i + 1]) + (h->range[i])) / 2;
|
|
Packit |
67cb25 |
double wi = h->bin[i];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (wi > 0) {
|
|
Packit |
67cb25 |
const long double delta = (xi - wmean);
|
|
Packit |
67cb25 |
W += wi ;
|
|
Packit |
67cb25 |
wvariance += (delta * delta - wvariance) * (wi / W);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sigma = sqrt (wvariance) ;
|
|
Packit |
67cb25 |
return sigma;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
sum up all bins of histogram
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_histogram_sum(const gsl_histogram * h)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum=0;
|
|
Packit |
67cb25 |
size_t i=0;
|
|
Packit |
67cb25 |
size_t n;
|
|
Packit |
67cb25 |
n=h->n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while(i < n)
|
|
Packit |
67cb25 |
sum += h->bin[i++];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return sum;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|