Gauss-Hermite Integration

One-dimensional Gauss-Hermite integration

Gauss-Hermite integration approximates a specific integral as follows

$$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}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$:

$$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}F(x)\frac{1}{\sqrt{2\pi \sigma }}\exp (-\frac{(y-\mu {)}^2}{2{\sigma }^2})\mathrm{d}y$$

$$=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}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 $\varphi (x,y;\rho )$ with correlation parameter $\rho$ :

$$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}F(x,y)\varphi (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(\sqrt{2}[a{x}_i+b{x}_j],\sqrt{2}[b{x}_i+a{x}_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\sqrt{2}[a{x}_i+b{x}_j],{\mu }_y+{\sigma }_y\sqrt{2}[b{x}_i+a{x}_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 $\varphi$ the density of the normal distribution

$$f(Y_{ij}=y_{ij}|y_{ij}>0,b_i,v_i)=H(\lambda +v_i)\varphi (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.