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 License along
15 * with this library; if not, write to the Free Software Foundation, Inc.,
16 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18 /***************************************************************
19 *
20 * File gsl_histogram_stat.c:
21 * Routines for statisticalcomputations on histograms.
22 * Need GSL library and header.
23 * Contains the routines:
24 * gsl_histogram_mean compute histogram mean
25 * gsl_histogram_sigma compute histogram sigma
26 *
27 * Author: S. Piccardi
28 * Jan. 2000
29 *
30 ***************************************************************/
31 #include <config.h>
32 #include <stdlib.h>
33 #include <math.h>
34 #include <gsl/gsl_errno.h>
35 #include <gsl/gsl_histogram.h>
36
37 /* FIXME: We skip negative values in the histogram h->bin[i] < 0,
38 since those correspond to negative weights (BJG) */
39
40 double
gsl_histogram_mean(const gsl_histogram * h)41 gsl_histogram_mean (const gsl_histogram * h)
42 {
43 const size_t n = h->n;
44 size_t i;
45
46 /* Compute the bin-weighted arithmetic mean M of a histogram using the
47 recurrence relation
48
49 M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
50 W(n) = W(n-1) + w(n)
51
52 */
53
54 long double wmean = 0;
55 long double W = 0;
56
57 for (i = 0; i < n; i++)
58 {
59 double xi = (h->range[i + 1] + h->range[i]) / 2;
60 double wi = h->bin[i];
61
62 if (wi > 0)
63 {
64 W += wi;
65 wmean += (xi - wmean) * (wi / W);
66 }
67 }
68
69 return wmean;
70 }
71
72 double
gsl_histogram_sigma(const gsl_histogram * h)73 gsl_histogram_sigma (const gsl_histogram * h)
74 {
75 const size_t n = h->n;
76 size_t i;
77
78 long double wvariance = 0 ;
79 long double wmean = 0;
80 long double W = 0;
81
82 /* Use a two-pass algorithm for stability. Could also use a single
83 pass formula, as given in N.J.Higham 'Accuracy and Stability of
84 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
gsl_histogram_sum(const gsl_histogram * h)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
141