Estimators of the correlation coefficient in the bivariate exponential distribution

TL;DR

This paper derives a lower bound on the estimation error of the correlation coefficient in bivariate exponential distributions and evaluates the efficiency of three nonlinear estimators, highlighting their performance across different correlation ranges.

Contribution

It introduces a finite-support parameter constraint to establish a lower bound and compares the optimality of three nonlinear estimators for the correlation coefficient.

Findings

The cosine similarity-based estimator is highly efficient for correlation > 0.35.

The transformed Pearson correlation performs better for smaller correlation values.

A lower bound on estimation error is derived under the finite-support constraint.

Abstract

A finite-support constraint on the parameter space is used to derive a lower bound on the error of an estimator of the correlation coefficient in the bivariate exponential distribution. The bound is then exploited to examine optimality of three estimators, each being a nonlinear function of moments of exponential or Rayleigh observables. The estimator based on a measure of cosine similarity is shown to be highly efficient for values of the correlation coefficient greater than 0.35; for smaller values, however, it is the transformed Pearson correlation coefficient that exhibits errors closer to the derived bound.

Equations

$${\mathcal{C}}_{XY}=\left[\begin{array}{cccc}{\sigma }_X^2& 0& {\sigma }_X{\sigma }_Y{\rho }_c& {\sigma }_X{\sigma }_Y{\rho }_s\\ 0& {\sigma }_X^2& -{\sigma }_X{\sigma }_Y{\rho }_s& {\sigma }_X{\sigma }_Y{\rho }_c\\ {\sigma }_X{\sigma }_Y{\rho }_c& -{\sigma }_X{\sigma }_Y{\rho }_s& {\sigma }_Y^2& 0\\ {\sigma }_X{\sigma }_Y{\rho }_s& {\sigma }_X{\sigma }_Y{\rho }_c& 0& {\sigma }_Y^2\end{array}\right]$$

$$V=\sqrt{X_I^2+X_Q^2}\mathrm{and}Z=\sqrt{Y_I^2+Y_Q^2}$$

$$p_{VZ}(v,z)=\frac{vz}{{\sigma }_X^2{\sigma }_Y^2(1-{\rho }^2)}\exp \left[-\frac{1}{2(1-{\rho }^2)}\left(\frac{v^2}{{\sigma }_X^2}+\frac{z^2}{{\sigma }_Y^2}\right)\right]$$

$$\times I_0\left[\frac{\rho vz}{{\sigma }_X{\sigma }_Y(1-{\rho }^2)}\right],v,z\ge 0,\rho \ge 0$$

$$\mathrm{E}\{V\}={\sigma }_X\sqrt{\pi /2},\mathrm{E}\{Z\}={\sigma }_Y\sqrt{\pi /2}$$

$$\mathrm{E}\{V^2\}=2{\sigma }_X^2,\mathrm{E}\{Z^2\}=2{\sigma }_Y^2$$

$$\mathrm{E}\{VZ\}={\sigma }_X{\sigma }_Y\left[2\text{mspace}1.0mu\mathsf{E}(\rho )-(1-{\rho }^2)\mathsf{K}(\rho )\right]$$

$$\mathsf{K}(0)=\mathsf{E}(0)=\pi /2,\mathsf{E}(1)=1$$

$$U=V^2\mathrm{and}W=Z^2$$

$$p_{UW}(u,w)=\frac{1}{4{\sigma }_X^2{\sigma }_Y^2(1-r)}\exp \left[-\frac{1}{2(1-r)}\left(\frac{u}{{\sigma }_X^2}+\frac{w}{{\sigma }_Y^2}\right)\right]$$

$$\times I_0\left[\frac{\sqrt{ruw}}{{\sigma }_X{\sigma }_Y(1-r)}\right],u,w\ge 0,r\ge 0$$

$$\mathrm{E}\{U\}=2{\sigma }_X^2,\mathrm{E}\{W\}=2{\sigma }_Y^2$$

$$\mathrm{E}\{U^2\}=8{\sigma }_X^4,\mathrm{E}\{W^2\}=8{\sigma }_Y^4$$

