ABACUS/src/FITTING/lin_reg.cc

/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  lin_reg.cc

Purpose:  Linear regression

***********************************************************/

#include "ABACUS.h"

using namespace std;

namespace ABACUS {

  void lin_reg (Vect_DP x, Vect_DP y, Vect_DP sigma, DP& a, DP& b, DP& chisq)
  {
    // Performs simple linear regression on data

    int Npts = x.size();

    DP S = 0.0;
    DP Sx = 0.0;
    DP Sy = 0.0;
    DP Stt = 0.0;
    Vect_DP t(0.0, Npts);

    for (int i = 0; i < Npts; ++i) {
      S += 1.0/(sigma[i] * sigma[i]);
      Sx += x[i]/(sigma[i] * sigma[i]);
      Sy += y[i]/(sigma[i] * sigma[i]);
    }

    for (int i = 0; i < Npts; ++i) {
      t[i] = (x[i] - Sx/S)/sigma[i];
      Stt += t[i] * t[i];
    }

    a = 0.0;
    b = 0.0;
    for (int i = 0; i < Npts; ++i) {
      b += t[i] * y[i]/sigma[i];
    }
    b /= Stt;
    a = (Sy - Sx * b)/S;

    chisq = 0.0;
    for (int i = 0; i < Npts; ++i) {
      chisq += (y[i] - a - b * x[i]) * (y[i] - a - b * x[i]) / (sigma[i] * sigma[i]);
    }

    return;
  }

  void lin_reg (Vect_DP x, Vect_DP y, DP& a, DP& b, DP& chisq)
  {
    // Assumes all sigma == 1

    Vect_DP sigma(1.0, x.size());

    lin_reg (x, y, sigma, a, b, chisq);

    return;
  }

}