Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I'd like to explain an apparent discrepancy between:

(1) The sample correlation distribution function between sample vectors for a bivariate, correlated random variable (correlation coefficient = $\rho$) and

(2) The sample correlation distribution function for two normally distributed random vectors that each contain a given signal with additive Gaussian noise. In this latter case, the correlation coefficient between these vectors $\rho$.

While (1) gives the standard Pearson Correlation Coefficient, (2) is more consistent with time series analysis and match filtering/correlation detection. I'd like to use a distribution function for (2) to identify the detection performance of such a correlation detector. Please read on below.

First sample correlation computation:

The "correct" way to view the correlation distribution function is as follows:

Suppose pairs $(x_{1}, s_{1})$, $(x_{2}, s_{2})$, $(x_{3}, s_{3})$, $\cdots$, $(x_{n}, s_{n})$ are a sample of size $n$ drawn from the distribution of the random vector $(S, \, X)^{\prime}$, which has a bivariate normal distribution with mean vector $(\mu_{S}, \, \mu_{X})^{\prime}$ and covariance:

\begin{equation} \Sigma = \begin{bmatrix} \sigma_{S}^{2} & \sigma_{S\,X} \\[0.3em] \sigma_{S\,X} & \sigma_{X}^{2} \end{bmatrix}, \end{equation}

where the population correlation coefficient is:

\begin{equation} \rho = \cfrac{\sigma_{S\,X}}{\sigma_{X}\, \sigma_{S}} \end{equation}

The estimator for $\rho$, call it $\hat{\rho}$, is:

\begin{equation} \hat{\rho}_{1} = \frac{ \boldsymbol{x}^{\text{T}} \boldsymbol{s} } { \left|\left|\, \boldsymbol{x} \,\right|\right| \,\left|\left| \,\boldsymbol{s} \,\right|\right| } \end{equation}

where $\boldsymbol{x}$ and $\boldsymbol{s}$ are concatenated $n$-dim vectors formed from from the samples $(x_{k}, s_{k})$ above ($k$ $=$ $1,2,\cdots,n$). Coefficient $\hat{\rho}$ has a Pearson's correlation coefficient distribution function, that I will call $p_{R}\left( r \,; \rho, n, \,\mathcal{H}_{1}\right)$. The expression $\mathcal{H}_{1}$ is borrowed from hypothesis decision theory. It indicates that the true correlation is nonzero, e.g., $\rho$ $\ne$ $0$.

Second sample correlation computation:

OK: Now take two vectors, $\boldsymbol{x}$ and $\boldsymbol{s}$. This time, assume:

$\boldsymbol{x}$ $\sim$ $\mathcal{N}(\boldsymbol{\mu}_{X}, \sigma_{X}^{2}\mathbf{I})$


$\boldsymbol{s}$ $\sim$ $\mathcal{N}(\boldsymbol{\mu}_{S}, \sigma_{S}^{2}\mathbf{I})$

Again, compute the sample correlation distribution function:

\begin{equation} \hat{\rho}_{2} = \frac{ \boldsymbol{x}^{\text{T}} \boldsymbol{s} } { \left|\left|\, \boldsymbol{x} \,\right|\right| \,\left|\left| \,\boldsymbol{s} \,\right|\right| } \end{equation}

Suppose now that $\rho$ is the same in both cases.


These coefficients DO NOT have the same distribution functions. That is, the distribution function for $\rho_{1}$ and $\rho_{2}$ do not overlap. In fact, the latter is lower variance, and shifted to the right of the first. The following figure illustrates that while $\rho_{1}$ has the "expected" Pearson correlation coefficient distribution (left panel), the second coefficient $\rho_{2}$ behaves as though it has more samples, and a higher correlation coefficient (right panel).

Correlation Distribution for Two Cases

A parameter fit for $\rho$ and $n$ using $p_{R}\left( r \,; \rho, n, \,\mathcal{H}_{1}\right)$ illustrates that the distribution for $\rho_{2}$ is parametrized differently (blue curve, right panel). Both $\rho$ and $n$ are larger than expected.