$$\mathrm{E}\{UW\}=4(r+1){\sigma }_X^2{\sigma }_Y^2.$$

$$m_{E\kappa \nu }\triangleq \frac{1}{n}\sum _{i=1}^n{u}_i^{\kappa }{w}_i^{\nu }\mathrm{or}{m}_{R\kappa \nu }\triangleq \frac{1}{n}\sum _{i=1}^n{v}_i^{\kappa }{z}_i^{\nu }$$

$$\mathbf{\theta }\triangleq ({\theta }_1,{\theta }_2,{\theta }_3{)}^T\equiv (r,{\sigma }_X^2,{\sigma }_Y^2{)}^T.$$

$${\left[\mathcal{I}(\mathbf{\theta })\right]}_{k,\ell }\triangleq \mathrm{E}\left\{\frac{\partial }{\partial {\theta }_k}\ln p_{UW}(u,w)\frac{\partial }{\partial {\theta }_{\ell }}\ln p_{UW}(u,w)\right\},$$

$$k,\ell =1,2,3$$

$$\mathrm{var}\{\hat{R}\}\le \frac{1}{n}{\left[{\mathcal{I}}^{-1}(\mathbf{\theta })\right]}_{1,1}\equiv {\sigma }_{CR}^2(r),$$

$$\hat{r}_{m}\,\,=\,\,\left\{\begin{array}[]{cl}{\!\!\hat{r},}&{\hat{r}\geq 0}\\ {\!\!0,}&{\hat{r}<0}.\end{array}\right.$$

$$p_{{\hat{R}}_m}({\hat{r}}_m;r)=\gamma \delta ({\hat{r}}_m)+p_{\hat{R}}({\hat{r}}_m;r)1({\hat{r}}_m)$$

$$\gamma ={\int }_{-\infty }^0{p}_{\hat{R}}(\hat{r};r)d\hat{r}$$

$$\mathbf{1}(\omega)\,\,\triangleq\,\,\left\{\begin{array}[]{cl}{\!\!1,}&{\omega\geq 0}\\ {\!\!0,}&{\omega<0}.\end{array}\right.$$

$${\epsilon }_{MS}^2(r)\triangleq \mathrm{E}\{({\hat{R}}_m-r{)}^2\}=\mathrm{var}\{{\hat{R}}_m\}+\underset{\mathrm{bias}\mathrm{squared}}{\underbrace{(\mathrm{E}\{{\hat{R}}_m\}-r)}}\text{mspace}-1.0mu{}^2.$$

$${\epsilon }_{MS}^2(r)=\underset{\mathrm{variance}}{\underbrace{{\sigma }_{CR}^{2}F(\mu )\left[(1-d)+F(-\mu )(\mu +h{)}^2\right]}}+$$

$$+\underset{\mathrm{bias}\mathrm{squared}}{\underbrace{[F(\mu )\text{mspace}1.0mu(r+h\text{mspace}1.0mu{\sigma }_{CR})-r]}}{}^2$$

$${\sigma }_{CR}^2(0)/2\le {\epsilon }_{MS}^2(r)<{\sigma }_{CR}^2(r)$$

$$r_P(U,W)\triangleq \frac{\mathrm{E}\{UW\}-\mathrm{E}\{U\}\mathrm{E}\{W\}}{\sqrt{\mathrm{var}\{U\}\mathrm{var}\{W\}}}.$$

$$\mathsf{s}(\mathbf{u},\mathbf{w})=\frac{m_{E11}-m_{E10}{m}_{E01}}{\sqrt{(m_{E20}-m_{E10}^2)(m_{E02}-m_{E01}^2)}}$$

$$\hat{r}_{1}\,\,=\,\,\left\{\begin{array}[]{cl}{\!\!\mathsf{s}(\mathbf{u},\mathbf{w}),}&{\mathsf{s}(\mathbf{u},\mathbf{w})\geq 0}\\ {\!\!0,}&{\text{otherwise.}}\end{array}\right.$$

$$r_P(V,Z)=\frac{2\left[2\text{mspace}1.0mu\mathsf{E}(\sqrt{r}\text{mspace}1.0mu)-(1-r)\mathsf{K}(\sqrt{r}\text{mspace}1.0mu)\right]-\pi }{4-\pi }.$$

$$\xi (\mathbf{v},\mathbf{z})=\mathsf{s}(\mathbf{v},\mathbf{z})\left\{1+g\text{mspace}1.0mu[1-\mathsf{s}(\mathbf{v},\mathbf{z})]\right\},g=49/500$$

$$\hat{r}_{2}\,\,=\,\,\left\{\begin{array}[]{cl}{\!\!\xi(\mathbf{v},\mathbf{z}),}&{\xi(\mathbf{v},\mathbf{z})\geq 0}\\ {\!\!0,}&{\text{otherwise.}}\end{array}\right.$$

$${\mathsf{c}}^2(\mathbf{v},\mathbf{z})\triangleq \frac{m_{R11}^2}{m_{R20}{m}_{R02}}.$$

Full text

I Introduction

The bivariate exponential and Rayleigh probability distributions, [1, pp. 401–475], [2], play a prominent role in the development of models of dependent nondeterministic phenomena in science and engineering. Such models include power of a random signal received at multiple sensors exploiting time/space/frequency diversity, weather radar returns observed at co-polar and cross-polar channels, weights of edges in random graphs being matched, time intervals between significant events occuring in parts of a complex biological or man-made system and many more.

The statistical association between observables of interest can be characterized by exploiting various measures of dependence, such as mutual information, copulas and parametric or nonparametric correlation coefficients [1, pp. 105–177], [3]. In practice, the correlation coefficient appears to be a preferred choice owing to its computational simplicity, and also the fact that it can be functionally related to copulas and mutual information [3], [4].

The problem of estimating the correlation coefficient between non-negative observables has been discussed in a number of publications [5]–[8]. However, since the finite-support constraint on the parameter space has been ignored, no conclusions regarding optimality of proposed estimators could be drawn. Therefore, it is of interest to establish a constrained lower bound on the estimator error and examine estimators that could attain this bound.

II Rayleigh and Exponential Distributions

Consider two complex Gaussian random variables (rvs), $\mathbf{X}\triangleq X_{I}+jX_{Q}$ and $\mathbf{Y}\triangleq Y_{I}+jY_{Q}$, where $j^{2}=-1$. The four jointly Gaussian components, $(X_{I},X_{Q},Y_{I},Y_{Q})$, have all zero means, $\mathrm{E}\{X_{I}\}=\mathrm{E}\{X_{Q}\}=\mathrm{E}\{Y_{I}\}=\mathrm{E}\{Y_{Q}\}~{}=0$, where $\mathrm{E}\{\cdot\}$ denotes expectation, and their covariance matrix is of the form [9]

[TABLE]

where $|\rho_{c}|\leq 1$ and $|\rho_{s}|\leq 1$ are correlation coefficients between respective rvs.

In signal processing, the complex Gaussian rvs, $\mathbf{X}$ and $\mathbf{Y}$, may be viewed as discrete-time samples of two dependent complex Gaussian processes ${\mathbf{X}}(t)$ and ${\mathbf{Y}}(t)$. The rvs, $\mathbf{X}$ and $\mathbf{Y}$, may also represent samples, taken at different time instants, say, $t$ and $t+\tau$, of a single stationary complex Gaussian process ${\mathbf{X}}(t)$; in such a case, ${\mathbf{Y}}(t)={\mathbf{X}}(t+\tau)$ and $\sigma^{2}_{Y}=\sigma^{2}_{X}$.

II-A Bivariate Rayleigh Distribution

Pairs of rvs, $(X_{I},X_{Q})$ and $(Y_{I},Y_{Q})$, can be used to construct two Rayleigh rvs, $V$ and $Z$, as follows

[TABLE]

The rvs $V$ and $Z$ represent magnitudes of the corresponding underlying complex Gaussian rvs $\mathbf{X}$ and $\mathbf{Y}$.

The joint probability density function (pdf) of $V$ and $Z$ is given by [2]

[TABLE]

where $\rho^{2}=\rho^{2}_{c}+\rho^{2}_{s}$, and $I_{0}(\cdot)$ denotes a modified Bessel function of the first kind of order zero. If $\rho=0$, then $p_{VZ}(v,z)\!=\!p_{V}(v)\mspace{2.0mu}p_{Z}(z)$, where $p_{V}(v)$ and $p_{Z}(z)$ are marginal Rayleigh pdfs of $V$ and $Z$, respectively. Therefore, in this case, zero correlation implies statistical independence.

Population joint moments, $\mathrm{E}\{V^{\kappa}Z^{\nu}\},\kappa\!+\!\nu=1,2$, of rvs $V$ and $Z$ are given by [9]

[TABLE]

where ${\mathsf{K}}(\cdot)$ and ${\mathsf{E}}(\cdot)$ are complete elliptic integrals of the first and second kind. In particular,

[TABLE]

and when $\rho$ approaches one, ${\mathsf{K}}(\rho)$ tends to infinity.

II-B Bivariate Exponential Distribution

The transformation

[TABLE]

converts two Rayleigh rvs, $V$ and $Z$, into two exponential rvs, $U$ and $W$. The joint pdf of $U$ and $W$ can be expressed as [2]

[TABLE]

where $r=\rho^{2}$. The parameter $r$ is, in fact, the correlation coefficient between exponential rvs $U$ and $W$ (see Section V). Also in this case, when $r=0$, rvs $U$ and $W$ are statistically independent.

Population joint moments, $\mathrm{E}\{U^{\kappa}W^{\nu}\},\kappa\!+\!\nu=1,2$, of rvs $U$ and $W$ are given by [9]

[TABLE]

III Problem Formulation

Assume that observations on rvs $U$ and $W$ are made in pairs, $(\mathbf{u},\mathbf{w})=\{(u_{i},w_{i}):i=1,2,\ldots,n\}$; alternatively, observations, $(\mathbf{v},\mathbf{z})=\{(v_{i},z_{i}):i=1,2,\ldots,n\}$, may be made on Rayleigh rvs $V$ and $Z$. Then, $n$ pairs of observations are used to determine sample joint moments,

[TABLE]

corresponding, respectively, to population moments (8) or (4).

This Letter addresses two associated problems:

1. Given the pdf (7) and the constraint, $0\leq r\leq 1$, derive a lower bound on the error of an estimator of the correlation coefficient $r$ appearing in (7).

2. Make use of sample moments (9) to construct estimators of $r$ and examine their optimality with respect to the derived lower bound.

IV Lower Bounds on Estimation Errors

In the case of a bivariate exponential distribution (7), allowed values of the correlation coefficient $r$ are restricted to the $(0,1)$-interval. If a statistic employed as an estimator of $r$ assumes values from a different, finite or infinite, interval, then the constraint, $0\leq r\leq 1$ must be taken into account when establishing a lower bound on the estimator error.

IV-A Cramér-Rao Bound (CRB)

It is known [10] that under suitable regularity conditions, the variance of any unbiased estimator can be bounded by the lower Cramér-Rao bound (CRB). Therefore, the CRB is a useful measure when examining optimality of several competing estimators of a parameter of interest.

Let a vector $\boldsymbol{\theta}$ of nonrandom parameters be defined by

[TABLE]

Then (neglecting any constraints on the parameters), the Fisher information matrix, $\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})$, is a $3\!\times\!3$ positive semidefinite symmetric matrix, comprising the elements

[TABLE]

Consequently, a lower bound on the variance of any unbiased estimator $\hat{R}$ of $r$ can be determined from

[TABLE]

where $\boldsymbol{\mathcal{I}}^{-1}$ is the inverse of $\boldsymbol{\mathcal{I}}$.

Elements of the Fisher information matrix $\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})$, for selected values of $r$, are given in [11]. Values of the Cramér-Rao bound, shown in Table 1, have been determined by selecting a first diagonal element of the inverse $\boldsymbol{\mathcal{I}}^{-1}(\boldsymbol{\theta})$ of $\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})$.

