8
$\begingroup$

I have a random variable $Y = \frac{e^{X}}{1 + e^{X}}$ and I know $X \sim N(\mu, \sigma^2)$.

Is there a way to compute $\mathbb{E}(Y)$? I have tried to work out the integral, but haven't made much progress. Is it even possible?

$\endgroup$
  • 3
    $\begingroup$ Apparently, an analytical solution is not known. A known approximation is given in this Math stackexchange link: math.stackexchange.com/questions/207861/… $\endgroup$ – Greenparker Nov 26 '18 at 15:07
  • $\begingroup$ If $\mu = 0$, then $E[Y] = \frac12$, for any $\sigma$. $\endgroup$ – wolfies Nov 26 '18 at 15:24
  • $\begingroup$ @wolfies Could you give a source/derivation of that? $\endgroup$ – Greenparker Nov 26 '18 at 15:26
  • 2
    $\begingroup$ @Greenparker The distribution of $Y-1/2$ is symmetric around $0$ in that case, QED. $\endgroup$ – whuber Nov 26 '18 at 16:11
  • 1
    $\begingroup$ I did it symbolically as a one-liner with mathStatica/Mathematica ... but an easy way to see why it must be so: ...... (i) if $X \sim N(0,\sigma^2)$, then its pdf is symmetric around 0. (ii) Consider the transformation $Z = Y-\frac12 = \frac12 \tanh(x/2)$. Then $Z$ is a symmetric S-shaped curve about $X = 0$, and E[Z] must equal 0 (by symmetry). Since $Y = Z + \frac12$, it follows that $E[Y] = \frac12$ $\endgroup$ – wolfies Nov 26 '18 at 16:11
5
+100
$\begingroup$

As mentioned already in the question comments and answer by @Martijn there doesn't appear to be an analytical solution for $E(Y)$ apart from the special case where $\mu = 0$ which gives $E(Y) = 0.5$.

Additionally by Jensen's inequality we have that $E(Y) = E(f(X)) < f(E(X))$ if $\mu > 0$ and conversely that $E(Y) = E(f(X)) > f(E(X))$ if $\mu < 0$. Since $f(x) = \frac{e^x}{1 + e^x}$ is convex when $x < 0$ and concave when $x > 0$ and most of the normal density mass will lie in those regions depending on the value of $\mu$.

There are many ways to approximate $E(Y)$, I have detailed a few I am familiar with and included some R code at the end.

Sampling

This is quite easy to understand/implement:

$$ E(Y) = \int_\infty^\infty f(x) \mathcal{N}(x|\mu, \sigma^2) dx \approx \frac{1}{n} \Sigma_{i = 1}^{n}f(x_i) $$

where we draw samples $x_1, \ldots, x_n$ from $\mathcal{N}(\mu, \sigma^2)$.

Numerical integration

This includes many methods of approximating the integral above - in the code I used R's integrate function which uses adaptive quadrature.

Unscented transform

See for example The Unscented Kalman Filter for Nonlinear Estimation by Eric A. Wan and Rudolph van der Merwe which describes:

The unscented transformation (UT) is a method for calculating the statistics of a random variable which undergoes a nonlinear transformation

The method involves calculating a small number of "sigma points" which are then transformed by $f$ and a weighted mean is taken. This is in contrast to randomly sampling many points, transforming them with $f$ and taking the mean.

This method is much more computationally efficient than random sampling. Unfortunately I couldn't find an R implementation online so haven't included it in the code below.

Code

The following code creates data with different values of $\mu$ and fixed $\sigma$. It outputs f_mu which is $f(E(X))$, and approximations of $E(Y) = E(f(X))$ via sampling and integration.

integrate_approx <- function(mu, sigma) {
    f <- function(x) {
        plogis(x) * dnorm(x, mu, sigma)
    }
    int <- integrate(f, lower = -Inf, upper = Inf)
    int$value
}

