Gauss-Hermite Integration

One-dimensional Gauss-Hermite integration

Gauss-Hermite integration approximates a specific integral as follows

\[\int\limits_{-\infty}^\infty f(x) \exp(-x^2) \mathrm{d}x \approx \sum_{j=1}^N w_{j} f(x_{j})\]

in which \(w_{j}\) and \(x_{j}\) are weights and abscissas for N-point Gauss-Hermite integration, see Abramowitz and Stegun (1972). N-point integration is exact for all polynomials \(f(x)\) of degree 2N-1, see Dahlquist and Bjorck (1974). This can for instance be used to approximate the mean of a function \(F(Y)\) of a normally distributed random variable \(Y\) with mean \(\mu\) and variance \(\sigma^2\):

\[\int\limits_{-\infty}^\infty F(x) \frac{1}{\sqrt{2\pi\sigma}} \exp \left( -\frac{(y - \mu)^2}{2\sigma^2} \right) \mathrm{d}y\]

\[= \int\limits_{-\infty}^\infty F(\mu + \sqrt 2 \sigma x) \frac{1}{\sqrt\pi} \exp(-x^2) \mathrm{d}x\]

\[= \frac{1}{\sqrt\pi} \sum_{j=1}^N w_{j} F (\mu + \sqrt 2 \sigma x_{j})\]

Two-dimensional Gauss-Hermite integration

One-dimensional Gauss-Hermite integration can readily be extended to two dimensions. The following principal result in two dimensions is more or less given in Jäckel (2005) for the standard bivariate normal distribution \(\phi(x,y; \rho)\) with correlation parameter \(\rho\) :

\[\int\limits_{-\infty}^\infty \int\limits_{-\infty}^\infty F(x,y) \phi(x,y; \rho) \mathrm{d}x \mathrm{d}y \approx \frac{1}{\pi} \sum_{i=1}^N \sum_{j=1}^N w_{i} w_{j} F(\sqrt2 [ax_{i} + bx_{j}], \sqrt 2 [bx_{i} + ax_{j}])\]

in which

\[a = \frac{\sqrt{1 + \rho} + \sqrt{1 - \rho}}{2}\]

and

\[b = \frac{\sqrt{1 + \rho} - \sqrt{1 - \rho}}{2}\]

as given in Jäckel (2005) .

Jäckel (2005) discusses other Gauss-Hermite approximations to the two-dimensional integral, but found that the approximation given above generally gives the most accurate results. For the general bivariate normal distribution with means \((\mu_{x},\mu_{y} )\) and variances \((\sigma_{x}^2, \sigma_{y}^2 )\) the integral can be approximated by means of

\[\frac{1}{\pi} \sum_{i=1}^N \sum_{j=1}^N w_{i} w_{j} F(\mu_{x} + \sigma_{x} \sqrt2 [ax_{i} + bx_{j}], \mu_{y} + \sigma_{y} \sqrt 2 [bx_{i} + ax_{j}])\]

The product \(w_{i} w_{j}\) can be very small, especially when many quadrature points are used, thus wasting possibly precious calculation time. This can be remedied by pruning, i.e. by dropping combinations of \((i,j)\) with very small values of the product \(w_{i} w_{j}\).

Maximum likelihood for the LNN model with two-dimensional Gauss-Hermite integration

Denote non-consumption on day \(j\) for individual \(i\) as \(Y_{ij} = 0\). The conditional likelihood, i.e. given random effects \(b_{i}\) and \(v_{i}\), of a non-consumption on day \(j\) equals, with \(H()\) the inverse of the logit function

\(P(Y_{ij}=0 | b_{i},v_{i} ) = 1 – H(\lambda + v_{i})\).

The conditional likelihood of a positive intake \(Y_{ij}\) > 0 equals, with \(\phi\) the density of the normal distribution

\[f(Y_{ij}=y_{ij} | y_{ij}>0,b_{i},v_{i} )=H(\lambda + v_{i} ) \phi (y_{ij}-\mu-b_{i};0,\sigma_{w}^2 )\]

The conditional likelihood contribution for individual \(i\) is the product of the individual contributions for each day. The marginal likelihood contribution for individual \(i\) is obtained by integrating over the possible values of \(b_{i}\) and \(v_{i}\). Since the pair \((b_{i},v_{i} )\) follows a bivariate normal distribution, the likelihood contribution for individual \(i\) can be approximated by means of two-dimensional Gauss-Hermite integration. Individually based covariables, such as sex or age, imply that \(\mu_{i}\) and \(\lambda_{i}\) must be used instead of \(\mu\) and \(\lambda\). The likelihood must be optimized by means of some general optimization routine.