11
$\begingroup$

I am trying to calculate the Pearson correlation coefficient according this formula over a large dataset:

$$r = \frac{n\left(\Sigma xy\right)-\left(\Sigma x\right)\left(\Sigma y\right)}{\sqrt{\left[n\Sigma x^2 - \left(x\right)^2\right]\left[n\Sigma y^2 - \left(y\right)^2\right]}}$$

Mostly, my values are between $-1$ and $1$, but sometimes I get weird numbers like $1.0000000002$ or $-3$.

And so on. Is it possible to have weird data that would result in this, or does this mean that I have an error in calculation?

For example, I notice that sometimes my summation of $x$ is $1$, and hence summation of $x^2$ would be $1$. This results in a value like $1.00000002$. Other times, I will have the summation of $xy$ as $0$, and then I will have the resulting calculation be $-3$. Is this statistically possible, or is there an error in my calculations?

$\endgroup$
7
  • 2
    $\begingroup$ What language or environment are you using? $\endgroup$ Commented Sep 14, 2016 at 18:22
  • 4
    $\begingroup$ It would be helpful to know a bit about the size of the numbers you were dealing with, how many of them there were, and the level of accuracy of your intermediate calculations e.g. $\sum xy$... there's clearly a numerical stability issue here that might be worth exploring. $\endgroup$
    – Silverfish
    Commented Sep 14, 2016 at 18:26
  • 1
    $\begingroup$ I second @Silverfish. Perhaps you can post an example that we can evaluate. N.b.1) you can access Chrome's JavaScript console with Ctrl+Shift+J. N.b.2) All JS numbers are 64 bit double w3schools.com/js/js_numbers.asp $\endgroup$ Commented Sep 14, 2016 at 18:29
  • 2
    $\begingroup$ Tongue in cheek answer: It is not possible to have $R>1$ or $R<-1$ mathematically (i.e. for $R\in\mathbb{R}$), but it is possible to have NOT((R>=-1)&(R<=1)) be True in IEEE arithmetic, if $x$ and/or $y$ are constant (as 0/0 equals NaN, which fails all comparisons). $\endgroup$
    – GeoMatt22
    Commented Sep 15, 2016 at 1:53
  • 1
    $\begingroup$ "For my Y dataset, the numbers are 0 < Y < 1 and generally anywhere from e-5 to e-350. For my X dataset, the numbers are anywhere from 0 to e7" O.k sports fans, such a wide range of orders of magnitudes of the numbers is not a recipe for success, particularly for numerically non-robust algorithms, but maybe not so great with them. $\endgroup$ Commented Sep 15, 2016 at 22:59

3 Answers 3

18
$\begingroup$

The formulas you're using have long been known to be numerically unstable. If the squared means are large compared to the variances and/or products-of-means are large compared to the covariances, then the difference in the numerator and in the bracketed terms in the denominator can have problems with catastrophic cancellation.

This can sometimes lead to calculated variances or covariances that don't even retain a single digit of precision (i.e. that are worse than useless).

Don't use these formulas. They made some sense when people calculated by hand, where you could see, and deal with such loss of precision when it happened – e.g., use of these formulas was normally preceded by eliminating the common digits, so numbers like this:

 8901234.567...
 8901234.575...
 8901234.412...

would first have 8901234 subtracted (at least) – which would save a lot of time in the working as well as avoid the cancellation issue. Means (and similar quantities) would then be adjusted back at the end, while variances and covariances could be used as-is.

Similar ideas (and other ideas) can be used with computers, but really you need to use them all the time, rather than trying to guess when you might need them.

Efficient ways to deal with this issue have been known for over half a century – e.g., see Welford's 1962 paper [1] (where he gives one-pass variance and covariance algorithms – stable two-pass algorithms were well know already). Chan et al [2] (1983) compare a number of variance algorithms and offer a way to decide when to use which (though in most implementations generally people use only one algorithm).

See Wikipedia's discussion on this issue in relation to variance and its discussion on variance algorithms.

Similar comments apply to covariance.

[1] B. P. Welford (1962),
"Note on a Method for Calculating Corrected Sums of Squares and Products",
Technometrics Vol. 4 , Iss. 3, 419-420
(link)

Also see https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm

[2] T.F. Chan, G.H. Golub and R.J. LeVeque (1983)
"Algorithms for Computing the Sample Variance: Analysis and Recommendations",
The American Statistician, Vol. 37, No. 3 (Aug.1983), pp. 242-247
Tech report version

$\endgroup$
5
  • $\begingroup$ I think for the numerical issues, the simple two-pass algorithm is usually sufficiently reliable: one pass for means, second pass for (co-)variances. My impression is the fancier variants (e.g. compensated sum) are rarely required for reliability if using double precision. (For efficiency, the online and/or parallel versions can be useful though.) I am definitely badly behaved sometimes, and use the unstable form though! (for things like moving-average filters over images) $\endgroup$
    – GeoMatt22
    Commented Sep 15, 2016 at 2:18
  • 1
    $\begingroup$ (+1) @Tim gives an R implementation of Welford's algorithm at stats.stackexchange.com/a/235151/919. $\endgroup$
    – whuber
    Commented Sep 15, 2016 at 15:16
  • 1
    $\begingroup$ Glen_b, I have removed the CiteSeer link for of late they aren't working. If you find any better link than what I provided, please do edit that. CiteSeer links aren't working anymore. $\endgroup$ Commented Jan 6, 2023 at 4:51
  • 1
    $\begingroup$ Thanks, I'll take a look. $\endgroup$
    – Glen_b
    Commented Jan 6, 2023 at 6:41
  • $\begingroup$ @Glen_b, just for the information, the problem with CiteSeerX has been addressed at Maths recently. I fixed those which I could find either by archived ones or other replacements. $\endgroup$ Commented Jan 6, 2023 at 7:02
8
$\begingroup$

The Pearson correlation coefficient is indeed between $-1$ and $+1$ (inclusive). This follows from the Cauchy-Schwarz inequality.

Getting a correlation coefficient of $1.0000000002$ is possibly (but unlikely) due to numerical error, while -3 almost certainly indicates an error in implementation (or a platform unsuitable for numerics! :).

$\endgroup$
1
  • 2
    $\begingroup$ Could you possibly consider extending your answer (i.e. show how it follows from C-S ineq. and saying in few words about numerical errors)..? $\endgroup$
    – Tim
    Commented Sep 16, 2016 at 7:45
1
$\begingroup$

The likely reason is the round-off error, that is the error due to the fact that a computer performs calculations to a finite precision - keeping only certain number of significant digits of a number (e.g., only 8 digits or only 16 digits - see , e.g., here.)

It can be swept under the carpet as a numerical error, since one encounters it when performing numerical calculations. However, it is different in nature: "numerical errors" arise from an approximate nature of numerical calculations, whereas round-off error is due to the finite number of digits stored in the computer memory. Numerical errors are usually defined by the number of digits after the floating point, whereas the round-off error concerns the total number of digits. Thus, if using, e.g., decimal32 format, the computer would be able to distinguish between $0.001$ and $0.002$, but not between $7000000.001$ and $7000000.002$.

In practice the rounding-off error often manifests itself as non-zero last digit(s) in a number after performing a sequence of algebraic operations, e.g., like getting $1.00000002$ instead of $1.00000000$. One thus should treat the results as correct only up to the precision of the number format used.

Some software provide common comparison operators robust to small error. E.g., numpy.allclose(a,b) allows comparing arrays up to specified tolerance, although they might not be equal when using usual comparison, like (a == b).all()

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.