IV-B Mean-Square-Error (MSE) Bound

When the parameter space is restricted, the Cramér-Rao approach appears to be inadequate [12]–[14]. Therefore, to determine a lower bound on the error of an estimator of parameter $r$ in (7), knowledge of the finite-support constraint, $0\leq r\leq 1$, should be suitably combined with Fisher information contained in available data.

Consider an unbiased estimator $\hat{R}$ of $r$ and let $p_{\hat{R}}(\hat{r};r)$ be a pdf of $\hat{R}$. Assume that the estimator $\hat{R}$ is so constructed that values of its realizations (estimates) $\hat{r}$ cannot exceed one. However, depending on a set of processed data, $\{(u_{i},w_{i}):i=1,2,\ldots,n\}$ or $\{(v_{i},z_{i}):i=1,2,\ldots,n\}$, some estimates $\hat{r}$ may assume negative, hence not allowed values.

Therefore, when such an aberrant estimate $\hat{r}$ is observed, its value must be set to zero, and so modified estimate, $\hat{r}_{m}$, will assume the form

[TABLE]

Consequently, the pdf of the modified estimator ${\hat{R}}_{m}$ will become a censored distribution [15],

[TABLE]

comprising a discrete probability mass and a continuous part. In (14), $\delta(\cdot)$ is an impulse (Dirac delta) function, $\gamma$ is the probability that ${\hat{R}}<0$,

