©This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.
© 2024 The Author(s) AIChE Journal published by Wiley Periodicals LLC on behalf of American Institute of Chemical Engineers.

1 Introduction↩︎

In Statistics, the logistic transformation is one of the most common used, foremost in the context of logistic regression for binary, binomial and multinomial data, but also implicitly in Bayesian versions of negative binomial regression. In this paper, we focus on problems where the inverse logit function has been applied to a normal random variable \(X\), with realisation \(x \sim \mathcal{N}(\mu,\sigma^2)\) and hence the random variable \(P := (1+e^{-X})^{-1}\) is distributed logit-normal with parameters \(\mu, \sigma^2\). In frequentist inference, such random variables, including powers of, appear in logistic random intercept models when the link function \(\boldsymbol{\theta}_i\) in cluster \(i\) has the form \(\boldsymbol{\theta}_i ={\boldsymbol{X}}_i\boldsymbol{\beta} + {\boldsymbol{u}}_i, {\boldsymbol{u}} \sim \mathcal{N}(0,\sigma^2_u{\boldsymbol{I}})\). In this setting, logit-normal random variables appear in the marginal likelihood of \(\boldsymbol{\beta}\) [1], [2]. Additionally, the (marginal) response probability in a frequentist treatment of any logistic mixed model is a mean of a logit-normal random variable. In Bayesian inference, logit-normal random variables, including powers of, appear implicitly in Variational Bayes [3]–[5] and Expectation Propagation [6], [7] algorithms used to determine parameters of approximate posteriors in logistic and negative binomial regression.

Despite the need to find a collection of logit-normal moments in many inferential problems, such moments are poorly understood, and assumed analytically unknown. As a result, literature on the logit-normal distribution that discusses moments splits into disconnected strands. One strand gives up on attempting to find analytic formulae for any logit-normal moments, instead promoting the use of standard numerical integration techniques [8], [9]. The second, usually buried in textbooks, are approximations of the first moment based on similarity of the logit and probit link or between the normal and logistic distributions [10], [11]. The last, pioneered in [12], [13] and updated in [14], starts from the expected value of a logit-normal random variable, \(E(P)\), being a special case of the Mordell integral [15] from complex analysis, but restricted to the real numbers. While the Mordell integral based representation of \(E(P)\) given in [14] can be combined with other properties of the logit-normal distribution to determine higher order moments, in practice Mordell integral based estimates of \(E(P)\), and thus higher moments, become numerically unstable when \(|\mu/\sigma^2| \rightarrow \infty\). Moreover, the repeated application of the quotient rule from Calculus needed to find higher order moments from the Mordell integral representation of \(E(P)\) rapidly becomes algebraically very complicated.

In this paper, we change tack from [14], and no longer attempt to estimate moments starting from the Mordell integral representation of \(E(P)\). Instead, we focus on constructing highly accurate approximations of the logistic function and its derivatives using an basis of exponential functions. This allows us to not just determine an approximation for \(E(P)\) but also easily develop a proposition to find the \(k^\text{th}\) non-central moment, \(k \in \mathbb{Z}^+\). We will show our approximations, unlike Mordell integral representations, are numerically stable when \(|\mu/\sigma^2| \rightarrow \infty\). We will show our method is faster to compute than standard numerical integration in the programming language R . However, we will also show the number of terms needed for accurate approximation increases with moment order. Unfortunately, in practice the number of moments our method can calculate has an upper bound around 8. This limits the utility of our method for evaluation of the logistic-normal integral. Nevertheless, getting to the \(8^\text{th}\) moment is more than sufficiently high, which we will illustrate with an example, for performing expectation propagation in logistic regression.

2 Inferential problems that require logit-normal moments↩︎

To justify our interest in the logit-normal distribution, we will provide two illustrative inferential problems, one frequentist, the other Bayesian, where a collection of logit-normal moments are needed in order to estimate parameters. First though, we need to define the logit-normal distribution and some of its properties.

2.1 The logit-normal distribution and some of its properties↩︎

A logit-normal distributed random variable \(P\) with realisation \(p \sim \text{logit}N(\mu,\sigma^2)\) has probability density (pdf), \[f(p|\mu,\sigma^2) = \frac{1}{\sigma p(1-p)\sqrt{2\pi}}e^{-\frac{(\text{logit}(p)-\mu)^2}{2\sigma^2}}, \quad p \in (0,1), \text{logit}(p) = \log\bigg(\frac{p}{1-p}\bigg) \label{eq146o}\tag{1}\]

To determine the moments \(k = 1, 2, \ldots\), of \(P\) we will work with the following representations of the \(k^\text{th}\) moment, \[\begin{align} E(P^k) &=& \frac{1}{\sigma\sqrt{2\pi}}\int_{-\infty}^{\infty} (1+e^{-x})^{-r}e^{-\frac{z^2}{2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}} {\rm d} x = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} (1+e^{-\mu-\sigma z})^{-r}e^{-\frac{z^2}{2}} {\rm d}z, \label{eq4} \end{align}\tag{2}\] where \(x \sim N(\mu,\sigma^2), z \sim N(0,1)\). However in practice, to find higher order moments of \(P\), \(E(P^2), E(P^3), \ldots\), it is easier to use the recurrence relation given in [14]. This is found by taking the derivative of \(E(P^k)\) with respect to \(\mu\), where \(E(P^k)\) is in the form of (2 ), \[\begin{align} \frac{{\rm d}E(P^k)}{{\rm d}\mu} &=& \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty \frac{ {\rm d} (1+e^{-\mu - \sigma z})^{-k} }{{\rm d}\mu} e^{-\frac{z^2}{2}} {\rm d}z \nonumber \\ &=& \frac{k}{\sqrt{2\pi}}\int_{-\infty}^\infty \{(1+e^{-\mu - \sigma z})^{-k}- (1+e^{-\mu - \sigma z})^{-(k+1)}\}e^{-\frac{z^2}{2}} {\rm d}z \tag{3} \nonumber \\ &=& k\{E(P^k) - E(P^{k+1})\}, \tag{4} \\ E(P^{k+1}) &=& E(P^k) - \frac{1}{k} \frac{{\rm d}E(P^k)}{{\rm d}\mu}. \tag{5} \end{align}\] Thus the main challenge for determining logit-normal moments is evaluating the first moment or expected value in a form such that applying the recurrence relation in (5 ) is straightforward for any choice of \(k\).

2.2 A frequentist problem: Some cases of the logistic mixed model↩︎

