Statistics/interface/LinearFit.h

0001 #ifndef LinearFit_H
0002 #define LinearFit_H
0003
0004 #include <vector>
0005
0006 /** Straight line fit for data with errors on one coordinate
0007  */
0008
0009 class LinearFit {
0010 public:
0011   /** x_i, y_i: measurements of y(x_i), sigy_i: error (sigma) on y_i
0012    *  slope, intercept: fitted parameters
0013    *  covss, covii, covsi: covariance matrix of fitted parameters,
0014    *  s denoting slope, i denoting intercept
0015    */
0016   void fit(const std::vector<float>& x,
0017            const std::vector<float>& y,
0018            int ndat,
0019            const std::vector<float>& sigy,
0020            float& slope,
0021            float& intercept,
0022            float& covss,
0023            float& covii,
0024            float& covsi) const;
0025 };
0026
0027 // template version, no std (error alrady double...)
0028 template <typename T>
0029 void linearFit(T const* __restrict__ x,
0030                T const* __restrict__ y,
0031                int ndat,
0032                T const* __restrict__ sigy2,
0033                T& slope,
0034                T& intercept,
0035                T& covss,
0036                T& covii,
0037                T& covsi) {
0038   T g1 = 0, g2 = 0;
0039   T s11 = 0, s12 = 0, s22 = 0;
0040   for (int i = 0; i != ndat; i++) {
0041     T sy2 = T(1) / sigy2[i];
0042     g1 += y[i] * sy2;
0043     g2 += x[i] * y[i] * sy2;
0044     s11 += sy2;
0045     s12 += x[i] * sy2;
0046     s22 += x[i] * x[i] * sy2;
0047   }
0048
0049   T d = T(1) / (s11 * s22 - s12 * s12);
0050   intercept = (g1 * s22 - g2 * s12) * d;
0051   slope = (g2 * s11 - g1 * s12) * d;
0052
0053   covii = s22 * d;
0054   covss = s11 * d;
0055   covsi = -s12 * d;
0056 }
0057
0058 #endif