Let's say you have a jointly gaussian vector random variable $\mathbf{x}$, with mean $\mathbf{M}$ and covariance $\mathbf{S}$. I now transform each scalar element of $\mathbf{x}$ , say $x_j$, with a sigmoid:

$$y_j = 1/(1+\exp(-x_j))$$

I am interested in the expectation between two variables $y_j$ and $y_j'$ of the resulting distribution, that is

$$E\{y_j y_j'\}$$

where $E\{\}$ is the expectation operator. Note:

  • The whole PDF of $y$ can be computed just by applying a change of variable. Unfortunately the integrals leading to the expectations are intractable (to the best of my knowledge).

  • I'm looking for closed-form approximations, no Markov Chain-Monte Carlo, no variational stuff. That is, approximations to the variable change, to the expectation integrals, to the sigmoid, to the resulting PDF, that allow computing $E\{y_j y_j'\}$.

  • Some dead ends: Taylor, often used in papers on the topic, is inaccurate by the mile. Gradshteyn and Ryzhik does not seem to contain the integrals.

  • $\begingroup$ What do you mean by "expectation between two variables"? $\endgroup$
    – bayerj
    Mar 2 '11 at 15:20
  • 1
    $\begingroup$ @user3509, you can get an easy lower bound by noting that $x^{-1}$ is convex and if $X$ is marginally Gaussian, then $e^{-X}$ is lognormal and so $\mathbb{E}(Y_i Y_j) > (1 + \mathbb{E} e^{-X_i} + \mathbb{E} e^{-X_j} + \mathbb{E} e^{-X_i-X_j})^{-1}$. The expectations on the right hand side are easy to compute. $\endgroup$
    – cardinal
    Mar 2 '11 at 15:53
  • $\begingroup$ Stationary phase (using a linear approximation to find the critical point) looks like it will work, but the resulting expressions--although they will be algebraic combinations of $\mathbf{M}$ and $\mathbf{S}$--will be extremely messy. How far are you prepared to go with this and what's the purpose of a "closed form" solution (when perhaps a numerical algorithm will in practice be more efficient and accurate)? $\endgroup$
    – whuber
    Mar 2 '11 at 16:40
  • $\begingroup$ @bayer, I mean E{y_1,y_2}=\int\int s(x_1)s(x_2) p(x_1,x_2) dx_1 dx_2 where y=s(x) is the sigmoid function, p(x,y) is the joint-PDF (gaussian). The integral is over RxR. $\endgroup$
    – user3509
    Mar 2 '11 at 17:05
  • 1
    $\begingroup$ @user3509 TeX is enabled in comments as well as questions and replies: surround it with $\$$, as in "$\$$E\{y_1,y_2\}=\int\int s(x_1)s(x_2) p(x_1,x_2) dx_1 dx_2$\$$" gives $E\{y_1,y_2\}=\int\int s(x_1)s(x_2) p(x_1,x_2) dx_1 dx_2$ $\endgroup$
    – whuber
    Mar 2 '11 at 17:12

The question really concerns pairs of normal variates. Let's call them $x_1$ and $x_2$ with means $\mu_i$, standard deviations $\sigma_i$, and correlation $\rho$. Whence their joint pdf is

$$\frac{1}{2 \pi \sqrt{1 - \rho^2} \sigma_1 \sigma_2} e^{-\frac{1}{1-\rho^2} \left(\frac{(x_1 - \mu_1)^2}{2 \sigma_1^2} + \frac{(x_2 - \mu_2)^2}{2 \sigma_2^2} - \frac{\rho (x_1 - \mu_1)(x_2 - \mu_2)}{\sigma_1 \sigma_2}\right)} dx_1 dx_2\text{.}$$

Let $f(x_1,x_2)$ be the product of this with the $y_i$ (as functions of the $x_i$). The first component of the gradient of $\log(f)$ is

