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 "gsl__config.h"
36 #include <math.h>
37 #include "gsl_errno.h"
38 #include "gsl_histogram2d.h"
39 
40 /*
41   sum up all bins of histogram2d
42  */
43 
44 double
gsl_histogram2d_sum(const gsl_histogram2d * h)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
gsl_histogram2d_xmean(const gsl_histogram2d * h)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
gsl_histogram2d_ymean(const gsl_histogram2d * h)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
gsl_histogram2d_xsigma(const gsl_histogram2d * h)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
gsl_histogram2d_ysigma(const gsl_histogram2d * h)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
gsl_histogram2d_cov(const gsl_histogram2d * h)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 }
267