Correlation estimation on Half-Normal Distribution

Question

Let $$(X,Y)\sim N\left(\begin{pmatrix}0\\0\end{pmatrix},\begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\right)$$ and let we are observing $(|X_1|,|Y_1|),\dots,(|X_n|,|Y_n|)$ independently. I wish to estimate $\rho$ given this observation. Is it possible? I worry I can't because what we observe are the absolute value.

Answer 1

Because taking the absolute values is an area-preserving transformation of the plane that maps four quadrants onto one, the likelihood of $\rho$ for data $(x,y)$ is the sum of the likelihoods under a Normal$\left((0,0), \pmatrix{1&\rho\\\rho&1}\right)$ distribution of the four data values $(\pm x, \pm y)$:

$$\lambda(\rho;x,y)=\frac{1}{2 \pi \sqrt{1-\rho ^2}} \left(2\exp\left(\frac{x^2-2 \rho x y+y^2}{2 \left(\rho ^2-1\right)}\right) + 2\exp\left(\frac{x^2+2 \rho x y+y^2}{2 \left(\rho ^2-1\right)}\right)\right).$$

The maximum likelihood estimator works well for extreme values of $\rho$ but is highly uncertain for small to moderate values, as attested by these simulations (of $1000$ iterations of a sample of $100$ values). Even visually the plot of $100$ points with $\rho=1/4$ is not readily distinguishable from the plot with $\rho=2/3$:

Each row depicts the situation with a different value of $\rho$. The vertical lines in the "Likelihood" plots show the locations of the maxima in the samples depicted at left; these are the maximum likelihood estimates (MLEs). The vertical lines in the "Estimation Distribution" plots of the MLEs for $1000$ additional simulated data sets are located at the true value of $\rho$ held constant throughout each simulation. Any value of $\rho$ for which the likelihood is larger than the dotted gray line can be considered within a two-sided $95\%$ confidence interval.

The R code that produced these results can readily be adapted for further study of the Maximum Likelihood estimator. In particular, the estimation distribution for small values of $\rho$ is bimodal and the MLE appears to be biased a little low: this is worth understanding well before relying on the MLE.

library(mvtnorm) # exports rmvnorm()
likelihood.log <- function(rho, xy) {
  const <- (dim(xy)[1]/2) * log((2*pi)^2 * (1-rho^2))
  sum(log(2*exp((xy[,1]^2 + xy[,2]^2 - 2*xy[,1]*xy[,2]*rho) / (2*(-1 + rho^2))) +
        2*exp((xy[,1]^2 + xy[,2]^2 + 2*xy[,1]*xy[,2]*rho) / (2*(-1 + rho^2))))) - const
}

half.normal <- function(rho.0=0, n=100, alpha=0.05, n.iter=1000) {
  # rho.0:  True value of rho
  # n    :  Sample size
  # alpha:  Determines gray line (for the confidence interval)
  # n.iter: Number of simulation iterations
  #
  # Simulate a dataset and plot it.
  #
  xy <- abs(rmvnorm(n, mean=c(0,0), sigma=matrix(c(1,rho.0,rho.0,1),2)))
  plot(xy, main=paste("Data (rho=", round(rho.0,2), ")", sep=""))
  #
  # Plot the log likelihood function for this dataset.
  #
  rho.min <- tanh(atanh(rho.0)-1/2)
  rho.max <- tanh(atanh(rho.0)+1/2)
  rho <- seq(rho.min, rho.max,length.out=511)
  l <- sapply(rho, function(rho) likelihood.log(rho, xy))
  i <- l >= max(l) - 8

  plot(rho[i], l[i],
       type="l", xlab="rho", ylab="Log likelihood",
       main="Likelihood")
  fit.ml <- optimize(function(rho) likelihood.log(rho, xy), interval=c(0,1),
            maximum=TRUE)
  abline(v = fit.ml$maximum, col="Red")
  abline(h = fit.ml$objective - qchisq(1-alpha, 1), lty=3, col="Gray")
  #
  # Simulate additional data to study the estimation distribution.
  #
  sim <- replicate(n.iter, {
    xy <- abs(rmvnorm(n, mean=c(0,0), sigma=matrix(c(1,rho.0,rho.0,1),2)))
    fit.ml <- optimize(function(rho) likelihood.log(rho, xy), interval=c(0,1),
              maximum=TRUE)
    fit.ml$maximum
  })
  abline(v = fit.ml$maximum, col="Red")
  hist(sim, main="Estimation Distribution")
  abline(v = abs(rho.0), col="Red")
  #
  # Return the average estimate (to assess bias).
  #
  mean(sim)
}

set.seed(17)
par(mfrow=c(3,3))
half.normal(1/4)
half.normal(2/3)
half.normal(0.99)

Answer 2

Here's a relatively simple* approach I've been thinking about.

You can work out that (in this particular instance) $\text{corr}(X^2,Y^2)=\rho^2$.

(I can give the derivation if needed).

So one quick and easy way to estimate $\rho$ would be to square your absolute values, and take the square root of the correlation: $\hat{\rho}=\sqrt{\text{corr}(X^2,Y^2)}$.

Note that this estimator won't be unbiased. Indeed, in with small correlations (as noted below), it has a problem, in that the correlation might be negative. In large samples, (at least with large correlation) it does moderately well.

Here's an example simulation:

With 10000 samples of size 100 with $\rho=0.8$, the mean correlation estimate was 0.79, and the distribution of correlation estimates looked like this:

This seems to perform adequately*, at least with reasonably large correlations.

* if not very efficiently compared to ML estimation.

As whuber notes in comments, there's a problem with smaller $\rho$ of negative sample correlation between the squares; indeed that becomes fairly likely. Since we have specified that $\rho$ will be non-negative, his suggestion of replacing negatives with 0 is a pretty good one, but if the correlations aren't reasonably large I wouldn't recommend this approach.