In order to perform inference in frequentist treatments of a logistic mixed model, the marginal likelihood, \(L(\boldsymbol{\beta}, \boldsymbol{\Sigma})\) is usually needed. This likelihood, where \(\boldsymbol{\beta}\) is a length \(p\) vector of fixed effects, and \(\boldsymbol{\Sigma}\) is a \(q \times q\) positive definite matrix of (co)-variance components is, \[L(\boldsymbol{\beta}, \boldsymbol{\Sigma}) \propto \int \prod_{i=1}^n {\boldsymbol{p}}_i^{{\boldsymbol{y}}_i}(1-{\boldsymbol{p}}_i)^{1-{\boldsymbol{y}}_i} |{\boldsymbol{\Sigma}}|^{-1/2}e^{-{\boldsymbol{u}}'{\boldsymbol{\Sigma}}^{-1}{\boldsymbol{u}}/2}{\rm d}{\boldsymbol{u}}\label{eqGLMM}\tag{6}\] In (6 ), \({\boldsymbol{p}}_i = (1+e^{-{\boldsymbol{x}}_i'\boldsymbol{\beta} - {\boldsymbol{z}}_i'{\boldsymbol{u}}})^{-1}\) [1], [16] are the elements of the \(n\) length vector of success probabilities, where \({\boldsymbol{x}}_i, {\boldsymbol{z}}_i\) correspond to the covariate values for observation \(i\) while \({\boldsymbol{u}}\) is a \(q\) length vector of random effects. In the special case where \({\boldsymbol{u}}\) corresponds to a random intercept, and all observations; \(j =1, \ldots {\boldsymbol{n}}_i\), within each random effect level share the same fixed effect covariate vector, [17], [18] then the marginal likelihood simplifies to, \[\begin{align} L(\boldsymbol{\beta}, \boldsymbol{\Sigma}) &=& \prod_{i=1}^q L_i(\boldsymbol{\beta}, \boldsymbol{\Sigma}), \quad \text{where } L_i(\boldsymbol{\beta}, \boldsymbol{\Sigma}) = \int {{{\boldsymbol{n}}_i}\choose{ {\boldsymbol{y}}_i}} {\boldsymbol{p}}_i^{{\boldsymbol{y}}_i}(1-{\boldsymbol{p}}_i)^{n_i-{\boldsymbol{y}}_i} \frac{1}{\sqrt{2\pi\sigma^2_u}}e^{-\frac{ {\boldsymbol{u}}_i^2}{2\sigma^2} }{\rm d}{\boldsymbol{u}}. \label{eqGLMM2} \end{align}\tag{7}\] By the binomial theorem, \({\boldsymbol{p}}_i^{{\boldsymbol{y}}_i}(1-{\boldsymbol{p}}_i)^{n_i-{\boldsymbol{y}}_i}\) can be rewritten as \[{\boldsymbol{p}}_i^{{\boldsymbol{y}}_i} \sum_{k=0}^{ {\boldsymbol{n}}_i-{\boldsymbol{y}}_i} {{{\boldsymbol{n}}_i-{\boldsymbol{y}}_i}\choose{k}} (-{\boldsymbol{p}}_i)^k = \sum_{k={\boldsymbol{y}}_i}^{{\boldsymbol{n}}_i} {{{\boldsymbol{n}}_i-{\boldsymbol{y}}_i}\choose{k -{\boldsymbol{y}}_i}} (-1)^{k-{\boldsymbol{y}}_i}{\boldsymbol{p}}_i^k. \label{pbinom}\tag{8}\] Substituting (8 ) into (7 ) implies that \(L_i(\boldsymbol{\beta}, \boldsymbol{\Sigma})\) as defined in (7 ) corresponds to a linear combination of logit-normal moments from order \({\boldsymbol{y}}_i\) to \({\boldsymbol{n}}_i\) with location parameter \({\boldsymbol{x}}_i\boldsymbol{\beta}\) and scale parameter \(\sigma_u\). As the value of \({\boldsymbol{n}}_i\) is arbitrary, this example provides our motivation to find a formula for a logit-normal moment of any positive integer order.

Note the model we are describing is also known as an observation level random effect model, often seen in ecology [19].

2.3 An approximate Bayesian problem: Expectation Propagation for logistic regression↩︎

Expectation Propagation (EP) is a method of approximate Bayesian inference introduced by Minka [6], [7]. In EP, it is assumed the posterior distribution \(p(\boldsymbol{\theta}|{\boldsymbol{y}})\) can be factorised as, \[p(\boldsymbol{\theta}|y_1, \ldots y_n) = \prod_{i=0}^n f_i(\boldsymbol{\theta}),\]

where \(i\) represents a data-point. As the posterior is proportional, up to a constant, to likelihood \(\times\) prior, the natural choice for \(f_i(\boldsymbol{\theta})\) is that when \(i=0\), \(f_i(\boldsymbol{\theta})\) will be proportional to the prior distribution \(p(\boldsymbol{\theta})\), while for \(i = 1, 2, \ldots n\), \(f_i(\boldsymbol{\theta})\) will be proportional to the likelihood of data-point \(i\), \(p(y_i|\boldsymbol{\theta})\). The true posterior \(p(\boldsymbol{\theta}|y_1, \ldots y_n)\) is then approximated with the density \(g(\boldsymbol{\theta})\) such that \(g(\boldsymbol{\theta}) = \prod_{i=0}^ng_i(\boldsymbol{\theta})\). The parameters of the approximate posterior \(g(\boldsymbol{\theta})\) are found by minimising the following Kullback-Leibler divergence \(D_{KL}\{cg_{-i}(\boldsymbol{\theta})f_i(\boldsymbol{\theta})||g(\boldsymbol{\theta})\} =\) \[\int cg_{-i}(\boldsymbol{\theta})f_i(\boldsymbol{\theta}) \{\log(cg_{-i}(\boldsymbol{\theta})f_i(\boldsymbol{\theta}))-\log(g(\boldsymbol{\theta}))\} {\rm d}\boldsymbol{\theta},\] where \(c\) is a normalising constant.

In a generalised linear model (GLM) with fixed dispersion parameter, such as logistic regression, the parameters \(\boldsymbol{\theta}\) are simply the vector of regression coefficients \(\boldsymbol{\beta}\). In this case, the obvious choice is to assume the approximate posterior is normal \(g(\boldsymbol{\beta}) := \mathcal{N}(\boldsymbol{\beta}|\hat{\boldsymbol{\beta}}, \hat{\boldsymbol{\Sigma}})\). In the resulting moment matching algorithm used to determine \(\hat{\boldsymbol{\beta}}, \hat{\boldsymbol{\Sigma}}\) (see [20] for details), the following univariate integrals must be evaluated, \[E_k = \int \boldsymbol{\eta}_i^kg_{-i}(\boldsymbol{\eta}_i)f_i(\boldsymbol{\eta}_i) {\rm d}\boldsymbol{\eta}_i, \quad k = 0, 1, 2, \label{eq46EP}\tag{9}\] where if \(i = 1, \ldots, n\), \(\boldsymbol{\eta}_i = {\boldsymbol{x}}_i'\boldsymbol{\beta}\), \(g_{-i}(\boldsymbol{\eta}_i)\) is a normal density with mean \(\mu_i={\boldsymbol{x}}_i'\hat{\boldsymbol{\beta}}\) and variance \(\sigma^2_i={\boldsymbol{x}}_i'\hat{\boldsymbol{\Sigma}}{\boldsymbol{x}}_i\) and \(f_i(\boldsymbol{\eta}_i)\) is the likelihood term associated with observation \(i\). In most cases, including logistic regression, \(E_k\) is evaluated using numerical integration methods. This is often emphasised for logistic regression, for example in [21], [22] and [23], because the alternative of probit regression is one case where \(E_k\) is analytically known.

To demonstrate that in binary logistic regression, (9 ) is equivalent to a collection of logit-normal moments, we need to re-write (9 ) using repeated application of Stein’s Lemma [24], which states:

This leads to an alternative set of univariate integrals for solve in the EP moment matching, as stated below.

For a binary logistic regression, \(f_i(\boldsymbol{\eta}_i)\) corresponds to \({\boldsymbol{p}}_i = (1+ e^{-{\boldsymbol{x}}_i'\boldsymbol{\beta}})^{-1}\), a logit-normal random variable with location parameter \({\boldsymbol{x}}_i'\hat{\boldsymbol{\beta}}\) and scale parameter \(\sqrt{{\boldsymbol{x}}_i'\hat{\boldsymbol{\Sigma}}{\boldsymbol{x}}_i}\) when \(y_i = 1\). When \(y_i = 0\), \(f_i(\boldsymbol{\eta}_i)= 1 - {\boldsymbol{p}}_i\), i.e. 1 minus a logit-normal variable. Lastly, we need an additional result, used to justify the reparametrisation trick used in machine learning [25],

Using Result [hxExpect], we can re-write \(I_k\) as \(\frac{\rm d}{{\rm d} \mu^k}\big\{ \int g_{-i}(\boldsymbol{\eta}_i)f_i(\boldsymbol{\eta}_i) {\rm d}\boldsymbol{\eta}_i \big\}; k = 1, 2\). For a binary logistic regression, this means determining the approximate EP posterior requires finding a set of \(n\) logit-normal first moments, \(E({\boldsymbol{p}}_i)\) along with its two derivatives w.r.t \(\mu _i = {\boldsymbol{x}}_i'\hat{\boldsymbol{\beta}}\). By the recurrence relation given in (4 ), this is equivalent to finding the first three moments of \(n\) logit normal random variables. Hence, even if analytic evaluation of logit-normal moments of arbitrary order proves impossible, the case of expectation propagation demonstrates determining formulae for a collection of low order logit-normal moments would remain valuable for inferential problems.

2.4 Numerical implementation issues with previous analytic attempts at finding \(E(P)\)↩︎

In theory, [14] provided a somewhat user friendly analytic expression for \(E(P)\). This states the first moment, or expected value is: \[E(P) = \frac{1}{2} + \frac{\sum_{n=1}^\infty e^{-{\sigma^2n^2}/{2}}\sinh(n\mu)\tanh(\frac{n\sigma^2}{2}) + {2\pi}/{\sigma^2}\sum_{n=1}^{\infty}\frac{e^{-{(2n-1)^2\pi^2}/{2\sigma^2}}\sin\{(2n-1)\pi\mu/{\sigma^2}\}}{\sinh\{{(2n-1)\pi^2}/{\sigma^2}\}} }{1+2\sum_{n=1}^\infty e^{-{\sigma^2n^2}/{2}}\cosh(n\mu)}. \label{EP46new}\tag{11}\]

In practice though, the inclusion of \(\cosh(n\mu)\) and \(\sinh(n\mu)\) terms in (11 ) makes this expression for \(E(P)\), and implicitly for all higher order moments unusable in some cases. As \(n \rightarrow \infty\), these terms, unless \(\mu =0\), will blow out to \(\pm \infty\) as \(\cosh(n\mu) \rightarrow \infty \forall \mu \in \mathbb{R}\backslash \{0\}\) and \(1\) if \(\mu = 0\) and \(\sinh(n\mu) \rightarrow \infty \text{ if } \mu > 0, -\infty \text{ if } \mu < 0\) and \(0\) if \(\mu = 0\). While the \(\cosh(n\mu), \sinh(n\mu)\) terms in (11 ) are multiplied by \(e^{-{\sigma^2n^2}/{2}}\), which goes to zero as \(n \rightarrow \infty\) as \(\sigma^2 > 0\), if \(|\mu/\sigma^2|\) is large \(|\cosh(n\mu)|,|\sinh(n\mu)| \rightarrow \infty\) at a lower value of \(n\) than \(e^{-{\sigma^2n^2}/{2}} \rightarrow 0\), causing numerical instability.

For the calculation of moments, [14] developed a truncated version of (11 ) based on the decay of the component series which defined the maximum index \(N\) in the truncation as \(N = \max(\text{ceiling}(|\mu/\sigma^2 -\sqrt{2\log(1/\epsilon)}/\sigma|),\text{ceiling}(|\mu/\sigma^2 +\sqrt{2\log(1/\epsilon)}/\sigma|), \text{ceiling}(1/2 + \sigma\sqrt{\log(2/\epsilon)/2}/\pi))\), with \(\epsilon = 1 \times 10^{-4}\). To illustrate why this truncation method will break down because of the inclusion of \(\cosh(n\mu), \sinh(n\mu)\) terms, consider the case where \(\mu =1 \Leftrightarrow \text{Median}(P) = 0.7311\) to 4 d.p., \(\forall \sigma^2 \in \mathbb{R}^+\). We wish to calculate \(E(P)\) for \(\sigma^2= (1, 0.1, 0.01, 0.001) \Leftrightarrow \sigma = (1, 0.3162, 0.1, 0.0316)\) and will use both numerical integration and the truncated sum version of (11 ). Results of the calculation in R , including the largest index used in the truncated sum are given in Table 1.

Importantly when \(\mu=1, \sigma = 0.0316\), where the highest index in the truncated sum is \(N =1136\), the representation of \(E(P)\) given in (11 ) could not give a meaningful evaluation. Changing the truncation value \(N\) used in the calculation of \(E(P)\), as shown in the later rows of Table 1, we found the breakdown value of \(N\) occurred before an accurate estimate of \(E(P)\) was found.

We could try to get around this by working, when \(\mu > 0\), with \[\begin{align} e^{-\frac{\sigma^2n^2}{2}+\log(\sinh(n\mu))+\log(\tanh(\frac{n\sigma^2}{2}))} &=& e^{-\frac{\sigma^2n^2}{2}-\log(2) +n\mu+\log(1-e^{-2n\mu})+\log(\tanh(\frac{n\sigma^2}{2}))} \nonumber \\ e^{-\frac{\sigma^2n^2}{2}+\log(\cosh(n\mu))} &=& e^{-\frac{\sigma^2n^2}{2}-\log(2)+n\mu+\log(1+e^{-2n\mu})}\nonumber \end{align}\] and \(e^{-\frac{\sigma^2n^2}{2}-\log(2) -n\mu+\log(1-e^{2n\mu})+\log(\tanh(\frac{n\sigma^2}{2}))}, e^{-\frac{\sigma^2n^2}{2}-\log(2)-n\mu+\log(1+e^{2n\mu})}\) when \(\mu < 0\). However, the high value of \(N\) needed when \(\sigma \rightarrow 0\) indicates (11 ) would no longer always out-perform numerical integration on speed. Hence we must conclude the Mordell integral based representation of \(E(P)\) in (11 ) is unsuitable for use in inferential problems, such as those discussed in sections 2.2 and especially 2.3. Hence we need to develop an alternative method for evaluating logit-normal moments.

3 A new approach to finding moments of a logit-normal random variable↩︎

Our new approach for determining the moments of the logit-normal distribution is motivated by [26]. This work included an approximation the evidence lower bound (ELBO), the expected value of the joint distribution with respect to the approximate posterior for a variational Bayes implementation of logistic regression. For Bernoulli/binomial likelihoods with a normal distribution for the approximate posterior this bound contains \[\begin{align} E(S) = E(\log(1+e^X)),\label{softplus46new} \end{align}\tag{12}\] where \(X\sim \mathcal{N}(\mu,\sigma^2)\), and \(S = \log(1+e^X)\) is also known as the softplus function [27].

Our interest in the softplus function, \(s(x)= \log(1+e^x)\), is that its first derivative is the logistic function, \[s'(x) = {e^x}/{(1+e^x)} = (1+e^{-x})^{-1} = p(x).\] If \(x\) is a realisation of a normal random variable with mean \(\mu\) and variance \(\sigma^2\), then \(p(x)\) is a realisation of a logit-normal random variable. This means that \(E(p(X)) = E(P)\), the expected value of a logit normal random variable with location and scale parameters \(\mu\) and \(\sigma\) respectively.

3.1 Re-writing \(s(x), p(x)\) into an expectation taking friendly form↩︎

The properties of the softplus function we need are those used to approximate \(E[s(X)]\) in [26]. First, we split \(s(x) = \log(1 +e^x)\) into even and odd functions. For this, start by rewriting \(s(x)\) as \[\begin{align} s(x) = \log(1+e^x) = \log[(e^{-x} + 1)e^x] = \log(1+e^{-x}) + \log(e^x) = x + \log(1 + e^{-x}), \label{otherSPdef} \end{align}\tag{13}\] then express \(s(x)\) as a piece-wise function using \(s(x) = \log(1+e^x)\) when \(x < 0\) and (13 ) when \(x \geq 0\), i.e. \[\begin{align} s(x) = \log(1+e^x) &=& \begin{cases} \log(1+e^x) \quad \text{If x < 0} \\ x+ \log(1+e^{-x}) \quad \text{If x \geq 0} \\ \end{cases} \end{align}\] or more compactly \[\begin{align} s(x) = x\mathbb{1}_{x > 0} + \log(1+e^{-|x|}), \label{rampSP} \end{align}\tag{14}\] where \(x\) is an odd, and \(\log(1+e^{-|x|})\) is an even function. Taking the derivative of \(s(x)\) as written in (14 ), we get an alternative representation of the logistic function, \[\begin{align} p(x) = (1+e^{-x})^{-1} = \mathbb{1}_{x>0} + (-1)^{\mathbb{1}_{x>0}}(1+e^{-|x|})^{-1}. \label{rampimpliedL} \end{align}\tag{15}\]

To complete the transformation of \(s(x), p(x)\) into a form where we can take approximate expectations analytically, we used the Maclaurin series of \(\log(1+y)\) [26], \[\begin{align} \log(1+y) = \sum_{i=1}^\infty \frac{(-1)^{i-1}y^i}{i}. \label{MSlog1y} \end{align}\tag{16}\]

The Maclaurin series for \(\log(1+y)\) is convergent when \(|y| \leq 1\), and as we assume \(y = e^{-|x|}, x \in \mathbb{R}\), (16 ) is a convergent series \(\forall x\). Hence we will work with the softplus and logistic functions in the following form, \[\begin{align} s(x) = \log(1+e^x) &\approx& x\mathbb{1}_{x > 0} + \sum_{i=1}^N \frac{(-1)^{i-1}e^{-i|x|}}{i}, \tag{17} \\ p(x) = (1+e^{-x})^{-1} &\approx& \mathbb{1}_{x > 0} + (-1)^{\mathbb{1}_{x>0}}\sum_{i=1}^N (-1)^{i-1}e^{-i|x|}.\tag{18} \end{align}\]

This will be our starting point for determining the moments of the logit-normal distribution.

3.1.1 Improving the approximation of \(p(x)\)↩︎

We now will establish a naive copying of [26], i.e. using (18 ), will not lead to a good approximation of \(p(x)\) for all \(x \in \mathbb{R}\). Consider \(p(x)\) when \(x = 0\), which we know is equal to 0.5. In this case the approximation of \(p(x)\) given in (18 ) can be simplified as \[\begin{align} p(x) \approx (-1)^{\mathbb{1}_{x>0}}\sum_{i=1}^N (-1)^{i-1},\label{ramplogisticKomo2} \end{align}\tag{19}\] where \(\sum_{i=1}^\infty (-1)^{i-1}\) is a divergent series whose partial sums oscillate between \(1\) if \(N\) is odd, and 0 if \(N\) is even. More generally, (18 ) contains a modified finite geometric series truncated at \(N-1\) with \(r = -e^{-|x|}\), \[\begin{align} \sum_{i=1}^N (-1)^{i-1}e^{-i|x|} = e^{-|x|}\sum_{i=0}^{N-1} (-e^{-|x|})^i = -r \sum_{i=0}^{N-1}r^i.\label{ramplogisticKomo3} \end{align}\tag{20}\] To determine \(E(P)\), we want to estimate \({-r}/{(1-r)}\), but \(-r \sum_{i=0}^{N-1}r^i\) estimates \({-r(1-r^N)}/{(1-r)}\). The resulting error in the approximation of \(p(x)\) in (18 ) is \({(-1)^{N}e^{-(N+1)|x|}}/{(1+e^{-|x|})}\). The absolute error function of \(p(x)\), \({e^{-(N+1)|x|}}/{(1+e^{-|x|})}\), when approximated by (18 ), is a product of two functions, \(e^{-N|x|}\), \({e^{-|x|}}/{(1+e^{-|x|})}=(1+e^{|x|})^{-1}\), that are monotonely decreasing with respect to \(|x|\). This means the absolute error in the approximation of \(p(x)\) accumulates in the region near \(x \rightarrow 0\) or equivalently \(p(x) \rightarrow 0.5\). In turn, estimation of \(E(P)\) with (18 ) will be particularly hard when \(\mu \approx 0\) and \(\sigma^2 \rightarrow 0\). This makes the implied estimate of \(E(P)\) we could construct by naive following of [26], like the Mordell integral representation of \(E(P)\) in [14], a particular risk if used in parameter estimation algorithms for large datasets.

To reduce the error in the approximation of \(p(x)\) in the region \((-L,L)\), we propose using the piece-wise function, \[\begin{align} p(x) \approx \mathbb{1}_{x > 0} + (-1)^{\mathbb{1}_{x>0}}\sum_{i=1}^{N} ((-1)^{i-1}\mathbb{1}_{|x| \geq L}+ic_i\mathbb{1}_{|x| < L})e^{-i|x|}, \label{ramplogistic46H} \end{align}\tag{21}\] where if \(|x| \geq L\), we use a Maclaurin series and if \(|x| < L\), we use a Chebyshev polynomial interpolation [28]. To calculate the coefficients of the Chebyshev interpolating polynomial, \(c_i\), we will use the R package pracma [29], but applied to the function \(\log(1+y)\) on the interval \((e^{-L},1)\) where \(y =e^x\), then take the derivative. By determining the interpolating polynomial for \(\log(1+y)\), rather than \(\frac{y}{1+y}\), we ensure the constant term \(c_0\) disappears in the approximation of \(p(x)\), leaving only \(c_i\) associated with \(e^{-i|x|}\).

Details on how to choose optimal \(L\) conditional on \(N\) will be presented in section 5. First though, we present formulae for approximating \(E(P)\) as well as higher order moments.

4 (Approximate) moments of a logit-normal random variable↩︎

4.1 The expected value of a logit-normal random variable↩︎

Unlike the Mordell integral based representation of \(E(P)\) in (11 ), the expression for \(E(P)\) given in Proposition [EP46LN2new] will not run into numerical calculation issues. While \(E(P)\) in (22 ) contains terms, \(e^{\pm \mu i+\sigma^2i^2/2}\), which \(\rightarrow \infty\) as \(i \uparrow\), these are multiplied by terms of the form \(\Phi(a-i\sigma)\), with \(a\) constant, which \(\rightarrow 0\) as \(i \uparrow\) at a faster rate than \(e^{\pm \mu i+\sigma^2i^2/2}\) \(\rightarrow \infty\). This means any numerical \(\infty\) that appears in calculation is multiplied by an effective zero.

4.2 Higher order moments↩︎

Higher order moments of \(P\) can be found as a function of \(E(P)\) and its derivatives with respect to \(\mu\), by applying the recursion formula given in (5 ) repeatedly.

By Result [hxExpect], we just need the derivatives of \(p(x)\) from (21 ) and take expectations to approximate higher order moments. As \(p(x)\) is a linear combination of exponentials, the derivatives of \(p(x)\) are highly repetitive, with \[\begin{align} \frac{{\rm d}^k p(x)}{{\rm d}x^k} \approx \begin{cases} \sum_{i=1}^{N} i^k((-1)^{i-1}\mathbb{1}_{|x| \geq L}+ic_i\mathbb{1}_{|x| < L})e^{-i|x|} \quad \text{ If k is odd,} \\ (-1)^{\mathbb{1}_{x>0}}\sum_{i=1}^{N} i^k((-1)^{i-1}\mathbb{1}_{|x| \geq L}+ic_i\mathbb{1}_{|x| < L})e^{-i|x|} \quad \text{If k is even.} \\ \end{cases} \label{ramplogistic46diff} \end{align}\tag{23}\]

To find \(\frac{{\rm d}^kE(P)}{{\rm d}\mu^k}\), we just change the leading coefficients of all components except \(\Phi(\mu/\sigma)\) in the approximation of \(E(P)\) from (22 ), \[\begin{align} \frac{{\rm d}^kE(P)}{{\rm d}\mu^k} &=& \sum_{i=1}^{N} i^k(-1)^{i-1}\left\{e^{\mu i+\sigma^2i^2/2}\Phi\left(\frac{-L-\mu-\sigma^2 i}{\sigma}\right) +(-1)^{k-1}e^{-\mu i+\sigma^2i^2/2}\Phi\left(\frac{-L+\mu-\sigma^2 i}{\sigma}\right) \right\} \nonumber \\ &&+\sum_{i=1}^{N} i^{k+1}c_ie^{\mu i+\sigma^2i^2/2}\left\{\Phi\left(\frac{-\mu-\sigma^2 i}{\sigma}\right) - \Phi\left(\frac{-L-\mu-\sigma^2 i}{\sigma}\right)\right\}\nonumber \\ &&+(-1)^{k-1}\sum_{i=1}^{N} i^{k+1}c_ie^{-\mu i+\sigma^2i^2/2}\left\{\Phi\left(\frac{\mu-\sigma^2 i}{\sigma}\right)-\Phi\left(\frac{-L+\mu-\sigma^2 i}{\sigma}\right) \right\}. \label{eq46EOdiff} \end{align}\tag{24}\]

Substituting (24 ) sequentially into (5 ), starting with \(k=1\) to find the second moment, leads to the following proposition for an approximation of positive integer moments,

A proof of Proposition [EP46LN2high] can be found in the appendix.

5 Determining optimal parameters for \(N\) and \(L\)↩︎

We have not yet established how accurate our approximations for \(E(P), E(P^2), \ldots\), will be, or what optimal choices for \(N, L\) are. While increasing \(N\) should result in improved accuracy of approximation, we still need to find the optimal \(L\) given \(N\). To determine optimal \(L\), we will work with the underlying functions \(p(x), p'(x), \ldots\), as if we can guarantee accurate approximation of the underlying functions for all \(x\), we can guarantee accuracy of moment estimates.

To choose \(N, L\), we consider accuracy of approximation up to the third derivative of \(p(x)\). This implies we are aiming to achieve accurate approximation up to at least the fourth moment. As the approximate component of \(p(x)\) in (21 ) is symmetric around zero, to determine optimal \(L\), we only need to look at \(x \in \mathbb{R}^+\) or \(x \in \mathbb{R}^-\). As at \(|x| = 5, N=5\), the absolute error of the Maclaurin series component of \(p(x)\) is \(e^{-30}/(1+e^{-5}) = 9.3 \times 10^{-14}\), we will focus approximation on \(x \in (-5,0)\) and \(N \geq 5\). We then carried optimisation of \(L\), using Brent’s method [30] as implemented using the optimize function in R. The cost function we choose to minimise is \(\sum_j |p^{(k)}(x_j)-\tilde{p}^{(k)}(x_j)|\), where \(p^{(k)}(x_j),\tilde{p}^{(k)}(x_j)\) is the \(k^\text{th}; k =0, 1, 2, 3\) derivative of the logistic function \(p(x)\), and the approximation of the logistic function given in (21 ) respectively, with \(x_{j+1} - x_{j} = 1 \times 10^{-4}\).

In Figure 1, we present approximation error for \(p(x)\) and its first three derivatives, while in Table 2, we present the optimal choice of \(L\) for \(N=5, \ldots, 11\) along minimised sum of absolute errors and maximum error for the same functions. Note allowing \(L\) to vary for each derivative of \(p(x)\) reduces the efficiency of any algorithm to calculate moments. This is because, 1. we would need to calculate a new set of evaluations of \(\Phi(\cdot)\) for the expected value of each derivative, and 2. we would not be able to use the result of proposition [EP46LN2high] for higher order moments. Hence, to optimise \(L\), we fix \(N=11\), as we can guarantee the maximum error in the approximation of the four underlying functions needed to evaluate the first four moments of the logit-normal distribution over all \(x\) is \(<2.7 \times 10^{-4}\), and find \(L\) working with the sum of \(p(x)\) and its first three derivatives. This resulted in our final choice of \(L\) being \(1.507306\), with the values of the Chebyshev polynomial coefficients \(c_i\) given in Table 3.

6 Testing accuracy and speed of moment calculation↩︎

6.1 Accuracy↩︎

Figure 1, by showing we have accurately approximated \(p(x)\) and its first three derivatives indicates our approach should be accurate for estimating \(E(P^k); k = 1, 2, 3, 4\) for any value of \(\mu, \sigma\). Nevertheless, to confirm the accuracy of our method, we will estimate \(E(P), E(P^2), E(P^3), E(P^4)\) using equations (22 ), (25 ) from section 4.1 and 4.2, with values of \(L, c_i\) fixed at those given in section 5 and compare to numerical integration estimates. As the logit-normal distribution has the symmetry \(P \sim \text{logit}N(\mu,\sigma^2)\) then \(1-P \sim \text{logit}N(-\mu,\sigma^2)\), we will only consider values \(\mu \in [0,6]\) in increments of 0.25, which corresponds to \(\text{Median}(P) \in [0.5,0.9975]\). For \(\sigma\) we will consider the values \((0.001, 0.01, 0.1, 0.25, 0.5, 1, 1.5, 2.5)\). This means we cover, when \(\mu\) is near zero, scenarios from near normality for \(P\) (low \(\sigma\)), to near uniform all the way to bi-modality (high \(\sigma\)) with modes at 0 and 1, as shown in Figure 2 (a).

In Tables 4-5, comparisons between estimates from our piece-wise approximations of order \(N=11\) and numerical integration are provided for the expected value and variance of a logit-normal random variable. Tables showing accuracy for the second to fourth moment, skewness and kurtosis are included in the appendix. For the first four moments and variance, we did obtain accurate estimates for all combinations of \(\mu, \sigma\) considered, with the average and worst case error being \(6.1 \times 10^{-9}\) and \(1.2 \times 10^{-6}\) for \(E(P)\), \(1.4 \times 10^{-8}\) and \(2.2 \times 10^{-6}\) for \(\text{Var}(P)\), \(2.6 \times 10^{-9}\) and \(7.1 \times 10^{-8}\) for \(E(P^2)\), \(1.2 \times 10^{-8}\) and \(2.2 \times 10^{-7}\) for \(E(P^3)\) and \(5.8 \times 10^{-7}\) and \(4.4 \times 10^{-5}\) for \(E(P^4)\). For skewness and kurtosis, the division by a power of \(\text{Var}(P)\) means when \(\sigma \leq 0.25\), estimates obtained are unreliable.

Our second test looked at the applicability of Proposition [EP46LN2high]. In practice, this is a test of how good the underlying approximating functions implied in Proposition [EP46LN2high] are. In Figure 2, we provide evidence of our ability to approximate \(p(x)^k\) for \(k = 5, 6, 7, 8\). For \(k=5, 6\), our choice of \(c_i, L\) remains good enough to enable accurate estimation of the corresponding logit-normal moment. For \(k = 7\) and more so for \(k=8\), we would start seeing systematic error appear when estimating the \(k^\text{th}\) moment for a logit-normal random variable with non-trivial probability in the region \(0.5 \pm \epsilon\) for some \(\epsilon > 0\). If we push \(k\) higher, approximations will get worse. To show this, we calculated \(E(P^{k}); k = 1, 2, \ldots, 10\) using both the proposition and numerical integration when \(\mu =1, \sigma = 0.8\). As shown in Table 6, the Proposition result does hold for \(k > 4\). However, the increasing error in approximation for higher order derivatives of \(p(x)\), shown in Figure 1, means approximation accuracy of \(E(P^k)\) declines noticeably when \(\mu =1, \sigma =0.8\) for \(k \geq 8\), if \(N =11, L = 1.507306\). This implies our method for finding logit-normal moments is not generally applicable to evaluating the integral needed to find the marginal likelihood of \(\boldsymbol{\beta}\), \(L(\boldsymbol{\beta})\) that appears in binomial logistic mixed models with only an observation level random effect.

6.2 Speed↩︎

To evaluate speed of calculation, we have chosen to focus on \(E(P), E(P^2), E(P^3), E(P^4), \text{Var}(P)\). Because of the highly sequential nature of the component functions needed to calculate our estimates of logit-normal moments, the marginal cost of calculating multiple moments simultaneously drops as the number of moments needed increases. This means we will provide two speed comparisons. The first comparison consists of calculating a single moment, \(E(P)\), using our method and numerical integration. The second will be a comparison of calculating \(E(P), E(P^2), E(P^3), E(P^4), \text{Var}(P)\) using our method and numerical integration, with results presented in Table 7. For speed testing, we implemented our approach in R but written in C++ code using the package Rcpp [31], while numerical integration was done using the integrate function in R. To determine the time distribution, the functions were evaluated 10,000 times and timed using the R package microbenchmark [32]. The values for \(\mu, \sigma\) were \(\mu = 0.2, 1.2, \sigma = 0.1, 1, 2.5\). This was done to see if speed of evaluation was affected by the choice of parameter values.

In Table 7, we see evaluation time for \(E(P)\) using numerical integration is somewhat sensitive to choice of \(\sigma\), with evaluation being faster for low \(\sigma\). This was less so for our approach based on piecewise approximation of \(p(x)\). Our approach for finding \(E(P)\) was roughly 5.5-6 times faster when \(\sigma = 0.1\) and was roughly 7-9 times faster when \(\sigma =1, 2.5\). When calculating the first four moments and variance simultaneously, the low marginal cost of obtaining higher order moments and the variance using our approach becomes apparent. On average, obtaining \(E(P^2), E(P^3), E(P^4), \text{Var}(P)\) in addition to \(E(P)\) increased run time by 2 to 3 %, or \(0.1-0.2 \times 10^{-6}\) s. Using standard numerical integration to obtain \(E(P^2), E(P^3), E(P^4), \text{Var}(P)\) in addition to \(E(P)\) caused a linear increase in coast, as numerical solving of four rather than one integral was needed. Our approach to finding \(E(P), E(P^2), E(P^3), E(P^4), \text{Var}(P)\) was on average between 23-25 times faster when \(\sigma =0.1\), 31-35 times faster when \(\sigma =1\) and 39-40 times faster when \(\sigma=2.5\). The low cost of finding higher order moments/expectation with respect to \(\mu\) of \(E(P)\) using our piecewise approximation to \(p(x)\) indicates our approach would be particularly useful in algorithms, such as logistic regression models where parameters are estimated using expectation propagation, where we need moments of multiple order simultaneously.

7 An inferential application that will work: Expectation propagation for logistic regression↩︎

Our method for approximating logit-normal moments has been shown to have a practical upper limit of 7/8 before losing accuracy. However there are many inferential problems where logit-normal moments under order 7/8 are needed. As an example, we will now demonstrate our results are sufficient to remove the need to use numerical integration when fitting logistic regression in using expectation propagation for many choices of prior. As the illustrative dataset, we chose the lbw dataset from the COUNT package in R, originally given in [33].

7.1 Dataset and variables↩︎

The lbw dataset consists of 189 observations with 10 variables. For the response variable, we have chosen to use low, an indicator (1 = yes, 0 = no) whether a birth was classed as low-weight or not . For our illustrative example, we have chosen to consider only 3 of the possible predictor variables. These are Smoke, a binary variable indicating if the mother had a history of smoking (1 = yes, 0 = no), lwt, a continuous variable giving maternal weight at last menstrual period in pounds, and Hypertension, a binary variable indicating if the mother had a history of hypertension (1 = yes, 0 = no). For the continuous variable lwt, we will standardise the covariate values to create Standardised Maternal Weight before running the algorithm,.

7.2 Estimation methods↩︎

Our primary interest in implementing expectation propagation (EP) is to demonstrate that for logistic regression, rewriting the moment matching step into the form of Proposition [prop1] leads to faster implementation, as we know the required logit-normal moments. This means we will implement two versions of EP, written in R code. The first, EP-NI, will implement sequential-EP, as described in [20] using numerical integration to solve the moment matching step. The second, EP-LN, will implement sequential-EP using our formula for \(E(P), \frac{{\rm d}E(P)}{{\rm d} \mu}, \frac{{\rm d}^2E(P)}{{\rm d} \mu^2}\) given in (22 ), (24 ) to avoid numerical integration. For more details on the particulars of the algorithm, particularly for accounting that \(f_i(\boldsymbol{\eta}_i)\) switches between \({\boldsymbol{p}}_i\) and \(1-{\boldsymbol{p}}_i\) when \({\boldsymbol{y}}_i\) changes from 1 to 0, see the appendix. While not of main interest, we will include a MCMC comparison to demonstrate the accuracy of expectation propagation at posterior approximation. This comparison will be with a Gibbs sampler using Pólya-Gamma data augmentation [34]. For all MCMC analyses, 3 chains will be run for 30,000 iterations with the first 15,000 iterations removed as burn-in. Convergence will be checked using Gelman-Rubin statistics [35].

For the prior on the regression coefficients, we choose \(p(\boldsymbol{\beta}) := \mathcal{N}(\boldsymbol{\beta}|\tilde{\boldsymbol{\beta}}, \tilde{\boldsymbol{\Sigma}})\), where \(\tilde{\boldsymbol{\beta}} = (0, -0.1, 0.1, 0.2)\), and \(\tilde{\boldsymbol{\Sigma}}_{ij} = 1.3\times 0.5^{|i-j|}\). There are two motivations for this. The first is with a normal prior, the moment matching step required in the EP algorithm can skip over the prior [21]. Second, in applications of EP to determine the marginal likelihood in frequentist probit mixed models [23], the prior equivalent is also normal. Moreover [23] state the extension of their algorithm to a logistic mixed model requires solving univariate integrals equivalent to (9 ) with \(f_i(\boldsymbol{\eta}_i) = (1-{\boldsymbol{p}}_i)\). Hence showing our approach works for logistic regression with a normal prior implies our work could be extended to the logistic examples mentioned in [23].

7.3 Results↩︎

In Figure 3, we present posterior density plots estimated using MCMC and EP-NI, EP-LN. As expected, very high accuracy of approximation was obtained. This was confirmed by calculating the accuracy diagnostic from [36], \[\text{accuracy}(g(\boldsymbol{\beta}_i)) = 100 \times \bigg( 1 - 0.5 \int_{-\infty}^\infty |g(\boldsymbol{\beta}_i)-p(\boldsymbol{\beta}_i|{\boldsymbol{y}})| d\boldsymbol{\beta}_i \bigg),\] where \(g(\boldsymbol{\beta})\) is the approximate posterior and \(p(\boldsymbol{\beta}|{\boldsymbol{y}})\) is the ‘true’ posterior as estimated with MCMC samples. For the 4 parameters in our test model, accuracy of EP approximation was \(0.9882\) for the Intercept, \(0.9926\) for Smoke, \(0.9788\) for Standardised Maternal Weight and \(0.9934\) for Hypertension.

Focusing on the comparison between EP-NI and EP-LN, we see, as expected, that the two algorithms are effectively equivalent, as shown in Table 8. The relative error in estimates of \(\hat{\boldsymbol{\beta}}, \hat{\boldsymbol{\Sigma}}\) between EP-NI and EP-LN are in the range \(10^{-8} - 10^{-9}\) for elements of \(\hat{\boldsymbol{\beta}}\), and \(10^{-6} - 10^{-8}\) for elements of \(\hat{\boldsymbol{\Sigma}}\).

Lastly, to benchmark the computational speed gain obtained by using EP-LN over EP-NI, we ran both versions of EP 1,000 times, with timings determined using the R package microbenchmark [32]. Results are provided in Table 9. This shows that for a sample of size 189, replacing numerical integration with our approximations of logit-normal moments reduced computation time by about 58 %.

8 Conclusions↩︎

We would argue this paper fulfils the conclusion of [14] that ‘the possibility of either analytic positive integer moments, or at least a highly accurate and simple approximation could be determined, through further study that utilises the similarity between the log-normal and logit-normal distributions’. Through examples, we have shown our estimates of moments up to order 4 do not suffer the numerical issues that the Mordell integral based representation of \(E(P)\) from [14] has when \(\sigma\) is small, and thus maintain a high level of accuracy in all settings considered.

We have shown our method for estimating low order moments, up to the fourth, of logit-normal random variables is faster than standard numerical integration in R, especially if multiple moments are needed simultaneously. As error in the underlying approximating function accumulates with increasing moment order, a finite limit exists in practice on the highest moment we can calculate accurately. Referring to the implied approximating function of the \(k^\text{th}\) power of the logistic function, this limit is around 7 or 8. However, the highest moment accurately approximated will also depend on \(\mu, \sigma\), as error in the underlying function \(p(x)^k\) first appears for \(x \approx 0\).

The inability of our approach to estimate moments of arbitrarily high order means our work has limited application as an alternative to quadrature to finding the marginal likelihood in frequentist treatments of the logistic mixed model. However, while not initially obvious, we have demonstrated our approach to finding logit-normal moments does enable faster, numerical integration free implementation of expectation propagation for logistic regression. This is because only knowledge of the equivalent to logit-normal moments of order \(1, 2, 3\) are required to solve the expectation propagation moment matching steps as described in [20] related to the likelihood in logistic regression. We do note expectation propagation algorithms for regression type problems with normal priors have been considered for finding the marginal likelihood for fairly general frequentist treatments of generalised linear mixed models [23]. This indicates our illustrative application may have indirect frequentist as well as Bayesian algorithm uses. Lastly, we suspect our results, apart from enabling faster implementation of expectation propagation, can be utilised in other versions of approximate Bayesian inference, such as variational Bayes methods for logistic and negative binomial regression. These applications are the focus of future work.

John Holmes proposed the study and potential illustrative applications. Ness Arps was involved with parameter testing. Marco Reale was involved in manuscript refinement.

None reported.

The authors declare no conflicts of interest.

9 Expectation Propagation as described in Bayesian Data Analysis↩︎

One of the inferential applications we considered as motivational for researching logit-normal moments was Expectation Propagation (EP) [6], [7] for logistic regression. Specifically, our motivation comes from the EP algorithm for a regression problem with fixed dispersion parameter described in [20], which we outline below.

9.1 EP for a regression problem with fixed dispersion parameter↩︎

As we are considering a regression type problem, we assume the approximate posterior is normal \(g(\boldsymbol{\beta}) := \mathcal{N}(\boldsymbol{\beta}|\hat{\boldsymbol{\beta}}, \hat{\boldsymbol{\Sigma}})\). This choice of \(g\) does admit the factorisation \(g(\boldsymbol{\beta}) = \prod_{i=0}^ng_i(\boldsymbol{\beta})\), which we need as the Kullback-Leibler divergence to be minimised is \[D_{KL}\{cg_{-i}(\boldsymbol{\beta})f_i(\boldsymbol{\beta})||g(\boldsymbol{\beta})\} = \int cg_{-i}(\boldsymbol{\beta})f_i(\boldsymbol{\beta}) \{\log(cg_{-i}(\boldsymbol{\beta})f_i(\boldsymbol{\beta}))-\log(g(\boldsymbol{\beta}))\} {\rm d}\boldsymbol{\beta},\] where \(c\) is a normalising constant.

The steps needed to perform this minimisation are given below. As we assume a normal prior, we let \(g_0(\boldsymbol{\beta})\) be equal to the prior. Using a normal prior also means \(g_0(\boldsymbol{\beta})\) will never need updating [21], so moment matching only needs to be applied to the \(n\) terms related to the likelihood.

• For \(t = 1, \ldots,\)

  • For \(i = 1, \ldots, n\)

    • Compute the (natural) parameters of the cavity distribution \(g_{-i}(\boldsymbol{\beta})\),

    • \(\hat{\boldsymbol{\Sigma}}_{-i}^{-1} = \hat{\boldsymbol{\Sigma}}^{-1} -\hat{\boldsymbol{\Sigma}}_i^{-1}\), \(\hat{\boldsymbol{\Sigma}}_{-i}^{-1}\hat{\boldsymbol{\beta}}_{-i} =\hat{\boldsymbol{\Sigma}}^{-1}\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\Sigma}}_i^{-1}\hat{\boldsymbol{\beta}}_i \Rightarrow \hat{\boldsymbol{\beta}}_{-i} = \hat{\boldsymbol{\Sigma}}_{-i}(\hat{\boldsymbol{\Sigma}}^{-1}\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\Sigma}}_i^{-1}\hat{\boldsymbol{\beta}}_i)\).

    • Compute the parameters of \(g_{-i}(\boldsymbol{\eta}_i)\), where \(\boldsymbol{\eta}_i = {\boldsymbol{x}}_i'\hat{\boldsymbol{\beta}}\) and \({\boldsymbol{x}}_i\) are the covariate values, including intercept, for the \(i^\text{th}\) observation,

    • \(M_{-i} = {\boldsymbol{x}}_i'\hat{\boldsymbol{\beta}}_{-i}\), \(V_{-i} = {\boldsymbol{x}}_i'\hat{\boldsymbol{\Sigma}}_{-i}{\boldsymbol{x}}_i\)

    • Construct the un-normalised tilted distribution and calculate, \[\begin{align} E_k = \int_{-\infty}^{\infty} \boldsymbol{\eta}_i^k g_{-i}(\boldsymbol{\eta}_i)f_i(\boldsymbol{\eta}_i) d\boldsymbol{\eta}_i \quad\text{for k = 0, 1, 2}. \nonumber \end{align}\] Then moment match and set \(M = \frac{E_1}{E_0}\) and \(V = \frac{E_2}{E_0}-M^2\). Note in practice, we will need to specify finite bounds on the integral. The simplest choice would be \(M_{-i} \pm \delta\sqrt{V_{-i}}\) for suitably large \(\delta\).

    • Determine the natural parameters of \(g_i(\boldsymbol{\eta}_i)\),

    • \(M_i/V_i = M/V - M_{-i}/V_{-i}\), \(1/V_i = 1/V - 1/V_{-i}\).

    • Transform the natural parameters found in step 4 to those of \(g_i(\boldsymbol{\beta})\),

    • \(\hat{\boldsymbol{\Sigma}}_{i}^{-1}\hat{\boldsymbol{\beta}}_{i} = {\boldsymbol{x}}_iM_i/V_i, \hat{\boldsymbol{\Sigma}}_{i}^{-1} = {\boldsymbol{x}}_i(1/V_i){\boldsymbol{x}}_i'\)

    • Update the natural parameters of \(g(\boldsymbol{\beta})\),

    • \(\hat{\boldsymbol{\Sigma}}^{-1} = \hat{\boldsymbol{\Sigma}}_i^{-1}+ \hat{\boldsymbol{\Sigma}}_{-i}^{-1}\) and \(\hat{\boldsymbol{\Sigma}}^{-1}\hat{\boldsymbol{\beta}}= \hat{\boldsymbol{\Sigma}}_i^{-1}\hat{\boldsymbol{\beta}}_i + \hat{\boldsymbol{\Sigma}}_{-i}^{-1}\hat{\boldsymbol{\beta}}_{-i}\).
      

  Stop once estimates have converged.

Some versions, known as parallel EP, skip step 6, updating \(g(\boldsymbol{\beta})\) only after \(g_i(\boldsymbol{\beta})\) has been updated for all \(i\). In our example, we will implement sequential EP, which does not skip step 6.

9.1.1 Proving Proposition 1↩︎

Our interest for speeding up inference is focused on avoiding the numerical integrations in step 3, \[\begin{align} E_k = \int_{-\infty}^{\infty} \boldsymbol{\eta}_i^k g_{-i}(\boldsymbol{\eta}_i)f_i(\boldsymbol{\eta}_i) d\boldsymbol{\eta}_i \quad\text{for k = 0, 1, 2}. \label{EQEK} \end{align}\tag{26}\] To achieve this for logistic regression, this requires re-writing (26 ) as shown in Proposition [prop1].

Lastly, we can use result [hxExpect] (see main paper),

This allows us to re-write \(E_1, E_2\) as \[\begin{align} E_1 &=& M_{-i}E(f_i(\boldsymbol{\eta}_i)) + V_{-i}\frac{ {\rm d}E(f_i(\boldsymbol{\eta}_i))}{{\rm d}M_{-i}} \\ E_2 &=& (M_{-i}^2+ V_{-i})E(f_i(\boldsymbol{\eta}_i)) + 2M_{-i}V_{-i}\frac{ {\rm d}E(f_i(\boldsymbol{\eta}_i))}{{\rm d}M_{-i}} + V_{-i}^2\frac{ {\rm d}^2E(f_i(\boldsymbol{\eta}_i))}{{\rm d}M_{-i}^2}. \label{eqlat} \end{align}\tag{30}\]

9.1.2 Implications in the binary logistic regression case↩︎

For a binary logistic regression, with \(\boldsymbol{\eta}_i\) being normal, \(f_i(\boldsymbol{\eta}_i)\) corresponds to \({\boldsymbol{p}}_i = (1+ e^{-\boldsymbol{\eta}_i})^{-1}\), a logit-normal random variable with location \({\boldsymbol{x}}_i\hat{\boldsymbol{\beta}}\) and scale \(({\boldsymbol{x}}_i'\hat{\boldsymbol{\Sigma}}{\boldsymbol{x}}_i)^{0.5}\) parameters when \(y_i = 1\). If \(y_i = 0\), \(f_i(\boldsymbol{\eta}_i)= 1 - {\boldsymbol{p}}_i\). In addition, the logit-normal distribution has the symmetry, if \(P \sim \text{logit}\mathcal{N}(\mu,\sigma^2)\) then \(1- P \sim \text{logit}\mathcal{N}(-\mu,\sigma^2)\). Exploiting this symmetry, for binary logistic regression the integrals in (26 ) needed for moment matching can be compactly written as \[\begin{align} E_0 = E({\boldsymbol{p}}_i), \quad E_1 = (2{\boldsymbol{y}}_i-1)E(\boldsymbol{\eta}_i{\boldsymbol{p}}_i), \quad E_2 = E(\boldsymbol{\eta}_i^2{\boldsymbol{p}}_i) \label{bineq3} \end{align}\tag{31}\] where \({\boldsymbol{p}}_i \sim \text{logit}\mathcal{N}((2{\boldsymbol{y}}_i-1)M_{-i},V_{-i}), \boldsymbol{\eta}_i \sim \mathcal{N}((2{\boldsymbol{y}}_i-1)M_{-i},V_{-i})\). Applying the alternative form of \(E_0, E_1, E_2\) given in Proposition [prop1], specifically (30 ), we have shown implementing EP moment matching for logistic regression just requires a logit-normal first moment along with its first two derivatives w.r.t. \((2{\boldsymbol{y}}_i-1)M_{-i}\).

10 The expected value (first moment) of a logit-normal random variable↩︎

We propose that the expected value of a logit-normal random variable with location parameter \(\mu\) and scale parameter \(\sigma\) is very accurately approximated by the result of Proposition [EP46LN2new]. This result is obtained by determining the expected value of our proposed approximation of the logistic function given in (21 )

11 The \(k^\text{th}\) moment of a logit-normal random variable↩︎

We propose the \(k^\text{th}\) moment of a logit-normal random variable can be found using the result of Proposition [EP46LN2high] (see main paper). We now will prove this result using induction. As part of the proof, we will use the formula for the expected value of the derivatives of \(E(P)\) with respect to \(\mu\) given in the main text. As a reminder, this showed that derivatives of the logistic function can be approximated by (23 ) and the expected value of \(\frac{{\rm d}^kE(P)}{{\rm d}\mu^k}\), can be approximated using the same results for component integrals from (35 38 ) as shown in (24 ).

12 Testing accuracy of approximation for higher order moments↩︎

In the main text, we showed that our approach provided accurate estimates of the expected value and variance of a logit-normal random variable for a wide range of values of \(\mu, \sigma\). In this section, we will show our approximations maintain accuracy up to at least the fourth moment. This is unsurprising as our approximation of the underlying functions is accurate up to the third derivative of the logistic function. We also include calculations for skewness and kurtosis. This have been put in to highlight when the numerator moments are divided by a variance near zero, numerical instability will result, even if the component moments are well estimated.

12.1 Second, third and fourth moments↩︎

12.1.1 Third moment↩︎

12.1.2 Fourth moment↩︎

12.2 Skewness and Kurtosis↩︎

12.2.1 Skewness↩︎

12.2.2 Kurtosis↩︎

13 Code needed in paper↩︎

13.1 R function for finding first to fourth moments and variance↩︎

In this section, we present functions for finding \(E(P)\) and \(E(P), E(P^2), E(P^3), E(P^4), \text{Var}(P)\). Note these are in C++ and users will need the R package Rcpp[31] for compiling the functions, before these can be used in R. The input parameters are \(\mu, \sigma\), but there is an additional C++ function included for the evaluation of the standard normal cdf \(\Phi(\cdot)\).

        #include <Rcpp.h>
        #include <cmath>
        #include <cstddef>
        #include <limits>
        #include <vector>
        
        // normal_cdf is a C++ function for calculating the standard normal CDF.
        
        // [[Rcpp::export]]
        static inline double normal_cdf(const double x) {
            // Standard normal CDF using erfc for numerical stability
            return 0.5 * std::erfc(-x / std::sqrt(2.0));
        }
        
        // EP_logitnormalkpartcplus(double mu, double sigma) is a C++ function for 
        // calculating the expectation of a logit normal random variable with 
        // location parameter mu and scale parameter sigma.
        
        // [[Rcpp::export]]
        double EP_logitnormalkpartcplus(double mu, double sigma) {
            // In part one we calculate the component functions and fixed inputs, coef_nearzero, L and k.
            const double sigma2 = sigma * sigma;
            
            double sum_EP_terms = 0.0;
            double sum_EdiffP = 0.0;
            double sum_E2diffP_terms = 0.0;
            double sum_E3diffP = 0.0;
            double L = 1.507306;
            std::size_t k = 11;
            std::vector<double> coef_nearzero = {0,0.999991820, -0.499903021,  0.332646462, -0.246745120,  
                0.189036359, -0.139389342,  0.091239940, -0.048514467,  0.018953113, -0.004725329, 0.000556454};
            
            for (std::size_t i = 1; i <= k; ++i) {
                const double N = static_cast<double>(i);
                
                double p0a = std::exp(N * mu + 0.5 * N * N * sigma2);
                double p0b = std::exp(-N * mu + 0.5 * N * N * sigma2);
                
                // if element is infinity, the cdf to be multiplied
                // is zero. Adding this line stops NaN.
                if (!std::isfinite(p0a)) p0a = 0.0;
                if (!std::isfinite(p0b)) p0b = 0.0;
                
                const double p0c = coef_nearzero[i] * N;    
                const double p1 = ((i               const double p2 = p0a * normal_cdf((-L - mu - N * sigma2) / sigma);
                const double p3 = p0b * normal_cdf((mu - L - N * sigma2) / sigma);
                const double p4 = p0c * (p0a * normal_cdf((-mu - N * sigma2) / sigma));
                const double p5 = p0c * (p0b * normal_cdf((mu - N * sigma2) / sigma));
                
                // In part two, we use the results of part one to find E(P) and the first three derivatives.
                sum_EP_terms += p1 * (p2 - p3) + p4 - p5;
                
            }
            
            const double EP = normal_cdf(mu / sigma) + sum_EP_terms;
            
            return EP;
        }
        
        // EP1to4_logitnormalkpartcplus(double mu, double sigma) is a C++ function 
        // for calculating the expectation of a logit normal random variable with 
        // location parameter mu and scale parameter sigma, along with the second, 
        // third, fourth moments and the variance. 
        
        // [[Rcpp::export]]
        std::vector<double> EP1to4_logitnormalkpartcplus(double mu, double sigma) {
            // In part one we calculate the component functions.
            const double sigma2 = sigma * sigma;
            
            double sum_EP_terms = 0.0;
            double sum_EdiffP = 0.0;
            double sum_E2diffP_terms = 0.0;
            double sum_E3diffP = 0.0;
            double L = 1.507306;
            std::size_t k = 11;
            std::vector<double> coef_nearzero = {0,0.999991820, -0.499903021,  0.332646462, -0.246745120, 
                0.189036359, -0.139389342,  0.091239940, -0.048514467,  0.018953113, -0.004725329, 0.000556454};
            
            for (std::size_t i = 1; i <= k; ++i) {
                const double N = static_cast<double>(i);
                
                double p0a = std::exp(N * mu + 0.5 * N * N * sigma2);
                double p0b = std::exp(-N * mu + 0.5 * N * N * sigma2);
                
                // if element is infinity, the cdf to be multiplied
                // is zero. Adding this line stops NaN.
                if (!std::isfinite(p0a)) p0a = 0.0;
                if (!std::isfinite(p0b)) p0b = 0.0;
                
                const double p0c = coef_nearzero[i] * N;    
                const double p1 = ((i               const double p2 = p0a * normal_cdf((-L - mu - N * sigma2) / sigma);
                const double p3 = p0b * normal_cdf((mu - L - N * sigma2) / sigma);
                const double p4 = p0c * (p0a * normal_cdf((-mu - N * sigma2) / sigma));
                const double p5 = p0c * (p0b * normal_cdf((mu - N * sigma2) / sigma));
                
                // In part two, we use the results of part one to find E(P) and the first three derivatives.
                sum_EP_terms += p1 * (p2 - p3) + p4 - p5;
                sum_EdiffP += N * (p1 * (p2 + p3) + p4 + p5);
                sum_E2diffP_terms += (N * N) * (p1 * (p2 - p3) + p4 - p5);
                sum_E3diffP += (N * N * N) * (p1 * (p2 + p3) + p4 + p5);
            }
            
            const double EP = normal_cdf(mu / sigma) + sum_EP_terms;
            const double EdiffP = sum_EdiffP;
            const double E2diffP = sum_E2diffP_terms;
            const double E3diffP = sum_E3diffP;
            
            // In part three, we find higher order moments.
            const double EP2 = EP - EdiffP;
            const double EP3 = EP2 - 0.5 * (EdiffP - E2diffP);
            const double EP4 = EP3 - (1.0 / 3.0) * (EdiffP - E2diffP - 0.5 * (E2diffP - E3diffP));
            
            // In part four, we calculate central moments (Var, skewness and kurtosis)
            const double VarP = EP2 - EP * EP;
            
            std::vector<double> results;
            results.reserve(5);
            results.push_back(EP);
            results.push_back(EP2);
            results.push_back(EP3);
            results.push_back(EP4);
            results.push_back(VarP);
            
            return results;
        }

13.2 Other code used in this paper↩︎

Additional code for implementing the EP examples to illustrate EP logistic regression is (1) requires the first-third moment of logit-normal random variables can be found at https://github.com/johnbholmes/Logit-normal-moments-with-EP-applications. This contains R code for implementing expectation propagation, and a Gibbs sampler for MCMC based comparison.

References↩︎