LinearFit.h

Go to the documentation of this file.

#ifndef LinearFit_H
#define LinearFit_H

#include <vector>

class LinearFit {
public:
  void fit(const std::vector<float>& x,
           const std::vector<float>& y,
           int ndat,
           const std::vector<float>& sigy,
           float& slope,
           float& intercept,
           float& covss,
           float& covii,
           float& covsi) const;
};

// template version, no std (error alrady double...)
template <typename T>
void linearFit(T const* __restrict__ x,
               T const* __restrict__ y,
               int ndat,
               T const* __restrict__ sigy2,
               T& slope,
               T& intercept,
               T& covss,
               T& covii,
               T& covsi) {
  T g1 = 0, g2 = 0;
  T s11 = 0, s12 = 0, s22 = 0;
  for (int i = 0; i != ndat; i++) {
    T sy2 = T(1) / sigy2[i];
    g1 += y[i] * sy2;
    g2 += x[i] * y[i] * sy2;
    s11 += sy2;
    s12 += x[i] * sy2;
    s22 += x[i] * x[i] * sy2;
  }

  T d = T(1) / (s11 * s22 - s12 * s12);
  intercept = (g1 * s22 - g2 * s12) * d;
  slope = (g2 * s11 - g1 * s12) * d;

  covii = s22 * d;
  covss = s11 * d;
  covsi = -s12 * d;
}

#endif