# How to calculate or approximate the integral $\int_{-\infty}^{\infty}\sigma(x)(1-\sigma(x))\mathcal{N}(x|0,1)dx$

$\sigma(x)$ is the sigmoid function, that is, $\sigma(x)=\frac{1}{1+e^{-x}}$. And $x$ is from a Normal distribution of 0 mean and 1 variance. Now I'd like to calculate the expectation of the derivative value of the sigmoid function, i.e,

$\int_{-\infty}^{\infty}\sigma(x)(1-\sigma(x))\mathcal{N}(x|0,1)dx$

• Monte Carlo will work. You can generate variables from the standard normal distribution, plug them into $\sigma(x) \left(1-\sigma(x)\right)$ and take the average. Commented Jun 7, 2017 at 7:26
• People are suggesting Monte Carlo which will work but is generally a terrible idea for a $1D$ integral. Use numerical integration or (adaptive) quadrature, it is much better/faster for low dimensions (up to $D=3$, if I remember correctly; from $D=4$, the scaling of the error of numerical integration and Monte Carlo becomes the same). Commented Jun 7, 2017 at 8:47
• @lacerbi Thanks for your comment! Indeed I also prefer numerical approach. Do you have any good idea? Thanks! Commented Jun 7, 2017 at 14:23
• Here I would just do a brutal integration via Simpson's rule: en.wikipedia.org/wiki/Simpson%27s_rule As for the bounds, go as far as you want in terms of SDs of the normal integral (presumably five or six SDs are more than enough). I would not recommend Gauss-Hermite quadrature here because the other piece of the integrand (besides the normal) is probably not well approximated by a polynomial. Commented Jun 7, 2017 at 19:46

Mathematica doesn't give an analytic result, but it does give the numerical result 0.206621.

If you want some sort of semi-analytic treatment, you can expand $$\sigma(x) ( 1 - \sigma(x) ) = \frac{1}{4} - \frac{x^2}{16} + \frac{x^4}{96} + \cdots$$ and then integrate term by term to get $$\frac{1}{4} - \frac{1}{16} + \frac{1}{32} + \cdots$$ The first three terms get you to 0.22. If you can write the $n$th term of the expansion in closed form, you might even be able to recognize the term-by-term integrated series as related to some established function or constant.

• The power series approach is mathematically invalid, because the radius of convergence of this power series for $\sigma$ is just $\pi$. I believe your method works approximately, and that only because most of the Normal probability is concentrated within the radius of convergence and $\sigma$ is bounded. In short, it looks like you got lucky. With more careful analysis, you could estimate the error and justify this approach as an approximate calculation.
– whuber
Commented Jun 7, 2017 at 13:24
• @whuber: It's a fair point, but it's more that just luck -- it's an asymptotic series. More exactly, if you introduce a width parameter for the Gaussian $s$, it's an asymptotic series in $s$. For $s=1$, the terms get down to $\sim 0.02$ before they start growing, so that's the accuracy limit. And the fact that it's asymptotic doesn't prevent it being Borel summable, so the route to a fully closed-form solution isn't closed by the divergence. Commented Jun 7, 2017 at 16:34
• I am unable to make much sense of those remarks, David, for two reasons. First, I don't see any connection between your "$s$" and the function $\sigma(x)(1-\sigma(x))$. Second, I understand "asymptotic" to mean that when you truncate the series at a fixed number of terms, then as $x$ grows large the truncated series grows arbitrarily close to the true value, relative to the true value. But the opposite holds: no matter where you truncate it, it will diverge for large $x$.
– whuber
Commented Jun 8, 2017 at 19:17

Expanding on JohnK's comment: You want to find the expectation of a function with respect to a standard normal distribution. Since it seems difficult to evaluate the integral, we can instead resort to estimating it using Monte Carlo.

If we sample $x_1, x_2, \dots x_N$ from a $N(0,1)$ distribution, then your quantity can be estimated by $$\dfrac{1}{n} \sum_{i=1}^{N} \sigma(x_i) (1 - \sigma(x_i)\,.$$

As you increase $N$, you can an increasingly better estimate of the quantity. The R code below returns the value .20669 with standard error .00015.

> set.seed(100)
> sig <- function(x)
+ {
+   1/(1+exp(-x))
+ }
>
> N <- 1e5
> samples <- rnorm(N, 0, 1)
> mean(sig(samples)*(1 - sig(samples)))
[1] 0.2066863

> # The standard error of the mean estimate
> sd(sig(samples)*(1 - sig(samples)))/sqrt(N)
[1] 0.0001463474


In addition, note that $$\sigma(x)\{1-\sigma(x)\} = \frac{e^{-x}}{(1+e^{-x})^2} = \left(\frac{1}{1+e^{-x}}\right)'$$ Thus, integrating by part, $$\int \sigma(x)\{1-\sigma(x)\}N(0,1)dx = \left[\frac{1}{1+e^{-x}} N(0,1)(x)\right] - \int \frac{1}{1+e^{-x}} N(0,1)' dx$$ The $[...]$ term equals 0. Then $$N(0,1)' = - x N(0,1)$$ Finally, if my calculations are correct, your integral is $$\int \frac{x}{1+e^{-x}}N(0,1)dx$$ Applying the Monte Carlo technique, generate $x_1,\ldots,x_n \sim N(0,1)$ and compute then average $$\frac{1}{n} \sum_{i=1}^n \frac{x_i}{1+e^{-x_i}}$$ It might be more stable [not even sure], though pretty much equivalent to JohnK and GreenPaker's suggestions.

• Given that the Monte Carlo experiment is a few lines of code, could you run the comparison in term of (estimated) standard error for both Monte Carlo approaches? Commented Jun 7, 2017 at 13:46
• The integration by parts is a nice idea. It doesn't seem to help though: it replaces an exponentially decreasing function by one that linearly increases. That likely makes the integral a little bit harder to evaluate numerically.
– whuber
Commented Jun 12, 2017 at 13:07