Correlation of the sigmoid function of normal random varaibles

Question

Suppose $X=(X_1,X_2)^T \sim N(\mu, \Sigma)$ and $p=(p_1,p_2)^T = (e^{X_1}/(1+e^{X_1}),e^{X_2}/(1+e^{X_2}))^T$. $Cov(X_1,X_2)=\sigma_{12}$. Does the following hold ?

1. $\sigma_{12}=0 \Leftrightarrow Cov(p1,p2)=0$
2. $\sigma_{12} > 0 \Leftrightarrow Cov(p1,p2) > 0$
3. $\sigma_{12} < 0 \Leftrightarrow Cov(p1,p2) < 0$

We only have the result that two uncorrelated normal random variables are independent, hence $\sigma_{12}=0 \Rightarrow Cov(p1,p2)=0$. I did several simulations and found no counter examples for the rest.

If we can claim that $p_1, p_2$ are independent when they are uncorrelated. Then first equivalence is immediate. So question is what additional conditions are needed for two uncorrelated r.v. to be independent ?

I am also trying to evaluate the following integrals to see if I can get some ideas(the first is solved, but the second is not).

$$\int_{-\infty}^{+\infty} \frac{1}{1+e^x} e^{-x^2/2} dx$$

$$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \frac{1}{1+e^x} \frac{1}{1+e^y}e^{-\frac{x^2+y^2-2\rho xy}{2(1-\rho^2)}} dx dy$$

EDIT I plot the covariance against $\rho$: $$Cov(p1,p2)= \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \frac{1}{1+e^x} \frac{1}{1+e^y} \frac{1}{2\pi\sqrt{1-\rho^2}}e^{-\frac{x^2+y^2-2\rho xy}{2(1-\rho^2)}} dx dy-(\int_{-\infty}^{+\infty} \frac{1}{\sqrt{2\pi}}\frac{1}{1+e^x} e^{-x^2/2} dx)^2 $$.

To my surprise, the curve is a linear function of $\rho$. And those three equivalence hold under this special case.

Answer 1

What you have could be called an multivariate logistic-normal distribution, but it is not the same as the multivariate distribution with that name discussed in https://en.wikipedia.org/wiki/Logit-normal_distribution If $X$ is multinormal, the multivariate logit normal distribution we discuss here is the distribution of $Y$, where each component of $Y$ is the logistic transform of the corresponding component of $X$, $$ Y_i = l(X_i) = \frac{e^{X_i}}{1+e^{X_i}}= \frac1{1+e^{-X_i}} $$ By calculating the relevant jacobian (which is simple since in this case is the determinant of a diagonal matrix) we can find the density of $Y$ as $$ f(y)=\frac{\exp\left(-\frac12 (\log\frac{y}{1-y}-\mu_1)^T\Sigma^{-1}(\log\frac{y}{1-y}-\mu_2) \right)}{(2\pi)^{d/2}|\Sigma |^{1/2}\prod_{j=1}^d(y_j(1-y_j))} $$ where the parameters are the parameters of the multinormal distribution, $d$ is the dimension of $X$ and $Y$, and the log function is applied componentwise. A direct answer to your questions would be easy if we could find the moments in explicit form, but that seems impossible. Let us have a look at an example of this density in the bivariate case, with

$\mu_1=\mu_2=0, \sigma_1=\sigma_2=1$ and $\rho=0.5$. This case is exchangeable, which shows in that the density is symmetric around the line $y_1=y_2$. One way to investigate the correlation is by numerical integration, and for this symmetric case, but with varying $\rho$, this is done below, and the result is that the correlation is remarkably little influenced by the transformation, so to a good approximation the correlation in $Y$ is the same as that in $X$. Some R code:

mu  <-  c(0, 0)
sigma1  <-  1.0
sigma2 <-  1.0
rho <-  0.5

