A Gaussian scale mixture representation of the logistic distribution Can the logistic distribution with density function
$$f(x) = \frac{e^{-x}}{\left(1 + e^{-x}\right)^2}$$
be represented as a Gaussian scale mixture? In other words, if
\begin{align*}
X|V &\sim N(0, V) \\
V &\sim g,
\end{align*}
is there a choice for $V$ which leads to $X$ having a marginal logistic distribution?
 A: This paper by Leonard A. Stefanski (1991) illustrates that the answer to this question is yes, although the mixing density $g$ is difficult to work with in practice. Mixing over a Gamma density with a particular choice of parameters can also be shown to provide an extremely good approximation to the logistic distribution.
Exact Result
Stefanski (1991) demonstrates that the logistic distribution can be represented as a Gaussian scale mixture by setting $\sqrt V \sim q$, where $q(x) = \frac{d}{dx}L(x/2)$ and $L$ is the Kolmogorov-Smirnov cumulative distribution function. In the present notation, the density corresponding to $g$ becomes
$$g(v) = \sum_{n=1}^\infty (-1)^{n+1}n^2\exp\left(-\frac{n^2v}{2}\right)$$
This distribution can be difficult to work with in practice. This infinite series must be truncated after a finite number terms, and may be expensive to compute when $v$ is small (in fact, $g(v)$ will be negative for small $v$, unless you remember to choose an odd number for the truncation parameter). Working with $\log g(v)$ offers no improvement. Moreover, it is not easy to generate random variates from $g(v)$ and full conditionals will be intractable when using Gibbs sampling.
Approximation via Gamma
For the reasons described in the previous paragraph, it would be convenient to have an alternative choice for $g$ which leads to a marginal logistic distribution for $X$. The gamma distribution is easy to compute, easy to simulate from and often leads to tractable distributions during Gibbs sampling. Thus we consider the scale mixture representation.
\begin{align*}
X|V &\sim N(0, V) \\
V &\sim \text{Gamma}(\alpha, \beta),
\end{align*}
Note the following.
\begin{align*}
\mu = E(X) &= 0 \\[1.5ex]
\sigma^2 = Var(X) = E(Var(X|V)) + Var(E(X|V)) = E(V) &= \frac{\alpha}{\beta} \\[1.5ex]
\kappa = E\left[\left(\frac{X-0}{\sqrt{\alpha/\beta}}\right)^4\right] = \frac{\beta^2}{\alpha^2} E\left[V^2\left(\frac{X}{\sqrt V}\right)^4\right] = \frac{3\beta^2}{\alpha^2}E(V^2) &= 3\left(\frac{1}{\alpha}+1\right)
\end{align*}
Using a simple moment-matching approach, we set
$$3\left(\frac{1}{\alpha}+1\right) = \frac{21}{6} \quad\quad\quad \frac{\alpha}{\beta} = \frac{\pi^2}{3},$$
where $21/6$ and $\pi^2/3$ are the kurtosis and variance, respectively, of the logistic distribution. Solving for $\alpha$ and $\beta$ gives
$$\alpha = 2.5 \quad\quad \beta = \frac{7.5}{\pi^2} = 0.75991\ldots.$$
Marginal Distribution of $X$
The joint distribution of $X$ and $V$ can be written as
$$f(x, v) = cv^{\left(\alpha-\frac{1}{2}\right) - 1}\exp\left(-\frac{1}{2}\left(2\beta v + x^2\frac{1}{v}\right)\right), v > 0, x \in \mathbb R.$$
By noticing that this function is proportional (in $v$) to a generalized inverse Gaussian distribution, we can derive the marginal distribution of $X$ as having density
$$f(x) = \frac{\beta^{\frac{\alpha}{2} + \frac{1}{4}}}{2^{\frac{\alpha}{2} - \frac{3}{4}}\sqrt{\pi}\Gamma(\alpha)}|x|^{\alpha-1/2}K_{\alpha-1/2}\left(\sqrt{2\beta}|x|\right),$$
where $K_p(x)$ is the modified Bessel function of the second kind. In the special case where $\alpha=2.5$ and $\beta = 7.5/\pi^2$, this density reduces to
$$f(x) = \frac{\sqrt{3375}}{3\pi^4}x^2K_2\left(\frac{\sqrt{15}}{\pi}|x|\right)$$
Note that this is a special case of the Bessel function distribution.
Numerical Simulations
#Simulate data
set.seed(1234)
N <- 1e5
V <- rgamma(N, 2.5, 7.5/pi^2)
X <- rnorm(N, 0, sqrt(V))

#Plot samples
hist(X, freq=F, breaks=100, xlim=c(-8, 10), main="")

# Plot logistic density
curve(dlogis(x), add=T, lwd=3, col='red', lty=2)

# Plot derived density
curve(sqrt(3375)/(3*pi^4)*x^2*besselK(sqrt(15)/pi*abs(x), 2), add=T, col='blue', lwd=3, lty=3)
legend('topright', c("Logistic", "Gamma-Mixture"), lwd=3, lty=c(2,3), col=c('red', 'blue'), bty='n', cex=1.2)

#Make QQ Plot
q_logis <- qlogis(seq(0, 1, length.out=N+2)[-c(1,N+2)])
qqplot(X, q_logis,
   xlab="Sample Quantiles", ylab="Theoretical Quantiles")
abline(0, 1, lwd=2, col='red', lty=3)



