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

enter image description here

Mostly, my values are between -1 and 1, but sometimes I get weird numbers like:


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?

  • 2
    $\begingroup$ What language or environment are you using? $\endgroup$ – P.Windridge Sep 14 '16 at 18:22
  • 3
    $\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 Sep 14 '16 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$ – P.Windridge Sep 14 '16 at 18:29
  • 1
    $\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 Sep 15 '16 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$ – Mark L. Stone Sep 15 '16 at 22:59

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:


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
(citeseer link)

[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

  • $\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 Sep 15 '16 at 2:18
  • 1
    $\begingroup$ (+1) @Tim gives an R implementation of Welford's algorithm at stats.stackexchange.com/a/235151/919. $\endgroup$ – whuber Sep 15 '16 at 15:16

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! :).

  • $\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 Sep 16 '16 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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