2 Existing multivariate logistic distributions

There are several existing multivariate distributions linked in some way to the univariate logistic distribution. One class of such distributions has marginal distributions that are logistic, but lacks rotational symmetry, including the bivariate case in [Reference Gumbel4] and more general forms in [Reference Ali, Mikhail and Haq1, Reference Malik and Abraham8], plus the special case of independence.

Circular symmetry of a multivariate distribution, as well as rotational symmetry of a whole family of distributions, can be achieved by considering elliptically symmetric distributions [Reference McGraw and Wagner9]. Typically, such distributions are written in the form

$$ \begin{align*} f_n({\mathbf{x}}) = |\Sigma|^{1/2}g_n(({\mathbf{x}}-\mu)^T\Sigma^{-1}({\mathbf{x}}-\mu)), \end{align*} $$

where $\Sigma $ is a positive definite $n\times n$ matrix, ${\mathbf {x}},\mu \in \mathbb {R}^n$ and $g_n(\cdot )$ is known as the kernel or density generator; often, $g_n$ depends on the dimension n only through a normalising constant $c_n$ [Reference Kelker6]. There are then two existing ways in which a link has been made with the logistic distribution.

The first, giving what is widely referred to as the elliptically symmetric logistic distribution, involves taking

$$ \begin{align*} g_n(u) = \frac{c_n\exp(-u)}{(1+\exp(-u))^2}, \end{align*} $$

that is, using a density generator based on the standard logistic distribution [Reference Jensen, Kotz, Johnson and Read5, Reference Yin, Wang and Sha17]. However, this approach does not have the property that the score function is bounded, even when $n=1$ , as this distribution does not reduce to the usual logistic in that case. See Section 5 for details.

The second approach is to require that $f_n({\mathbf{x}})$ is elliptically symmetric with marginal distributions given by the univariate logistic distribution. Volodin [Reference Volodin15] discusses this case (see references therein for its origins) and gives a general construction for elliptically symmetric distributions with marginals $f(\cdot )$ , taking very different forms in odd and even dimensions. For even dimension $n=2m$ , the density involves m-fold differentiation with respect to $r = \sqrt {({\mathbf {x}}-\mu )^T\Sigma ^{-1}({\mathbf {x}}-\mu )}$ ) of $\int _0^\infty f(\sqrt {z^2+r})\,dz$ [Reference Volodin15, (2)]. In the logistic case, Volodin [Reference Volodin15] shows that the integral can be rewritten as an infinite series. For odd dimension $n=2m+1$ , the density is somewhat simpler, but still involves m-fold differentiation of $f(\cdot )$ . Though this distribution has logistic marginals and elliptical symmetry, it lacks the required tractability, particularly when $n=2$ .

3 The proposed distribution

Staying with the elliptically symmetric framework, instead of taking the density generator to be the univariate logistic, we can simply take the square root first, to give

(3.1) $$ \begin{align} g_n(u) = \frac{c_n\exp(-\sqrt{u})}{(1+\exp(-\sqrt{u}))^2}. \end{align} $$

Then, we have

(3.2) $$ \begin{align} f_n({\mathbf{x}} \mid \mu,\Sigma) ={|\Sigma|^{-1/2}}c_n\frac{\exp(-r({\mathbf{x}}))}{(1+\exp(-r({\mathbf{x}})))^2}, \end{align} $$

where

(3.3) $$ \begin{align} r({\mathbf{x}}) = \sqrt{({\mathbf{x}}-\mu)^T\Sigma^{-1}({\mathbf{x}}-\mu)}. \end{align} $$

When $n=1$ , this reduces to the usual univariate logistic distribution, with $c_1=1$ . In two dimensions, it is straightforward to get the normalising constant $c_2$ . Since the Jacobian for transforming to polar co-ordinates $(r,\theta )$ is just r, we need

$$ \begin{align*} c^{-1}_2 &= \int_{\theta=-\pi}^{\pi}\int_{r=0}^\infty r f_n(r,\theta) \,d r \,d \theta\\ &= 2\pi \int_{r=0}^\infty r f_n(r) \,d r\\ &= 2\pi \int_{r=0}^\infty \frac{r \exp(-r)}{(1 + \exp(-r))^2} \,dr\\ &= 2\pi \log(2). \end{align*} $$

Explicitly,

$$ \begin{align*} f_2({\mathbf{x}} \mid \mu,\Sigma) = \frac{|\Sigma|^{-1/2}}{2\pi\log(2)}\frac{\exp(-r({\mathbf{x}}))}{(1+\exp(-r({\mathbf{x}})))^2}. \end{align*} $$

Normalising constants for $n\geq 3$ are given in Section 4.

Clearly, this proposed distribution is both simple in form and elliptically symmetric. It occurs as a particular case of the very general family defined by Wang and Yin [Reference Wang and Yin16]; the next section makes that link explicitly and builds on results of [Reference Wang and Yin16] to obtain the normalising constant $c_n$ in general. In later sections, I will show that its score function is bounded (Section 5) and explore its marginal distributions (Section 6).

4 Normalising constants

