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)45gsl_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)58gsl_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)98gsl_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)139gsl_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)184gsl_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)228gsl_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