My Question:

Can we develop the correct distribution function for $\rho_{2}$? Or, can we "correct" $\rho$ and $n$ and simply use $p_{R}\left( r \,; \rho, n, \,\mathcal{H}_{1}\right)$?

To generate these figures, Matlab Code follows below:

%NOTE! Function corrdist.m is commented out at bottom of code. Save
%the function as it's own .m file and comment out to run this code.

%STEP 0: 

clear h;
%length of time series/dimension of vector (make N even for convenience)
N       = 1e2;
%number of samples drawn from population (to make histogram). Choose a very
%large sample to make the histogram smooth.
Nsim    = 2e4;

%the true correlation coefficient (||n||^2 = sigma^2*(N-1) = 1*(N-1) ):
r0  = 1/( 1 + (N-1)./(s'*s));

%now add N(0,1) noise to signal to make Nsim sample-vectors from same
%distribution. Call result x.

%signal definition
s   = sqrt(2)*randn(N,1);
x   = repmat(s,1,Nsim) + randn(N, Nsim);

%STEP 1: 

%Define a bivariate, normally distributed random variable z.
R   = chol([1,r0;r0,1]);
mu  = [0,0];

%Compute a correlation coefficient between each random vector
%This is not efficient, but it shows the reader (you) what I am doing:
c   = [];
for k = 1:Nsim,

    z       = (repmat(mu,N,1) +  randn(N,2)*R);
    c       = cat(1, c, z(:,1)'*z(:,2)/(norm(z(:,1))*norm(z(:,2))));    

[Nb, b]     = hist(c, floor(sqrt(length(cc))));
Nb          = Nb./trapz(b, Nb);

h(1) = bar(b, Nb,'facecolor','k','edgecolor','none');
hold on;
r       = linspace(0,1,1e4);
h(2)    = plot(r, corrdist(r, r0, N),'-r','linewidth',4);
f1      = gcf;
a1      = gca;

legend(h,{'Sample Correl. Hist.','Theor. Distr.'},'Location','Best');
legend boxoff;
title('Using Correlated Normal Random Variables');
ylabel('Normalized Density');

%STEP 2: 

%normalize data vectors (columns) for correlation coefficient computation. 
%Take care not to include norm-computations that are below eps to 
%avoid numerical blow-up (divide by zero).
n           = sqrt( sum( x.*conj(x),1) );
m           = n(n>eps);
x(:,n<eps)  = [];
Nsim        = length(m);
temp        = repmat(m.^(-1),size(x,1),1);

%normalize x = signal + noise for correlation computation. 
x           = x.*temp;

%now compute the sample correlation between the first sample vector and the
%following independent and identically distributed column vectors in x.
cc      = x(:,1)'*x(:,2:end);

%compute histogram, and normalized histogram.
[Nb, b]     = hist(cc, floor(sqrt(length(cc))));
Nb          = Nb./trapz(b, Nb);

%NOTE: The correlation distribution function cc and the histogram do not
%agree in this case.
h(3) = bar(b, Nb,'facecolor','k','edgecolor','none');
hold on;
r   = linspace(0,1,1e4);
h(4) = plot(r, corrdist(r, r0, N),'-r','linewidth',4);

%BUT: The EFFECTIVE distribution function for cc behaves as though it has
%many more degrees of freedom, and a different true correlation:

%Find the best fit correlation value and sample number to fit histogram.
dof     = fminsearch( @(p) norm( corrdist(b,p(1),p(2)) - Nb), [r0, N]);
h(5)    = plot(r, corrdist(r, dof(1), dof(2)),'-b','linewidth',4);
f2      = gcf;
a2      = gca;
ylimRef = ylim;
set(a1,'xlim',[0.27, 0.92]);
set(a2,'xlim',[0.27, 0.92]);

legend(h(3:end),{'Correl. Histogram','Theor. Distr.','Best-Fit Distr.'},'Location','Best');
legend boxoff;
title('Using Signal + Noise Vectors');
%Note that the latter,  best-fit correlation distribution fits the 
%histogram MUCH better than the correlation distribution function using the
%assumed parameters, r0 and N.

mydir   = pwd;


sizeX= get(gcf,'Position');
sizeX= sizeX(3:4);


print -dpdf  -r300 CorrelationParadox;

% function y = corrdist(r, ro, n) 
% %This function computes the probability density function for the
% %correlation coefficient of a bivariate random variable.
% %
% % y = corrdist(r, ro, n)
% %
% % r:    Vector of possible correlation random variables, i.e. the values at
% %       which the pdf is evaluated at.
% % ro:   The given (true) correlation coefficient, i.e. the population
% %       correlation coefficient. length(ro) > 1 supported.
% % n:    The number of samples in the correlated data. Only length(n) = 1
% %       supported.
% % 
% % y:    The probability density function for r, given ro, for n data
% %       samples of a bivariate normal distribution.
% %
% %-----------------------------------------------------------------------
% % Latest Edit: 11.June.2012
% % Joshua D Carmichael
% % josh.carmichael@gmail.com
% %
% % Original Author: Xu Cui, Stanford University (retrieved 11.June.2012)
% %-----------------------------------------------------------------------
% %accept vectorized inputs.
% if(length(ro)> 1.5),
%     r   = repmat(r(:),1,length(ro));
%     ro  = repmat(ro(:)', length(r),1);
% end;
% if( n < 120 ),
%     y = (n-2) * gamma(n-1) * ((1-ro.^2).^((n-1)/2)).* (1-r.^2).^((n-4)/2);
%     y = y./ (sqrt(2*pi) * gamma(n-1/2) * (1-ro.*r).^(n-3/2));
%     y = y.* (1+ 1/4*(ro.*r+1)/(2*n-1) + 9/16*(ro.*r+1).^2 / (2*n-1)/(2*n+1));
% else
%     y = (n-2) * (1-ro.^2)^((n-1)/2) * (1-r.^2).^((n-4)/2);
%     y = y./ (sqrt(2*pi) * (1-ro.*r).^(n-3/2)) * n.^(-1/2);
%     y = y.* (1+ 1/4*(ro.*r+1)/(2*n-1) + 9/16*(ro.*r+1).^2 / (2*n-1)/(2*n+1));
% end;
% y(r>1)              = 0;
% y(r<-1)             = 0;
% y(~isfinite(y))     = 0;
% return;
share|improve this question

1 Answer 1

I think I have found the discrepancy.The short answer is that the distribution for case (2) is a doubly noncentral beta distribution. The bulk correlation is the same under each case.

In the solution below, I've used slightly different notation. In particular, $\boldsymbol{s }$ $\rightarrow$ $\boldsymbol{w}$. I've also used $\hat{A}$ to describe the maximum liklehood estimate for the amplitude of the sample vector $\boldsymbol{w}$ that maximally correlates with $\boldsymbol{x}$. This notation and formulation is targeted at detection theory applications.

The relative square error in approximating approximating a data stream $ \boldsymbol{x} $ using a waveform template $\boldsymbol{w }$ is the ratio of the least-squares error $\vert\vert \boldsymbol{e} \vert\vert^{2}$ to measure signal energy $ \vert \vert \boldsymbol{x} \vert \vert^{2}$. This error may be re-written as a ratio of quadratic forms that has a doubly noncentral $\beta$ distribution:

\begin{equation} \begin{split} \cfrac{\vert \vert \boldsymbol{e} \vert \vert^{2}}{\vert \vert \boldsymbol{x} \vert \vert^{2}} &= \cfrac{\vert \vert \boldsymbol{x} - \hat{A} \boldsymbol{w } \vert \vert^{2} }{ \vert \vert \boldsymbol{x} \vert \vert^{2}} \\ &= \cfrac{\vert \vert \boldsymbol{x} - \cfrac{\langle \boldsymbol{x},\, \boldsymbol{w} \rangle}{ \vert \vert \boldsymbol{w } \vert \vert^{2}} \boldsymbol{w } \vert \vert^{2}} {\vert \vert \boldsymbol{x} \vert \vert^{2}} \\ &= \cfrac{ \vert \vert P_{\boldsymbol{w}}^{\perp} \left( \boldsymbol{x} \right) \vert \vert ^{2}} { \vert \vert P_{\boldsymbol{w}}^{\perp}\left( \boldsymbol{x} \right) \vert \vert ^{2} + \vert \vert P_{\boldsymbol{w}} \left( \boldsymbol{x} \right) \vert \vert ^{2}} \\ &\overset{d}{=} \cfrac{ \chi_{1}^{2}( \lambda^{\perp} )} { \chi_{1}^{2}( \lambda ) + \chi_{N_{E} - 1}^{2}( \lambda^{\perp} ) }, \end{split} \end{equation}

where the noncentrality parameters are defined by $\lambda$ $=$ $\cfrac{\vert \vert P_{\boldsymbol{w}} \left( \boldsymbol{x} \right) \vert \vert ^{2}}{\sigma^{2}}$ and $\lambda^{\perp}$ $=$ $\cfrac{\vert \vert P_{\boldsymbol{w}}^{\perp} \left( \boldsymbol{x} \right) \vert \vert ^{2}}{\sigma^{2}}$, and where $\overset{d}{=}$ indicates distributional equality. This ratio is also related to the sample correlation coefficient $r$:

\begin{equation} \begin{split} \cfrac{\vert \vert \boldsymbol{x} - \cfrac{\langle \boldsymbol{x},\, \boldsymbol{w} \rangle}{ \vert \vert \boldsymbol{w } \vert \vert^{2}} \boldsymbol{w } \vert \vert^{2}} {\vert \vert \boldsymbol{x} \vert \vert^{2}} &= 1- \cfrac{\langle \boldsymbol{x},\, \boldsymbol{w} \rangle^{2} }{ \vert \vert \boldsymbol{w } \vert \vert^{2} \vert \vert \boldsymbol{x } \vert \vert^{2}} \\ &= 1 - r^{2} \end{split} \end{equation}


\begin{equation} \begin{split} r^{2} &\overset{d}{=} \cfrac{ \chi_{1}^{2}( \lambda )} { \chi_{1}^{2}( \lambda ) + \chi_{N_{E} - 1}^{2}( \lambda^{\perp} ) } \\ &\sim \text{Beta} \left( \frac{1}{2}, \frac{1}{2}N_{E} ; \lambda, \lambda^{\perp} \right) \end{split} \end{equation}

Distinct hypotheses regarding the distribution for $\boldsymbol{x}$ simplify the form for this distribution. When the data stream contains only noise, the hypothesis $\mathcal{H}_{0}$ is satisfied and $r^{2}$ has a central Beta distribution, where $\lambda^{\perp}$ $=$ $\lambda$ $=$ $0$. In the presence of signal, a data stream $\boldsymbol{x}$ will generally have a non-zero projection $P_{\boldsymbol{w}}^{\perp}\left( \boldsymbol{x} \right)$ orthogonal to the noise-contaminated template data vector $\boldsymbol{w}$. In this case $\lambda$, $\lambda^{\perp}$ $\ne$ $0$, and $r^{2}$ has doubly noncentral Beta distribution. If the template signal has a very large SNR, then $ \boldsymbol{x}$ $\cong$ $A \boldsymbol{w}$ $+$ $\boldsymbol{n}$, $\lambda^{\perp}$ $=$ $0$, and $r^{2}$ is reasonably approximated by a noncentral Beta distribution. The noncentral Beta distribution therefore provides an absolute upper bound on the detection performance of a correlation detector.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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