[TABLE]

and $\mathbf{1}(\cdot)$ denotes the Heaviside step function,

[TABLE]

Fig. 1 illustrates the effect of transforming the pdf of an estimator $\hat{R}$ into its censored version, $p_{{\hat{R}}_{m}}({\hat{r}}_{m};r)$, when the value of the correlation coefficient $r$ being estimated decreases from $r_{\beta}$ to $r_{\alpha}$.

The mean-square error (MSE), $\varepsilon^{2}_{MS}(r)$, of the modified estimator $\hat{R}_{m}$ can be expressed as

[TABLE]

In order to determine the lower MSE bound, assume that $\hat{R}$ is a maximum-likelihood (ML) estimator. Since ML estimators are known to be asymptotically unbiased, efficient and Gaussian [10], let $\hat{R}\sim\mathcal{N}(r,\sigma^{2}_{CR})$. The MS error of a modified estimator $\hat{R}_{m}$ can be evaluated by exploiting moments of a censored Gaussian distribution [15].

Let $\phi(\lambda)=(1/\sqrt{2\pi})\exp(-\lambda^{2}/2)$ be the pdf of a standard Gaussian rv $\Lambda$, and $F(\mu)=\Pr\{\Lambda\leq\mu\}$ its cumulative distribution function. Then, the MSE of the estimator $\hat{R}_{m}$ can be expressed as follows

