htswanalysis/MACS/lib/gsl/gsl-1.11/statistics/wvariance_source.c

   1 /* statistics/wvariance_source.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Jim Davies, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 static double
  21 FUNCTION(compute,wvariance) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean);
  22
  23 static double
  24 FUNCTION(compute,wtss) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean);
  25
  26
  27 static double
  28 FUNCTION(compute,factor) (const BASE w[], const size_t wstride, const size_t n);
  29
  30 static double
  31 FUNCTION(compute,wvariance) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
  32 {
  33   /* takes a dataset and finds the weighted variance */
  34
  35   long double wvariance = 0 ;
  36   long double W = 0;
  37
  38   size_t i;
  39
  40   /* find the sum of the squares */
  41   for (i = 0; i < n; i++)
  42     {
  43       BASE wi = w[i * wstride];
  44
  45       if (wi > 0) {
  46         const long double delta = (data[i * stride] - wmean);
  47         W += wi ;
  48         wvariance += (delta * delta - wvariance) * (wi / W);
  49       }
  50     }
  51
  52   return wvariance ;
  53 }
  54
  55 static double
  56 FUNCTION(compute,wtss) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
  57 {
  58   /* takes a dataset and finds the weighted sum of squares about wmean*/
  59
  60   long double wtss = 0 ;
  61   size_t i;
  62
  63   /* find the sum of the squares */
  64   for (i = 0; i < n; i++)
  65     {
  66       BASE wi = w[i * wstride];
  67
  68       if (wi > 0) {
  69         const long double delta = (data[i * stride] - wmean);
  70         wtss += wi * delta * delta;
  71       }
  72     }
  73
  74   return wtss ;
  75 }
  76
  77
  78 static double
  79 FUNCTION(compute,factor) (const BASE w[], const size_t wstride, const size_t n)
  80 {
  81   /* Find the factor ``N/(N-1)'' which multiplies the raw std dev */
  82
  83   long double a = 0 ;
  84   long double b = 0;
  85   long double factor;
  86
  87   size_t i;
  88
  89   /* find the sum of the squares */
  90   for (i = 0; i < n; i++)
  91     {
  92       BASE wi = w[i * wstride];
  93
  94       if (wi > 0)
  95         {
  96           a += wi ;
  97           b += wi * wi ;
  98         }
  99     }
 100
 101   factor = (a*a) / ((a*a) - b);
 102
 103   return factor ;
 104 }
 105
 106 double
 107 FUNCTION(gsl_stats,wvariance_with_fixed_mean) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
 108 {
 109   const double wvariance = FUNCTION(compute,wvariance) (w, wstride, data, stride, n, wmean);
 110   return wvariance;
 111 }
 112
 113 double
 114 FUNCTION(gsl_stats,wsd_with_fixed_mean) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
 115 {
 116   const double wvariance = FUNCTION(compute,wvariance) (w, wstride, data, stride, n, wmean);
 117   const double wsd = sqrt (wvariance);
 118
 119   return wsd;
 120 }
 121
 122
 123 double
 124 FUNCTION(gsl_stats,wvariance_m) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
 125 {
 126   const double variance = FUNCTION(compute,wvariance) (w, wstride, data, stride, n, wmean);
 127   const double scale = FUNCTION(compute,factor)(w, wstride, n);
 128
 129   return scale * variance;
 130 }
 131
 132 double
 133 FUNCTION(gsl_stats,wsd_m) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
 134 {
 135   const double variance = FUNCTION(compute,wvariance) (w, wstride, data, stride, n, wmean);
 136   const double scale = FUNCTION(compute,factor)(w, wstride, n);
 137   const double wsd = sqrt(scale * variance) ;
 138
 139   return wsd;
 140 }
 141
 142 double
 143 FUNCTION(gsl_stats,wsd) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n)
 144 {
 145   const double wmean = FUNCTION(gsl_stats,wmean) (w, wstride, data, stride, n);
 146   return FUNCTION(gsl_stats,wsd_m) (w, wstride, data, stride, n, wmean) ;
 147 }
 148
 149 double
 150 FUNCTION(gsl_stats,wvariance) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n)
 151 {
 152   const double wmean = FUNCTION(gsl_stats,wmean) (w, wstride, data, stride, n);
 153   return FUNCTION(gsl_stats,wvariance_m)(w, wstride, data, stride, n, wmean);
 154 }
 155
 156 double
 157 FUNCTION(gsl_stats,wtss_m) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
 158 {
 159   const double wtss = FUNCTION(compute,wtss) (w, wstride, data, stride, n, wmean);
 160   return wtss;
 161 }
 162
 163 double
 164 FUNCTION(gsl_stats,wtss) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n)
 165 {
 166   const double wmean = FUNCTION(gsl_stats,wmean) (w, wstride, data, stride, n);
 167   return FUNCTION(gsl_stats,wtss_m)(w, wstride, data, stride, n, wmean);
 168 }
 169
 170