Why sigmoid function instead of anything else? Why is the de-facto standard sigmoid function, $\frac{1}{1+e^{-x}}$, so popular in (non-deep) neural-networks and logistic regression? 
Why don't we use many of the other derivable functions, with faster computation time or slower decay (so vanishing gradient occurs less). Few examples are on Wikipedia about sigmoid functions. One of my favorites with slow decay and fast calculation is $\frac{x}{1+|x|}$.
EDIT
The question is different to Comprehensive list of activation functions in neural networks with pros/cons as I'm only interested in the 'why' and only for the sigmoid.
 A: Since the original question mentioned the decaying gradient problem, I'd just like to add that, for intermediate layers (where you don't need to interpret activations as class probabilities or regression outputs), other nonlinearities are often preferred over sigmoidal functions. The most prominent are rectifier functions (as in ReLUs), which are linear over the positive domain and zero over the negative. One of their advantages is that they're less subject to the decaying gradient problem, because the derivative is constant over the positive domain. ReLUs have become popular to the point that sigmoids probably can't be called the de-facto standard anymore.

Glorot et al. (2011). Deep sparse rectifier neural networks

A: Quoting myself from this answer to a different question:

In section 4.2 of Pattern Recognition and Machine Learning (Springer 2006), Bishop shows that the logit arises naturally as the form of the posterior probability distribution in a Bayesian treatment of two-class classification. He then goes on to show that the same holds for discretely distributed features, as well as a subset of the family of exponential distributions. For multi-class classification the logit generalizes to the normalized exponential or softmax function.

This explains why this sigmoid is used in logistic regression.
Regarding neural networks, this blog post explains how different nonlinearities including the logit / softmax and the probit used in neural networks can be given a statistical interpretation and thereby a motivation. The underlying idea is that a multi-layered neural network can be regarded as a hierarchy of generalized linear models; according to this, activation functions are link functions, which in turn correspond to different distributional assumptions.
A: I have asked myself this question for months. The answers on CrossValidated and Quora all list nice properties of the logistic sigmoid function, but it all seems like we cleverly guessed this function. What I missed was the justification for choosing it. I finally found one in section 6.2.2.2 of the "Deep Learning" book by Bengio (2016). In my own words:
In short, we want the logarithm of the model's output to be suitable for gradient-based optimization of the log-likelihood of the training data.
Motivation


*

*We want a linear model, but we can't use $z = w^T x + b$ directly as $z \in (-\infty, +\infty)$. 

*For classification, it makes sense to assume the Bernoulli distribution and model its parameter $\theta$ in $P(Y=1) = \theta$.

*So, we need to map $z$ from $(-\infty, +\infty)$ to $[0, 1]$ to do classification.


Why the logistic sigmoid function?
Cutting off $z$ with $P(Y=1|z) = max\{0, min\{1, z\}\}$ yields a zero gradient for $z$ outside of $[0, 1]$. We need a strong gradient whenever the model's prediction is wrong, because we solve logistic regression with gradient descent. For logistic regression, there is no closed form solution.
The logistic function has the nice property of asymptoting a constant gradient when the model's prediction is wrong, given that we use Maximum Likelihood Estimation to fit the model. This is shown below:
For numerical benefits, Maximum Likelihood Estimation can be done by minimizing the negative log-likelihood of the training data. So, our cost function is:
$$
\begin{align}
J(w, b) &= \frac{1}{m} \sum_{i=1}^m -\log P(Y = y_i | x_i; w, b) \\ 
&= \frac{1}{m} \sum_{i=1}^m - \big(y_i \log P(Y=1 | z) + (y_i-1)\log P(Y=0 | z)\big)
\end{align}$$
Since $P(Y=0 | z) = 1-P(Y=1|z)$, we can focus on the $Y=1$ case. So, the question is how to model $P(Y=1 | z)$ given that we have $z = w^T x + b$.
The obvious requirements for the function $f$ mapping $z$ to $P(Y=1 | z)$ are:


*

*$\forall z \in \mathbb{R}: f(z) \in [0, 1]$

*$f(0) = 0.5$

*$f$ should be rotationally symmetrical w.r.t. $(0, 0.5)$, i.e. $f(-x) = 1-f(x)$, so that flipping the signs of the classes has no effect on the cost function.

*$f$ should be non-decreasing, continuous and differentiable.


These requirements are all fulfilled by rescaling sigmoid functions. Both $f(z) = \frac{1}{1 + e^{-z}}$ and $f(z) = 0.5 + 0.5 \frac{z}{1+|z|}$ fulfill them. However, sigmoid functions differ with respect to their behavior during gradient-based optimization of the log-likelihood. We can see the difference by plugging the logistic function $f(z) = \frac{1}{1 + e^{-z}}$ into our cost function.
Saturation for $Y=1$
For $P(Y=1|z) = \frac{1}{1 + e^{-z}}$ and $Y=1$, the cost of a single misclassified sample (i.e. $m=1$) is:
$$
\begin{align}
J(z) &= -\log(P(Y=1|z)) \\
&= -\log(\frac{1}{1 + e^{-z}}) \\
&= -\log(\frac{e^z}{1+e^z}) \\
&= -z + \log(1 + e^z)
\end{align}
$$
We can see that there is a linear component $-z$. Now, we can look at two cases:


*

*When $z$ is large, the model's prediction was correct, since $Y=1$. In the cost function, the $\log(1 + e^z)$ term asymptotes to $z$ for large $z$. Thus, it roughly cancels the $-z$ out leading to a roughly zero cost for this sample and a weak gradient. That makes sense, as the model is already predicting the correct class.

*When $z$ is small (but $|z|$ is large), the model's prediction was not correct, since $Y=1$. In the cost function, the $\log(1 + e^z)$ term asymptotes to $0$ for small $z$. Thus, the overall cost for this sample is roughly $-z$, meaning the gradient w.r.t. $z$ is roughly $-1$. This makes it easy for the model to correct its wrong prediction based on the constant gradient it receives. Even for very small $z$, there is no saturation going on, which would cause vanishing gradients.


Saturation for $Y=0$
Above, we focussed on the $Y=1$ case. For $Y=0$, the cost function behaves analogously, providing strong gradients only when the model's prediction is wrong.
This is the cost function $J(z)$ for $Y=1$:

It is the horizontally flipped softplus function. For $Y=0$, it is the softplus function.
Alternatives
You mentioned the alternatives to the logistic sigmoid function, for example $\frac{z}{1+|z|}$. Normalized to $[0,1]$, this would mean that we model $P(Y=1|z) = 0.5 + 0.5 \frac{z}{1+|z|}$.
During MLE, the cost function for $Y=1$ would then be 
$J(z) = - \log (0.5 + 0.5 \frac{z}{1+|z|})$,
which looks like this:

You can see, that the gradient of the cost function gets weaker and weaker for $z \rightarrow - \infty$.
A: One reason this function might seem more "natural" than others is that it happens to be the inverse of the canonical parameter of the Bernoulli distribution:
\begin{align}
f(y) &= p^y (1 - p)^{1 - y} \\
&= (1 - p) \exp \left \{ y \log \left ( \frac{p}{1 - p} \right ) \right \} .
\end{align}
(The function of $p$ within the exponent is called the canonical parameter.)
Maybe a more compelling justification comes from information theory, where the sigmoid function can be derived as a maximum entropy model.  Roughly speaking, the sigmoid function assumes minimal structure and reflects our general state of ignorance about the underlying model.