[TABLE]

where $\sigma^{2}_{CR}\equiv\sigma^{2}_{CR}(r),\,\,\mu=r/\sigma_{CR},\,\,h=\phi(\mu)/F(\mu)$ and $d=h(h+\mu)$.

The constrained error bound (17) differs from the CR bound (12), when $r$ is less than approximately $3\,\sigma_{CR}(r)\,$. In the region, $0\leq r<3\,\sigma_{CR}(r)$, the estimator $\hat{R}_{m}$ becomes biased, and its MS error,

[TABLE]

remains below the CR bound. The bound reduction has resulted from incorporating knowledge of the constraint.

V Estimators of the Correlation Coefficient

Consider the population Pearson product-moment correlation coefficient defined by

[TABLE]

By inserting moments (8) into (19), it can be verified that $r_{P}(U,W)=r$. Therefore, the sample Pearson correlation coefficient, i.e. the statistic

[TABLE]

can be used to construct a censored estimate $\hat{r}_{1}$ of $r$ as follows

[TABLE]

When the number $n$ of observations tends to infinity, sample moments converge to population moments, and the sample correlation coefficient $\mathsf{s}(\mathbf{u},\mathbf{w})$ will approach $r$.

The use of sample correlation coefficient to estimate a population correlation coefficient is a standard practice. However, such an approach may not necessarily lead to an efficient estimator (an estimator whose variance attains the Cramér-Rao bound), especially in small or moderate sample sizes.

V-A Estimator Based on Correlation of Rayleigh Variables

