Statistics/src/LinearFitErrorsIn2Coord.cc

0001 #include "CommonTools/Statistics/interface/LinearFitErrorsIn2Coord.h"
0002 #include "CommonTools/Statistics/interface/LinearFit.h"
0003 #include <cmath>
0004
0005 float LinearFitErrorsIn2Coord::slope(const std::vector<float>& x,
0006                                      const std::vector<float>& y,
0007                                      int ndat,
0008                                      const std::vector<float>& sigx,
0009                                      const std::vector<float>& sigy) const {
0010   // scale y and sigy, compute scaled errors
0011   float scale = sqrt(variance(x, ndat) / variance(y, ndat));
0012   std::vector<float> yScaled = y;
0013   std::vector<float> sigyScaled = sigy;
0014   std::vector<float> sig(ndat);
0015   for (int i = 0; i != ndat; i++) {
0016     yScaled[i] *= scale;
0017     sigyScaled[i] *= scale;
0018     sig[i] = sqrt(sigx[i] * sigx[i] + sigyScaled[i] * sigyScaled[i]);
0019   }
0020
0021   // usual linear fit
0022   LinearFit lf;
0023   float fs, fi, covss, covii, covsi;
0024   lf.fit(x, yScaled, ndat, sig, fs, fi, covss, covii, covsi);
0025
0026   // unscale result
0027   fs /= scale;
0028   return fs;
0029 }
0030
0031 float LinearFitErrorsIn2Coord::intercept(const std::vector<float>& x,
0032                                          const std::vector<float>& y,
0033                                          int ndat,
0034                                          const std::vector<float>& sigx,
0035                                          const std::vector<float>& sigy) const {
0036   float fs = slope(x, y, ndat, sigx, sigy);
0037   float fi = 0;
0038   float sumWi = 0;
0039   for (int i = 0; i != ndat; i++) {
0040     float wi = 1. / (sigy[i] + fs * fs * sigx[i]);
0041     fi += wi * (y[i] - fs * x[i]);
0042     sumWi += wi;
0043   }
0044
0045   return fi / sumWi;
0046 }
0047
0048 float LinearFitErrorsIn2Coord::variance(const std::vector<float>& x, int ndat) const {
0049   double m1 = 0., m2 = 0.;
0050   for (int i = 0; i != ndat; i++) {
0051     m1 += x[i];
0052     m2 += x[i] * x[i];
0053   }
0054   m1 /= ndat;
0055   m2 /= ndat;
0056
0057   return float(m2 - m1 * m1);
0058 }