C++: Correlation / LeastSquaresCoefs [closed]

// 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?

-

closed as off-topic by Malachi, Simon André Forsberg, Jamal♦Nov 21 '13 at 19:26

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "Your question must contain working code for us to review it here. For questions regarding specific problems encountered while coding, try Stack Overflow. After getting your code to work, you may edit this question seeking a review of your working code." – Malachi, Simon André Forsberg, Jamal
If this question can be reworded to fit the rules in the help center, please edit the question.

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 coeﬃcient.
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

-