sampling_approx <- function(mu, sigma, n = 1e6) {
    x <- rnorm(n, mu, sigma)
    mean(plogis(x))
}

mu <- seq(-2.0, 2.0, by = 0.5)

data <- data.frame(mu = mu,
                   sigma = 3.14,
                   f_mu = plogis(mu),
                   sampling = NA,
                   integration = NA)

for (i in seq_len(nrow(data))) {
    mu <- data$mu[i]
sigma <- data$sigma[i]
    data$sampling[i] <- sampling_approx(mu, sigma)
data$integration[i] <- integrate_approx(mu, sigma)
}

output:

    mu sigma      f_mu  sampling integration
1 -2.0  3.14 0.1192029 0.2891102   0.2892540
2 -1.5  3.14 0.1824255 0.3382486   0.3384099
3 -1.0  3.14 0.2689414 0.3902008   0.3905315
4 -0.5  3.14 0.3775407 0.4450018   0.4447307
5  0.0  3.14 0.5000000 0.4999657   0.5000000
6  0.5  3.14 0.6224593 0.5553955   0.5552693
7  1.0  3.14 0.7310586 0.6088106   0.6094685
8  1.5  3.14 0.8175745 0.6613919   0.6615901
9  2.0  3.14 0.8807971 0.7105594   0.7107460

EDIT

I actually found an easy to use unscented transformation in the python package filterpy (although it is actually quite quick to implement from scratch):

import filterpy.kalman as fp
import numpy as np
import pandas as pd


def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))


m = 9
n = 1
z = 1_000_000
alpha = 1e-3
beta = 2.0
kappa = 0.0
means = np.linspace(-2.0, 2.0, m)
sigma = 3.14
points = fp.MerweScaledSigmaPoints(n, alpha, beta, kappa)
ut = np.empty_like(means)
sampling = np.empty_like(means)

for i, mean in enumerate(means):
    sigmas = points.sigma_points(mean, sigma**2)
    trans_sigmas = sigmoid(sigmas)
    ut[i], _ = fp.unscented_transform(trans_sigmas, points.Wm, points.Wc)

    x = np.random.normal(mean, sigma, z)
    sampling[i] = np.mean(sigmoid(x))

print(pd.DataFrame({"mu": means,
                    "sigma": sigma,
                    "ut": ut,
                    "sampling": sampling}))

which outputs:

    mu  sigma        ut  sampling
0 -2.0   3.14  0.513402  0.288771
1 -1.5   3.14  0.649426  0.338220
2 -1.0   3.14  0.716851  0.390582
3 -0.5   3.14  0.661284  0.444856
4  0.0   3.14  0.500000  0.500382
5  0.5   3.14  0.338716  0.555246
6  1.0   3.14  0.283149  0.609282
7  1.5   3.14  0.350574  0.662106
8  2.0   3.14  0.486598  0.710284

So the unscented transform seems to perform quite poorly for these values of $\mu$ and $\sigma$. This is maybe not surprising since the unscented transform attempts to find the best normal approximation to $Y = f(X)$ and in this case it is far from normal:

import matplotlib.pyplot as plt

x = np.random.normal(means[0], sigma, z)
plt.hist(sigmoid(x), bins=50)
plt.title("mu = {}, sigma = {}".format(means[0], sigma))
plt.xlabel("f(x)")
plt.show()

histogram

For smaller values of $\sigma$ it seems to be OK.

$\endgroup$
3
$\begingroup$

The variable $Y$ has a logit normal or logistic normal distribution whose moments have no known analytical description. You can obtain the values computationally.

More about these distributions is described in a freely available article: Atchison, J., and Sheng M. Shen. "Logistic-normal distributions: Some properties and uses." Biometrika 67.2 (1980): 261-272.

In that text they do not give any expressions for limits, approximations or behavior of the moments (except mentioning that they exist). But, they do continue with expressions for the expectation value for the ratio of two components in a multivariate logistic normal distributed variable.

$\endgroup$

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.