htswanalysis/MACS/lib/gsl/gsl-1.11/fit/linear.c

   1 /* fit/linear.c
   2  *
   3  * Copyright (C) 2000, 2007 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 #include <config.h>
  21 #include <gsl/gsl_errno.h>
  22 #include <gsl/gsl_fit.h>
  23
  24 /* Fit the data (x_i, y_i) to the linear relationship
  25
  26    Y = c0 + c1 x
  27
  28    returning,
  29
  30    c0, c1  --  coefficients
  31    cov00, cov01, cov11  --  variance-covariance matrix of c0 and c1,
  32    sumsq   --   sum of squares of residuals
  33
  34    This fit can be used in the case where the errors for the data are
  35    uknown, but assumed equal for all points. The resulting
  36    variance-covariance matrix estimates the error in the coefficients
  37    from the observed variance of the points around the best fit line.
  38 */
  39
  40 int
  41 gsl_fit_linear (const double *x, const size_t xstride,
  42                 const double *y, const size_t ystride,
  43                 const size_t n,
  44                 double *c0, double *c1,
  45                 double *cov_00, double *cov_01, double *cov_11, double *sumsq)
  46 {
  47   double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0;
  48
  49   size_t i;
  50
  51   for (i = 0; i < n; i++)
  52     {
  53       m_x += (x[i * xstride] - m_x) / (i + 1.0);
  54       m_y += (y[i * ystride] - m_y) / (i + 1.0);
  55     }
  56
  57   for (i = 0; i < n; i++)
  58     {
  59       const double dx = x[i * xstride] - m_x;
  60       const double dy = y[i * ystride] - m_y;
  61
  62       m_dx2 += (dx * dx - m_dx2) / (i + 1.0);
  63       m_dxdy += (dx * dy - m_dxdy) / (i + 1.0);
  64     }
  65
  66   /* In terms of y = a + b x */
  67
  68   {
  69     double s2 = 0, d2 = 0;
  70     double b = m_dxdy / m_dx2;
  71     double a = m_y - m_x * b;
  72
  73     *c0 = a;
  74     *c1 = b;
  75
  76     /* Compute chi^2 = \sum (y_i - (a + b * x_i))^2 */
  77
  78     for (i = 0; i < n; i++)
  79       {
  80         const double dx = x[i * xstride] - m_x;
  81         const double dy = y[i * ystride] - m_y;
  82         const double d = dy - b * dx;
  83         d2 += d * d;
  84       }
  85
  86     s2 = d2 / (n - 2.0);        /* chisq per degree of freedom */
  87
  88     *cov_00 = s2 * (1.0 / n) * (1 + m_x * m_x / m_dx2);
  89     *cov_11 = s2 * 1.0 / (n * m_dx2);
  90
  91     *cov_01 = s2 * (-m_x) / (n * m_dx2);
  92
  93     *sumsq = d2;
  94   }
  95
  96   return GSL_SUCCESS;
  97 }
  98
  99
 100 /* Fit the weighted data (x_i, w_i, y_i) to the linear relationship
 101
 102    Y = c0 + c1 x
 103
 104    returning,
 105
 106    c0, c1  --  coefficients
 107    s0, s1  --  the standard deviations of c0 and c1,
 108    r       --  the correlation coefficient between c0 and c1,
 109    chisq   --  weighted sum of squares of residuals */
 110
 111 int
 112 gsl_fit_wlinear (const double *x, const size_t xstride,
 113                  const double *w, const size_t wstride,
 114                  const double *y, const size_t ystride,
 115                  const size_t n,
 116                  double *c0, double *c1,
 117                  double *cov_00, double *cov_01, double *cov_11,
 118                  double *chisq)
 119 {
 120
 121   /* compute the weighted means and weighted deviations from the means */
 122
 123   /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */
 124
 125   double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0;
 126
 127   size_t i;
 128
 129   for (i = 0; i < n; i++)
 130     {
 131       const double wi = w[i * wstride];
 132
 133       if (wi > 0)
 134         {
 135           W += wi;
 136           wm_x += (x[i * xstride] - wm_x) * (wi / W);
 137           wm_y += (y[i * ystride] - wm_y) * (wi / W);
 138         }
 139     }
 140
 141   W = 0;                        /* reset the total weight */
 142
 143   for (i = 0; i < n; i++)
 144     {
 145       const double wi = w[i * wstride];
 146
 147       if (wi > 0)
 148         {
 149           const double dx = x[i * xstride] - wm_x;
 150           const double dy = y[i * ystride] - wm_y;
 151
 152           W += wi;
 153           wm_dx2 += (dx * dx - wm_dx2) * (wi / W);
 154           wm_dxdy += (dx * dy - wm_dxdy) * (wi / W);
 155         }
 156     }
 157
 158   /* In terms of y = a + b x */
 159
 160   {
 161     double d2 = 0;
 162     double b = wm_dxdy / wm_dx2;
 163     double a = wm_y - wm_x * b;
 164
 165     *c0 = a;
 166     *c1 = b;
 167
 168     *cov_00 = (1 / W) * (1 + wm_x * wm_x / wm_dx2);
 169     *cov_11 = 1 / (W * wm_dx2);
 170
 171     *cov_01 = -wm_x / (W * wm_dx2);
 172
 173     /* Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2 */
 174
 175     for (i = 0; i < n; i++)
 176       {
 177         const double wi = w[i * wstride];
 178
 179         if (wi > 0)
 180           {
 181             const double dx = x[i * xstride] - wm_x;
 182             const double dy = y[i * ystride] - wm_y;
 183             const double d = dy - b * dx;
 184             d2 += wi * d * d;
 185           }
 186       }
 187
 188     *chisq = d2;
 189   }
 190
 191   return GSL_SUCCESS;
 192 }
 193
 194
 195
 196 int
 197 gsl_fit_linear_est (const double x,
 198                     const double c0, const double c1,
 199                     const double cov00, const double cov01, const double cov11,
 200                     double *y, double *y_err)
 201 {
 202   *y = c0 + c1 * x;
 203   *y_err = sqrt (cov00 + x * (2 * cov01 + cov11 * x));
 204   return GSL_SUCCESS;
 205 }
 206
 207
 208 int
 209 gsl_fit_mul (const double *x, const size_t xstride,
 210              const double *y, const size_t ystride,
 211              const size_t n,
 212              double *c1, double *cov_11, double *sumsq)
 213 {
 214   double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0;
 215
 216   size_t i;
 217
 218   for (i = 0; i < n; i++)
 219     {
 220       m_x += (x[i * xstride] - m_x) / (i + 1.0);
 221       m_y += (y[i * ystride] - m_y) / (i + 1.0);
 222     }
 223
 224   for (i = 0; i < n; i++)
 225     {
 226       const double dx = x[i * xstride] - m_x;
 227       const double dy = y[i * ystride] - m_y;
 228
 229       m_dx2 += (dx * dx - m_dx2) / (i + 1.0);
 230       m_dxdy += (dx * dy - m_dxdy) / (i + 1.0);
 231     }
 232
 233   /* In terms of y =  b x */
 234
 235   {
 236     double s2 = 0, d2 = 0;
 237     double b = (m_x * m_y + m_dxdy) / (m_x * m_x + m_dx2);
 238
 239     *c1 = b;
 240
 241     /* Compute chi^2 = \sum (y_i -  b * x_i)^2 */
 242
 243     for (i = 0; i < n; i++)
 244       {
 245         const double dx = x[i * xstride] - m_x;
 246         const double dy = y[i * ystride] - m_y;
 247         const double d = (m_y - b * m_x) + dy - b * dx;
 248         d2 += d * d;
 249       }
 250
 251     s2 = d2 / (n - 1.0);        /* chisq per degree of freedom */
 252
 253     *cov_11 = s2 * 1.0 / (n * (m_x * m_x + m_dx2));
 254
 255     *sumsq = d2;
 256   }
 257
 258   return GSL_SUCCESS;
 259 }
 260
 261
 262 int
 263 gsl_fit_wmul (const double *x, const size_t xstride,
 264               const double *w, const size_t wstride,
 265               const double *y, const size_t ystride,
 266               const size_t n,
 267               double *c1, double *cov_11, double *chisq)
 268 {
 269
 270   /* compute the weighted means and weighted deviations from the means */
 271
 272   /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */
 273
 274   double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0;
 275
 276   size_t i;
 277
 278   for (i = 0; i < n; i++)
 279     {
 280       const double wi = w[i * wstride];
 281
 282       if (wi > 0)
 283         {
 284           W += wi;
 285           wm_x += (x[i * xstride] - wm_x) * (wi / W);
 286           wm_y += (y[i * ystride] - wm_y) * (wi / W);
 287         }
 288     }
 289
 290   W = 0;                        /* reset the total weight */
 291
 292   for (i = 0; i < n; i++)
 293     {
 294       const double wi = w[i * wstride];
 295
 296       if (wi > 0)
 297         {
 298           const double dx = x[i * xstride] - wm_x;
 299           const double dy = y[i * ystride] - wm_y;
 300
 301           W += wi;
 302           wm_dx2 += (dx * dx - wm_dx2) * (wi / W);
 303           wm_dxdy += (dx * dy - wm_dxdy) * (wi / W);
 304         }
 305     }
 306
 307   /* In terms of y = b x */
 308
 309   {
 310     double d2 = 0;
 311     double b = (wm_x * wm_y + wm_dxdy) / (wm_x * wm_x + wm_dx2);
 312
 313     *c1 = b;
 314
 315     *cov_11 = 1 / (W * (wm_x * wm_x + wm_dx2));
 316
 317     /* Compute chi^2 = \sum w_i (y_i - b * x_i)^2 */
 318
 319     for (i = 0; i < n; i++)
 320       {
 321         const double wi = w[i * wstride];
 322
 323         if (wi > 0)
 324           {
 325             const double dx = x[i * xstride] - wm_x;
 326             const double dy = y[i * ystride] - wm_y;
 327             const double d = (wm_y - b * wm_x) + (dy - b * dx);
 328             d2 += wi * d * d;
 329           }
 330       }
 331
 332     *chisq = d2;
 333   }
 334
 335   return GSL_SUCCESS;
 336 }
 337
 338 int
 339 gsl_fit_mul_est (const double x,
 340                  const double c1, const double cov11,
 341                  double *y, double *y_err)
 342 {
 343   *y = c1 * x;
 344   *y_err = sqrt (cov11) * fabs (x);
 345   return GSL_SUCCESS;
 346 }