C++: Correlation / LeastSquaresCoefs

Question

// X.size() == Y.size() > 0.  returns correlation between X and Y
double Correlation(const vector<double>& X, const vector<double>& Y)
{
    size_t N = X.size();
    double EX(0), EY(0), EXY(0), EX2(0), EY2(0);

    for (size_t i = 0; i < N; i++)
    {
        EX += X[i];
        EY += Y[i];
        EXY += X[i]*Y[i];
        EX2 += X[i]*X[i];
        EY2 += Y[i]*Y[i];
    }

    return (N*EXY - EX*EY) / sqrt((N*EX2 - EX*EX) * (N*EY2 - EY*EY));
}

// X.size() == Y.size() > 0.  returns {a,b}, where y = a*x + b is
// line of best fit with least mean squared error.
pair<double, double> LeastSquaresCoefs(const vector<double>& X, const vector<double>& Y)
{
    size_t N = X.size();
    double EX(0), EY(0), EXY(0), EX2(0), EY2(0);

    for (size_t i = 0; i < N; i++)
    {
        EX += X[i];
        EY += Y[i];
        EXY += X[i]*Y[i];
        EX2 += X[i]*X[i];
        EY2 += Y[i]*Y[i];
    }

    double b = (EX2*EY - EX*EXY) / (N * EX2 - EX*EX);
    double a = (EXY - b * EX) / EX2;

    return {a, b};
}

Do these functions work? Are they correct?

See also: http://math.stackexchange.com/questions/120941/simple-least-squares-regression

Answer 1

First, your code for Pearson correlation coeff. seems to be incorrect:

return (EXY - N*EX*EY) / sqrt((EX2 - N*EX*EX) * (EY2 - N*EY*EY))

is better.

Second, if you have big data set, you may meet "precision lost" problem, so it is better to use original formula, where you substract mean value from every set value

for (size_t i = 0; i < n; i++) { Find the means.
 ex += x[i];
 ey += y[i];
}
ex /= n;
ey /= n;
for (size_t i = 0; i < n; i++) { Compute the correlation coefficient.
  xt = x[i] - ex;
  yt = y[i] - ey;
  sxx += xt * xt;
  syy += yt * yt;
  sxy += xt * yt;
}
return sxy/(sqrt(sxx*syy)+TINY_VALUE);

where TINY_VALUE could be something like 1e-20 and is used to "compensate" perfect correlation case (and avoid special verification).

Same ideas relate to least squares implementation.

When it comes to numerical methods, algorithms should take care about precision and robustness of calculation routines.

You can find more at http://mathworld.wolfram.com/CorrelationCoefficient.html