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 "gsl__config.h"
33 #include <stdlib.h>
34 #include <math.h>
35 #include "gsl_errno.h"
36 #include "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
gsl_histogram_mean(const gsl_histogram * h)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
gsl_histogram_sigma(const gsl_histogram * h)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
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