dmvlogist_norm  <-  function(x, y, mu, sigma1, sigma2, rho) {
    const  <- 2*pi*sigma1*sigma2*sqrt(1-rho^2)*y*x*(1-y)*(1-x)
    d <- exp(-(1/(2*(1-rho^2)))*( ((log(x/(1-x))-mu[1])/sigma1)^2 +
                                  ((log(y/(1-y))-mu[2])/sigma2)^2 -
                                  2*rho*(log(x/(1-x))-mu[1])*
                                  (log(y/(1-y))-mu[2])/(sigma1*sigma2) ) )/const
    return( d )
}

# Code for the contour plot:

x <- seq(0.0001,  0.9999,  length=101)
xy <- as.matrix(expand.grid(x, x))
d  <-  dmvlogist_norm(xy[, 1], xy[, 2], mu, sigma1, sigma2, rho)

d  <-  matrix(d, 101, 101)

image(x, x, d)
contour(x, x, d, nlevels=20, add=TRUE)
title("bivariate logit normal density")

Since there is little hope of finding the correlation via symbolic integration, we use numerical methods:

library(pracma)

calc_corr  <-  function(rho) {# we can simplify code using exchangeability
    m1  <-  m2 <- integral2(function(x, y) x*dmvlogist_norm(x, y, c(0, 0),
                                                       1, 1, rho), 0, 1, 0, 1)$Q
    s1 <- s2 <-  integral2(function(x, y) (x-m1)^2*dmvlogist_norm(x, y, c(0, 0),
                                                              1, 1, rho), 0, 1, 0, 1)$Q
    s12 <-  integral2(function(x, y) (x-m1)*(y-m2)*dmvlogist_norm(x, y, c(0, 0),
                                                                   1, 1, rho), 0, 1, 0, 1)$Q
    rhoout  <-  s12/s1 # using exchangeability
    rhoout
}

rho <-  seq(-0.99, 0.99, length=31)
rho_logist <- rho
for (r in seq(along=rho)) rho_logist[r] <- calc_corr(rho[r])

The results as a table:

> cbind(rho, rho_logist, rho_logist-rho)
         rho    rho_logist              
 [1,] -0.990 -9.320447e-01  5.795527e-02
 [2,] -0.924 -9.218907e-01  2.109328e-03
 [3,] -0.858 -8.541739e-01  3.826148e-03
 [4,] -0.792 -7.871561e-01  4.843921e-03
 [5,] -0.726 -7.203866e-01  5.613378e-03
 [6,] -0.660 -6.539386e-01  6.061385e-03
 [7,] -0.594 -5.877733e-01  6.226669e-03
 [8,] -0.528 -5.218579e-01  6.142147e-03
 [9,] -0.462 -4.561610e-01  5.838953e-03
[10,] -0.396 -3.906522e-01  5.347779e-03
[11,] -0.330 -3.253032e-01  4.696777e-03
[12,] -0.264 -2.600864e-01  3.913566e-03
[13,] -0.198 -1.949740e-01  3.025960e-03
[14,] -0.132 -1.299396e-01  2.060400e-03
[15,] -0.066 -6.495576e-02  1.044244e-03
[16,]  0.000  3.711409e-08  3.711409e-08
[17,]  0.066  6.495857e-02 -1.041433e-03
[18,]  0.132  1.299417e-01 -2.058303e-03
[19,]  0.198  1.949761e-01 -3.023938e-03
[20,]  0.264  2.600891e-01 -3.910861e-03
[21,]  0.330  3.253074e-01 -4.692627e-03
[22,]  0.396  3.906578e-01 -5.342230e-03
[23,]  0.462  4.561665e-01 -5.833508e-03
[24,]  0.528  5.218657e-01 -6.134344e-03
[25,]  0.594  5.877817e-01 -6.218258e-03
[26,]  0.660  6.539465e-01 -6.053502e-03
[27,]  0.726  7.203924e-01 -5.607636e-03
[28,]  0.792  7.871625e-01 -4.837515e-03
[29,]  0.858  8.542809e-01 -3.719148e-03
[30,]  0.924  9.218908e-01 -2.109192e-03
[31,]  0.990  9.909956e-01  9.956272e-04

In the third column the difference which is mainly in the third decimal.