# Pearson correlation of data sets with possibly zero standard deviation?

I am having a problem computing the pearson correlation coefficient of data sets with possibly zero standard deviation (i.e. all data has the same value).

Suppose that I have the following two data sets:

float x[] = {2, 2, 2, 3, 2};
float y[] = {2, 2, 2, 2, 2};


The correlation coefficient "r", would be computed using the following equation:

float r = covariance(x, y) / (std_dev(x) * std_dev(y));


However, because all data in data set "y" has the same value, the standard deviation std_dev(y) would be zero and "r" would be undefined.

Is there any solution for this problem? Or should I use other methods to measure data relationship in this case?

-
There is no "data relationship" in this example because y does not vary. Assigning any numerical value to r would be a mistake. –  whuber Apr 2 '11 at 4:42
@whuber - it is true that the $r$ is undefined, but not necessarily that the "true" unknown correlation $\rho$ cannot be estimated. Just have to use something different to estimate it. –  probabilityislogic Apr 2 '11 at 7:26
@probability You presuppose this is a problem of estimation and not simply one of characterization. But accepting that, what estimator would you propose in the example? No answer can be universally correct because it depends on how the estimator will be used (a loss function, in effect). In many applications, such as PCA, it seems likely that using any procedure that imputes a value to $\rho$ may be worse than other procedures that recognize $\rho$ cannot be identified. –  whuber Apr 2 '11 at 15:57
@whuber - estimate is a bad choice of words for me (you may have noticed I'm not the best wordsmith), what I meant was that although $\rho$ may not be uniquely identified, this does not mean that the data are useless in telling us about $\rho$. My answer gives an (ugly) demonstration of this from an algebraic point of view. –  probabilityislogic Apr 2 '11 at 16:59
@Probability It seems your analysis is contradictory: if indeed y is modeled with a normal distribution, then a sample of five 2's shows this model is inappropriate. Ultimately, you don't get something for nothing: your results depend strongly on the assumptions made about the priors. The original problems in identifying $\rho$ are still there but have been hidden by all these additional assumptions. That seems IMHO just to obscure the issues rather than clarify them. –  whuber Apr 2 '11 at 17:30

## migrated from stackoverflow.comApr 2 '11 at 0:11

This question came from our site for professional and enthusiast programmers.

The "sampling theory" people will tell you that no such estimate exists. But you can get one, you just need to be reasonable about your prior information, and do a lot harder mathematical work.

If you specified a Bayesian method of estimation, and the posterior is the same as the prior, then you can say the data say nothing about the parameter. Because things may get "singular" on us, then we cannot use infinite parameter spaces. I am assuming that because you use Pearson correlation, you have a bivariate normal likelihood:

$$p(D|\mu_x,\mu_y,\sigma_x,\sigma_y,\rho)=\left(\sigma_x\sigma_y\sqrt{2\pi(1-\rho^2)}\right)^{-N}exp\left(-\frac{\sum_{i}Q_i}{2(1-\rho^2)}\right)$$ where $$Q_i=\frac{(x_i-\mu_x)^2}{\sigma_x^2}+\frac{(y_i-\mu_y)^2}{\sigma_y^2}-2\rho\frac{(x_i-\mu_x)(y_i-\mu_y)}{\sigma_x\sigma_y}$$

Now to indicate that one data set may be the same value, write $y_i=y$, and then we get:

$$\sum_{i}Q_i=N\left[\frac{(y-\mu_y)^2}{\sigma_y^2}+\frac{s_x^2 + (\overline{x}-\mu_x)^2}{\sigma_x^2}-2\rho\frac{(\overline{x}-\mu_x)(y-\mu_y)}{\sigma_x\sigma_y}\right]$$ where $$s_x^2=\frac{1}{N}\sum_{i}(x_i-\overline{x})^2$$

And so your likelihood depends on four numbers, $s_x^2,y,\overline{x},N$. So you want an estimate of $\rho$, so you need to multiply by a prior, and integrate out the nuisance parameters $\mu_x,\mu_y,\sigma_x,\sigma_y$. Now to prepare for integration, we "complete the square" $$\frac{\sum_{i}Q_i}{1-\rho^2}=N\left[\frac{\left(\mu_y-\left[y-(\overline{x}-\mu_x)\frac{\rho\sigma_y}{\sigma_x}\right]\right)^2}{\sigma_y^2(1-\rho^{2})}+\frac{s_x^2}{\sigma_{x}^{2}(1-\rho^{2})} + \frac{(\overline{x}-\mu_x)^2}{\sigma_x^2}\right]$$

Now we should err on the side of caution and ensure a properly normalised probability. That way we can't get into trouble. One such option is to use a weakly informative prior, which just places restriction on the range of each. So we have $L_{\mu}<\mu_x,\mu_y<U_{\mu}$ for the means with flat prior and $L_{\sigma}<\sigma_x,\sigma_y<U_{\sigma}$ for the standard deviations with jeffreys prior. These limits are easy to set with a bit of "common sense" thinking about the problem. I will take an unspecified prior for $\rho$, and so we get (uniform should work ok, if not truncate the singularity at $\pm 1$):

$$p(\rho,\mu_x,\mu_y,\sigma_x,\sigma_y)=\frac{p(\rho)}{A\sigma_x\sigma_y}$$

Where $A=2(U_{\mu}-L_{\mu})^{2}[log(U_{\sigma})-log(L_{\sigma})]^{2}$. This gives a posterior of:

$$p(\rho|D)=\int p(\rho,\mu_x,\mu_y,\sigma_x,\sigma_y)p(D|\mu_x,\mu_y,\sigma_x,\sigma_y,\rho)d\mu_y d\mu_x d\sigma_x d\sigma_y$$

$$=\frac{p(\rho)}{A[2\pi(1-\rho^2)]^{\frac{N}{2}}}\int_{L_{\sigma}}^{U_{\sigma}}\int_{L_{\sigma}}^{U_{\sigma}}\left(\sigma_x\sigma_y\right)^{-N-1}exp\left(-\frac{N s_x^2}{2\sigma_{x}^{2}(1-\rho^{2})}\right) \times$$ $$\int_{L_{\mu}}^{U_{\mu}}exp\left(-\frac{N(\overline{x}-\mu_x)^2}{2\sigma_x^2}\right)\int_{L_{\mu}}^{U_{\mu}}exp\left(-\frac{N\left(\mu_y-\left[y-(\overline{x}-\mu_x)\frac{\rho\sigma_y}{\sigma_x}\right]\right)^2}{2\sigma_y^2(1-\rho^{2})}\right)d\mu_y d\mu_x d\sigma_x d\sigma_y$$

Now the first integration over $\mu_y$ can be done by making a change of variables $z=\sqrt{N}\frac{\mu_y-\left[y-(\overline{x}-\mu_x)\frac{\rho\sigma_y}{\sigma_x}\right]}{\sigma_y\sqrt{1-\rho^{2}}}\implies dz=\frac{\sqrt{N}}{\sigma_y\sqrt{1-\rho^{2}}}d\mu_y$ and the first integral over $\mu_y$ becomes:

$$\frac{\sigma_y\sqrt{2\pi(1-\rho^{2})}}{\sqrt{N}}\left[\Phi\left( \frac{U_{\mu}-\left[y-(\overline{x}-\mu_x)\frac{\rho\sigma_y}{\sigma_x}\right]}{\frac{\sigma_y}{\sqrt{N}}\sqrt{1-\rho^{2}}} \right)-\Phi\left( \frac{L_{\mu}-\left[y-(\overline{x}-\mu_x)\frac{\rho\sigma_y}{\sigma_x}\right]}{\frac{\sigma_y}{\sqrt{N}}\sqrt{1-\rho^{2}}} \right)\right]$$

And you can see from here, no analytic solutions are possible. However, it is also worthwhile to note that the value $\rho$ has not dropped out of the equations. This means that the data and prior information still have something to say about the true correlation. If the data said nothing about the correlation, then we would be simply left with $p(\rho)$ as the only function of $\rho$ in these equations.

It also shows how that passing to the limit of infinite bounds for $\mu_y$ "throws away" some of the information about $\rho$, which is contained in the complicated looking normal CDF function $\Phi(.)$. Now if you have a lot of data, then passing to the limit is fine, you don't loose much, but if you have very scarce information, such as in your case - it is important keep every scrap you have. It means ugly maths, but this example is not too hard to do numerically. So we can evaluate the integrated likelihood for $\rho$ at values of say $-0.99,-0.98,\dots,0.98,0.99$ fairly easily. Just replace the integrals by summations over a small enough intervals - so you have a triple summation

-
@probabilityislogic: Wow. Simply wow. After seen some of your answers I really wonder: what should a doofus like me do to reach such a flexible bayesian state of mind ? –  steffen Apr 7 '11 at 15:02
@steffen - lol. Its not that difficult, you just need to practice. And always always always remember that the product and sum rules of probability are the only rules you will ever need. They will extract whatever information is there - whether you see it or not. So you apply product and sum rules, then just do the maths. That is all I have done here. –  probabilityislogic Apr 7 '11 at 15:07
@steffen - and the other rule - more a mathematical one than stats one - don't pass to an infinite limit too early in your calculations, your results may become arbitrary, or little details may get thrown out. Measurement error models are a perfect example of this (as is this question). –  probabilityislogic Apr 7 '11 at 15:11
@probabilityislogic: Thank you, I'll keep this in mind... as soon as I am done working through my "Bayesian Analysis"-copy ;). –  steffen Apr 7 '11 at 15:33
@probabilityislogic: If you could humor a nonmathematical statistician/researcher...would it be possible to summarize or translate your answer to a group of dentists or high school principals or introductory statistics students? –  rolando2 Apr 30 '11 at 16:53

I agree with sesqu that the correlation is undefined in this case. Depending on your type of application you could e.g. calculate the Gower Similarity between both vectors, which is: $gower(v1,v2)=\frac{\sum_{i=1}^{n}\delta(v1_i,v2_i)}{n}$ where $\delta$ represents the kronecker-delta, applied as function on $v1,v2$.

So for instance if all values are equal, gower(.,.)=1. If on the other hand they differ only in one dimension, gower(.,.)=0.9. If they differ in every dimension, gower(.,.)=0 and so on.

Of course this is no measure for correlation, but it allows you to calculate how close the vector with s>0 is to the one with s=0. Of course you can apply other metrics,too, if they serve your purpose better.

-
+1 That's a creative idea. It sounds like the "Gower Similarity" is a scaled Hamming distance. –  whuber Apr 2 '11 at 18:04
@whuber: Indeed it is ! –  steffen Apr 3 '11 at 8:30