Wang and Yin [Reference Wang and Yin16] define a very general family of multivariate densities that includes the proposed density of (3.2) as a special case. They consider density generators of the form

$$ \begin{align*} g_n(u) = \frac{c_n u^{N-1}\exp(-au^{s_1})}{(1+\exp(-bu^{s_2})^{2r}}, \end{align*} $$

(their parameter r is not to be confused with the function $r(\cdot )$ used here for (scaled) distance from the origin) which includes the proposed distribution when $N=1, s_1=s_2=1/2, f=1,a=b=1, r=1$ . Wang and Yin [Reference Wang and Yin16] do not give much information about the particular distribution proposed in (3.1); this case is not included in any of the eight sub-families that they explore in more detail. However, their Theorem 2.1 relates the normalising constants $c_n$ to the generalised Hurwitz–Lerch zeta function $\Phi ^*$ , defined in [Reference Wang and Yin16, Section 2.2], in cases where $s_1=s_2$ . For the case proposed here, the expression for $c_n$ simplifies in one of two ways. For $n\geq 4$ , the sufficient conditions to express $\Phi ^*$ in summation form are met. For $n=2,3$ , an integral expression for $\Phi ^*$ is given, but direct calculation shows that the summation form still holds even though Wang and Yin’s sufficient conditions are not met. For all n, the necessary value of $\Phi ^*$ simplifies to $\eta (n-1)$ , where

$$ \begin{align*} \eta(t) = \sum_{i=1}^\infty\frac{(-1)^{i-1}}{i^t} \end{align*} $$

is the Dirichlet eta function or alternating zeta function, and

(4.1) $$ \begin{align} c_n = \frac{\Gamma(n/2)}{2\pi^{n/2}\Gamma(n)\eta(n-1)}. \end{align} $$

5 Score functions

The linear score function,

$$ \begin{align*} \nabla\log(f({\mathbf{x}})) = \bigg(\frac{\partial\log f}{\partial x_1},\ldots,\frac{\partial\log f}{\partial x_n}\bigg)^T, \end{align*} $$

is of particular interest because of its role in constructing models in movement ecology that are tailored to specific stationary distributions, which can be thought of as analogous to the construction of Markov Chain Monte Carlo algorithms for specific target distributions [Reference Michelot, Blackwell, Chamaillé-Jammes and Matthiopoulos10–Reference Michelot, Gloaguen, Blackwell and Etienne12]. For example, the Langevin algorithm for a target $\pi ({\mathbf {x}})$ has a drift term proportional to $-\nabla \log (\pi ({\mathbf {x}}))$ (see [Reference Michelot, Gloaguen, Blackwell and Etienne12]), and the Bouncy Particle Sampler [Reference Bouchard-Côté, Vollmer and Doucet3] involves ‘bounce’ events with a rate that depends on $\nabla \log (\pi ({\mathbf {x}}))$ . Boundedness of the score function is relevant to both the realism and the practical implementation of such models.

For a general elliptically symmetric distribution with density generator $g(\cdot )$ , and defining

$$ \begin{align*} q(\cdot) : q({\mathbf{x}}) &= ({\mathbf{x}}-\mu)^T\Sigma^{-1}({\mathbf{x}}-\mu),\\ h_n(\cdot) : h_n(u) &= \frac{d \log(g_n(u))}{du}, \end{align*} $$

we have

(5.1) $$ \begin{align} \nabla\log(f_n({\mathbf{x}})) = 2h_n(q({\mathbf{x}}))\Sigma^{-1}({\mathbf{x}}-\mu). \end{align} $$

For the existing elliptically symmetric multivariate logistic, with $g_n(\cdot )$ proportional to the univariate logistic,

$$ \begin{align*} h_n(u) = -\frac{1-\exp(-u)}{1+\exp(-u)} \end{align*} $$

and so $h_n(u) \rightarrow -1$ as $u \rightarrow \infty $ and $h_n(u) \rightarrow 1$ as $u \rightarrow -\infty $ . Substituting into (5.1), clearly the score is not bounded.

However, if we take $g_n(\cdot )$ as in (3.1), the proposed ‘square root’ version, we have

$$ \begin{align*} h_n(u) = -\frac{1}{2\sqrt{u}}\cdot\frac{1-\exp(-\sqrt{u})}{1+\exp(-\sqrt{u})}. \end{align*} $$

To see that the score is bounded in this case, it is simplest to consider the case $\Sigma =I,\mu ={\mathbf {0}}$ . Then,

$$ \begin{align*} \nabla\log(f_n({\mathbf{x}})) &= \frac{-{\mathbf{x}}}{\sqrt{{\mathbf{x}}^T{\mathbf{x}}}}\cdot\frac{1-\exp(-\sqrt{{\mathbf{x}}^T{\mathbf{x}}})}{1+\exp(-\sqrt{{\mathbf{x}}^T{\mathbf{x}}})},\\ |\nabla\log(f_n({\mathbf{x}}))| &= \frac{1-\exp(-\sqrt{{\mathbf{x}}^T{\mathbf{x}}})}{1+\exp(-\sqrt{{\mathbf{x}}^T{\mathbf{x}}})}. \end{align*} $$

So the score is bounded in this case and, by linear transformation, also for general $\Sigma $ and $\mu $ .

6 Marginal distributions

The marginal distributions for the distribution proposed in Section 3 can be evaluated by considering first the standard circular case $\Sigma =I_n,\mu ={\mathbf {0}}$ . There, the marginal (for any $x_i$ ) is not exactly a logistic distribution, except when $n=1$ ; empirically, it is close to one, although not with scale 1. Before making this relationship more precise, it is useful to assess the scale by evaluating the marginal at 0, which we can do exactly.

For any standard elliptically symmetric distribution where the density generator is of the form

$$ \begin{align*} g_n(u) = c_n g_0(u), \end{align*} $$

we have

$$ \begin{align*} f_{x_n}(x \mid \mu={\mathbf{0}},\Sigma=I_n) = \int_{{\mathbf{x}}_{(n)}} g_n(({\mathbf{x}}_{(n)},x_n)^T({\mathbf{x}}_{(n)},x_n)) \,d\mathbf{x}_{(n)}, \end{align*} $$

where, with a slight abuse of notation, we write $f_{x_n}(\cdot )$ for the univariate marginal for $x_n$ derived from $f_{n}(\cdot )$ , and ${{\mathbf {x}}_{(n)}}$ for $(x_1,\ldots ,x_{n-1})$ . Hence,

$$ \begin{align*} f_{x_n}(0 \mid \mu={\mathbf{0}},\Sigma=I_n) &= \int_{{\mathbf{x}}_{(n)}} g_n({\mathbf{x}}_{(n)}^T{\mathbf{x}}_{(n)}) \,d\mathbf{x}_{(n)}\\&= \int_{{\mathbf{x}}_{(n)}} c_n g_0({\mathbf{x}}_{(n)}^T{\mathbf{x}}_{(n)}) \,d\mathbf{x}_{(n)}\\&= \frac{c_n}{c_{n-1}}\int_{{\mathbf{x}}_{(n)}} c_{n-1} g_0({\mathbf{x}}_{(n)}^T{\mathbf{x}}_{(n)}) \,d\mathbf{x}_{(n)}\\&= \frac{c_n}{c_{n-1}}. \end{align*} $$

A univariate logistic with scale $\sigma $ has a density at zero of $1/4\sigma $ ; re-arranging means that the marginal density at zero matches that of a univariate logistic with scale

$$ \begin{align*} \tau_n = \frac{c_{n-1}}{4c_n}. \end{align*} $$

In particular, for $g_n(\cdot )$ given by (3.1), we have $c_1=1, c_2=(2\pi \log (2))^{-1}, c_3=(3/2)\pi ^{-3}$ , and so $\tau _2 = \pi \log (2)/2 \approx 1.088793, \tau _3 = \pi ^2/12\log (2) \approx 1.186569$ .

In both these cases, calculation of the marginal distribution by numerical integration (using the standard R routine [14]) shows that the main parts of the marginal densities are well approximated by logistic distributions with scale parameters $\tau _2$ and $\tau _3$ , respectively. However, in the extreme tails, there is some divergence in shape, with densities decaying more quickly than $\exp (-\tau _n x)$ but less quickly than $\exp (-x)$ .

For $n=2$ and $x_2\in [-6.5,6.5]$ , containing more than the central 99% of mass, the marginal density $f_{x_2}(x \mid \mu ={{\mathbf {0}}},\Sigma =I_2)$ and the univariate logistic density with scale parameter $\tau _2$ agree to within a factor of 1.01. A more precise expression for the density in the extreme tails is suggested by comparison with the standard bivariate Laplace distribution [Reference Kotz, Kozubowski and Podgorski7], which has joint density

$$ \begin{align*} \pi^{-1} K_0\Big(\sqrt{2(x_1^2+x_2^2)}\Big), \end{align*} $$

where $K_0(\cdot )$ is the modified Bessel function of the second kind, and marginal density

$$ \begin{align*} 2^{-1/2}\exp(-x\sqrt{2}). \end{align*} $$

The asymptotic form

$$ \begin{align*} K_0(z) \sim \sqrt{\pi/(2z)}\exp(-z) \end{align*} $$

(see, for example, [13]) suggests

$$ \begin{align*} \sqrt{x/\pi}\exp(-x) \end{align*} $$

for the marginal of the proposed logistic with $n=2$ , and numerical experiments show that it is indeed extremely close for sufficiently large x.

7 Implementation

Calculation of the proposed density is straightforward once the normalising constants are known. In R [14], the Dirichlet eta function is available in the pracma package [Reference Borchers2], and so calculation of $c_n$ from (4.1) is trivial. Some improvement in speed can be achieved by taking advantage of the existing function for the univariate logistic, $f_1(\cdot )$ , writing

$$ \begin{align*} f_n({\mathbf{x}} \mid \mu,\Sigma) = {|\Sigma|^{-1/2}} c_n f_1(r({\mathbf{x}})\mid 0,1), \end{align*} $$

where $r({\mathbf {x}})$ is as defined in (3.3).