Consider now the bivariate Rayleigh distribution (3) and the population Pearson correlation coefficient $r_{P}(V,Z)$, given by a formula analogous to (19). The correlation coefficient $r_{P}(V,Z)$ can be expressed in terms of moments (4) as follows

[TABLE]

In this case, $r_{P}(V,Z)=r$, only when $r=0$ or $r=1$; otherwise, $r_{P}(V,Z)$ is a nonlinear function of $r$.

When $n\rightarrow\infty$, the sample correlation coefficient $\mathsf{s}(\mathbf{v},\mathbf{z})$ will approach (22). By employing the nonlinear transformation

[TABLE]

a censored estimate $\hat{r}_{2}$ of $r$ is obtained as

[TABLE]

V-B Approximate Maximum-Likelihood Estimator

It has been shown [16] that in a case of highly correlated Rayleigh rvs, and when $\sigma_{X}=\sigma_{Y}$, an approximate ML estimator of the correlation coefficient $\rho$ is of the form, $[\mspace{1.0mu}2m_{R11}/(m_{R20}+m_{R02})]^{2}$. The constraint, $\sigma_{X}=\sigma_{Y}$, can be removed by employing the geometric mean rather than the arithmetic mean. Consequently, the following statistic of cosine-similarity-squared is obtained

[TABLE]

The statistic (25) asymptotically converges to

[TABLE]

For $r=0$ and $r=1$, the respective limits are $\pi^{2}/16$ and $1$.

When the nonlinear transformation

[TABLE]

is applied, a censored estimate $\hat{r}_{3}$ of $r$ assumes the form

[TABLE]

Owing to its origin, the estimator $\hat{R}_{3}$ is expected to be asymptotically efficient, at least for larger values of $r$.

V-C Performance of the Estimators

Computer simulations were employed to examine the performance of the three estimators, $\hat{R}_{1},\,\hat{R}_{2}$ and $\hat{R}_{3}$, of the correlation coefficient $r$. Three sample sizes, $n=10,\,n=50$ and $n=200$, were chosen, somewhat arbitrarily, to represent the cases of small, moderate, and large sample sizes. Values of the correlation coefficient $r$ to be estimated varied from $r=0$ to $r=0.98$, in steps of $0.02$. For each combination of $n$ and $r$, $10^{6}$ Monte Carlo experiment replications were carried out to determine the MS error, $\varepsilon^{2}$, for each of the three estimators.

Results of the study are shown in Fig. 2 along with the MSE bound (17) and the Cramér-Rao bound (12); values of the MSE bound are only shown when they differ from those of the CR bound.

The results can be summarized as follows:

1. The derived MSE lower bound is superior to the standard CRB when predicting errors of estimators of the correlation coefficient $r$; the MSE lower bound is more precise when the sample size is moderate or large.

2. When $r$ is greater than $r^{*}\approx 0.35$, the estimator $\hat{R}_{3}$ is better than the other two estimators, and its estimated MS error, $\hat{\varepsilon}^{2}$, differs only slightly from the derived lower bound.

3. In the region, $r<r^{*}$, the estimator $\hat{R}_{2}$ is superior to the estimator $\hat{R}_{3}$.

4. When $r\approx 0$, the estimator $\hat{R}_{1}$ exhibits the smallest MS error; this observation supports the conclusion in [8] that the sample correlation coefficient (20) is an asymptotically most powerful test of the hypothesis $r=0$ against the alternative $r>0$.

5. When $r<0.1$, the MS error of the estimator $\hat{R}_{3}$ markedly exceeds those of the other two estimators; this effect can partly be attributed to the approximate nature of the nonlinearity (27).

VI Conclusion

The non-negativity constraint has been incorporated into the standard CR bounding technique by utilizing moments of a censored Gaussian distribution. The resulting MSE bound establishes a lower bound on the MS error of any estimator of the correlation coefficient of exponentially distributed variables.

The simulation study has shown that MS errors associated with two of the examined estimators are close to the derived lower bound in two subintervals that jointly cover the entire (0,1)-interval. Each of the two estimators is a nonlinear function of a measure of either cosine similarity or centred cosine similarity (i.e. the sample correlation coefficient) between Rayleigh variables.