$$\frac{\partial \log(f)}{\partial x_1} = \frac{1}{1 + e^{x_1}} + \frac{\rho(\mu_2 - x_2) \sigma_1 + (x_1 - \mu_1)\sigma_2}{(\rho^2-1)\sigma_1^2 \sigma_2},$$

with a similar expression for the second component (via the symmetry achieved by exchanging the subscripts 1 and 2). There will be a unique global maximum, which we can detect by setting the gradient to zero. This pair of nonlinear equations has no closed form solution. It is rapidly found by a few Newton-Raphson iterations. Alternatively, we can linearize these equations. Indeed, through second order, the first component equals

$$\frac{1}{2} + x_1\left(\frac{-1}{4} + \frac{1}{(\rho^2-1)\sigma_1^2}\right) + \frac{-\rho x_2 \sigma_1 + \rho \mu_2 \sigma_1 - \mu_1 \sigma_2}{(\rho^2 -1)\sigma_1^2 \sigma_2}.$$

This gives a pair of linear equations in $(x_1, x_2)$, which therefore do have a closed form solution, say $\hat{x}_i(\mu_1, \mu_2, \sigma_1, \sigma_2, \rho)$, which obviously are rational polynomials.

The Jacobian at this critical point has 1,1 coefficient

$$\frac{e^\hat{x_1}\left(2 - (\rho^2-1)\sigma_1^2 + 2\cosh(\hat{x_1})\right)}{(1+e^\hat{x_1})^2(\rho^2-1)\sigma_1^2},$$

1,2 and 2,1 coefficients

$$\frac{\rho}{\sigma_1 \sigma_2(1 - \rho^2)},$$

and 2,2 coefficient obtained from the 1,1 coefficient by symmetry. Because this is a critical point (at least approximately), we can substitute

$$e^\hat{x_1} = \frac{(\rho^2-1)\sigma_1^2 \sigma_2}{(\mu_2 - \hat{x_2})\rho \sigma_1 + (\hat{x_1} - \mu_1)\sigma_2} - 1$$

and use that also to compute $\cosh(\hat{x_1}) = \frac{e^\hat{x_1} - e^{-\hat{x_1}}}{2}$, with a similar manipulation for $e^\hat{x_2}$ and $\cosh(\hat{x_2})$. This enables evaluation of the Hessian (the determinant of the Jacobian) as a rational function of the parameters.

The rest is routine: the Hessian tells us how to approximate the integral as a binormal integral (a saddlepoint approximation). The answer equals $\frac{1}{2\pi}$ times a rational function of the five parameters: that's your closed form (for what it's worth!).

  • $\begingroup$ I am a bit concerned over the linearization but it is definitely worth trying. Do you know of some available literature on stationary phase aproximation of integrals?. $\endgroup$
    – user3509
    Mar 2 '11 at 19:20
  • $\begingroup$ I was concerned about the linearization, too, but (a) the series actually is linear through second order and (b) extensive plotting of the pdf for means between -3 and 3, sds from 0.1 to 3, and correlations from -1 to 1, indicate it should generally work well (the linearized critical point appears quite close to the global maximum). I learned this technique from an article in The American Statistician perhaps some 10-15 years ago, but I can't locate it at this moment. $\endgroup$
    – whuber
    Mar 2 '11 at 19:27
  • $\begingroup$ @user3509 I found the TAS reference and linked to it in the text: it goes by the name "saddlepoint approximation." I apologize for confusing it with stationary phase (they are essentially the same thing: one with a real variable, the other with an imaginary variable). $\endgroup$
    – whuber
    Mar 2 '11 at 19:42
  • $\begingroup$ Yes, that was exactly what I could not understand I was looking for the complex exponential somewhere. Thank you very much. $\endgroup$
    – user3509
    Mar 2 '11 at 20:13
  • 1
    $\begingroup$ +1 for the nice answer and the link to the article. Did not know that it is feasible to use such approximations. I would have resorted to numerical integration without even trying to get the closed form solution. $\endgroup$
    – mpiktas
    Mar 3 '11 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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