R/lav_integrate.R

In lavaan: Latent Variable Analysis

# routines for numerical intregration

# integrate (-infty to +infty) a product of univariate Gaussian densities
# with givens means (mus) and standard deviations (sds) (or variances, vars)
lav_integration_gaussian_product <- function(mus = NULL, sds = NULL, vars = NULL) {
  n <- length(mus)
  if (is.null(vars)) {
    vars <- sds^2
  }

  # variance product
  var.prod <- 1 / sum(1 / vars)

  # mean product
  mu.prod <- sum(mus / vars) * var.prod

  # normalization constant
  const <- 1 / sqrt((2 * pi)^(n - 1)) * sqrt(var.prod) * sqrt(1 / prod(vars)) * exp(-0.5 * (sum(mus^2 / vars) - mu.prod^2 / var.prod))

  const
}



# return Gauss-Hermite quadrature rule for given order (n)
# return list: x = nodes, w = quadrature weights
#

# As noted by Wilf (1962, chapter 2, ex 9), the nodes are given by
# the eigenvalues of the Jacobi matrix; weights are given by the squares of the
# first components of the (normalized) eigenvectors, multiplied by sqrt(pi)
#
# (This is NOT identical to Golub & Welsch, 1968: as they used a specific
#  method tailored for tridiagonal symmetric matrices)
#
# TODO: look at https://github.com/ajt60gaibb/FastGaussQuadrature.jl/blob/master/src/gausshermite.jl
# featuring the work of Ignace Bogaert (UGent)
#
# approximation of the integral of 'f(x) * exp(-x*x)' from -inf to +inf
# by sum( f(x_i) * w_i )
#
# CHECK: sum(w_i) should be always sqrt(pi) = 1.772454
lav_integration_gauss_hermite_xw <- function(n = 21L, revert = FALSE) {
  # force n to be an integer
  n <- as.integer(n)
  stopifnot(n > 0L)

  if (n == 1L) {
    x <- 0
    w <- sqrt(pi)
  } else {
    # construct symmetric, tridiagonal Jacobi matrix
    # diagonal = 0, -1/+1 diagonal is sqrt(1:(n-1)/2)
    u <- sqrt(seq.int(n - 1L) / 2) # upper diagonal of J
    Jn <- matrix(0, n, n)
    didx <- lav_matrix_diag_idx(n)
    Jn[(didx + 1)[-n]] <- u
    # Jn[(didx-1)[-1]] <- u # only lower matrix is used anyway

    # eigen decomposition
    # FIXME: use specialized function for tridiagonal symmetrix matrix
    ev <- eigen(Jn, symmetric = TRUE)
    x <- ev$values
    tmp <- ev$vectors[1L, ]
    w <- sqrt(pi) * tmp * tmp
  }

  # revert? (minus to plus)
  if (revert) {
    x <- -x
  }

  list(x = x, w = w)
}

# generate GH points + weights
lav_integration_gauss_hermite <- function(n = 21L,
                                          dnorm = FALSE,
                                          mean = 0, sd = 1,
                                          ndim = 1L,
                                          revert = TRUE,
                                          prune = FALSE) {
  XW <- lav_integration_gauss_hermite_xw(n = n, revert = revert)

  # dnorm kernel?
  if (dnorm) {
    # scale/shift x
    x <- XW$x * sqrt(2) * sd + mean

    # scale w
    w <- XW$w / sqrt(pi)
  } else {
    x <- XW$x
    w <- XW$w
  }

  if (ndim > 1L) {
    # cartesian product
    x <- as.matrix(expand.grid(rep(list(x), ndim), KEEP.OUT.ATTRS = FALSE))
    w <- as.matrix(expand.grid(rep(list(w), ndim), KEEP.OUT.ATTRS = FALSE))
    w <- apply(w, 1, prod)
  } else {
    x <- as.matrix(x)
    w <- as.matrix(w)
  }

  # prune?
  if (is.logical(prune) && prune) {
    # always divide by N=21
    lower.limit <- XW$w[1] * XW$w[floor((n + 1) / 2)] / 21
    keep.idx <- which(w > lower.limit)
    w <- w[keep.idx]
    x <- x[keep.idx, , drop = FALSE]
  } else if (is.numeric(prune) && prune > 0) {
    lower.limit <- quantile(w, probs = prune)
    keep.idx <- which(w > lower.limit)
    w <- w[keep.idx]
    x <- x[keep.idx, , drop = FALSE]
  }

  list(x = x, w = w)
}

# backwards compatibility
lav_integration_gauss_hermite_dnorm <- function(n = 21L, mean = 0, sd = 1,
                                                ndim = 1L,
                                                revert = TRUE,
                                                prune = FALSE) {
  lav_integration_gauss_hermite(
    n = n, dnorm = TRUE, mean = mean, sd = sd,
    ndim = ndim, revert = revert, prune = prune
  )
}

# plot 2-dim
# out <- lavaan:::lav_integration_gauss_hermite_dnorm(n = 20, ndim = 2)
# plot(out$x, cex = -10/log(out$w), col = "darkgrey", pch=19)

# integrand g(x) has the form g(x) = f(x) dnorm(x, m, s^2)
lav_integration_f_dnorm <- function(func = NULL, # often ly.prod
                                    dnorm.mean = 0, # dnorm mean
                                    dnorm.sd = 1, # dnorm sd
                                    XW = NULL, # GH points
                                    n = 21L, # number of nodes
                                    adaptive = FALSE, # adaptive?
                                    iterative = FALSE, # iterative?
                                    max.iter = 20L, # max iterations
                                    ...) { # optional args for 'f'

  # create GH rule
  if (is.null(XW)) {
    XW <- lav_integration_gauss_hermite_xw(n = n, revert = TRUE)
  }

  if (!adaptive) {
    w.star <- XW$w / sqrt(pi)
    x.star <- dnorm.sd * (sqrt(2) * XW$x) + dnorm.mean
    out <- sum(func(x.star, ...) * w.star)
  } else {
    # Naylor & Smith (1982, 1988)
    if (iterative) {
      mu.est <- 0
      sd.est <- 1

      for (i in 1:max.iter) {
        w.star <- sqrt(2) * sd.est * dnorm(sqrt(2) * sd.est * XW$x + mu.est, dnorm.mean, dnorm.sd) * exp(XW$x^2) * XW$w
        x.star <- sqrt(2) * sd.est * XW$x + mu.est
        LIK <- sum(func(x.star, ...) * w.star)

        # update mu
        mu.est <- sum(x.star * (func(x.star, ...) * w.star) / LIK)

        # update sd
        var.est <- sum(x.star^2 * (func(x.star, ...) * w.star) / LIK) - mu.est^2
        sd.est <- sqrt(var.est)

        if (lav_verbose()) {
          cat(
            "i = ", i, "LIK = ", LIK, "mu.est = ", mu.est,
            "sd.est = ", sd.est, "\n"
          )
        }
      }
      out <- LIK

      # Liu and Pierce (1994)
    } else {
      # integrand g(x) = func(x) * dnorm(x; m, s^2)
      log.g <- function(x, ...) {
        ## FIXME: should we take the log right away?
        log(func(x, ...) * dnorm(x, mean = dnorm.mean, sd = dnorm.sd))
      }
      # find mu hat and sd hat
      mu.est <- optimize(
        f = log.g, interval = c(-10, 10),
        maximum = TRUE, tol = .Machine$double.eps, ...
      )$maximum
      H <- as.numeric(numDeriv::hessian(func = log.g, x = mu.est, ...))
      sd.est <- sqrt(1 / -H)

      w.star <- sqrt(2) * sd.est * dnorm(sd.est * (sqrt(2) * XW$x) + mu.est, dnorm.mean, dnorm.sd) * exp(XW$x^2) * XW$w
      x.star <- sd.est * (sqrt(2) * XW$x) + mu.est

      out <- sum(func(x.star, ...) * w.star)
    }
  }

  out
}

# integrand g(z) has the form g(z) = f(sz+m) dnorm(z, 0, 1)
lav_integration_f_dnorm_z <- function(func = NULL, # often ly.prod
                                      f.mean = 0, # f mean
                                      f.sd = 1, # f sd
                                      XW = NULL, # GH points
                                      n = 21L, # number of nodes
                                      adaptive = FALSE, # adaptive?
                                      iterative = FALSE, # iterative?
                                      max.iter = 20L, # max iterations
                                      ...) { # optional args for 'f'

  # create GH rule
  if (is.null(XW)) {
    XW <- lav_integration_gauss_hermite_xw(n = n, revert = TRUE)
  }

  if (!adaptive) {
    w.star <- XW$w / sqrt(pi)
    x.star <- sqrt(2) * XW$x
    out <- sum(func(f.sd * x.star + f.mean, ...) * w.star)
  } else {
    # Naylor & Smith (1982, 1988)
    if (iterative) {
      mu.est <- 0
      sd.est <- 1

      for (i in 1:max.iter) {
        w.star <- sqrt(2) * sd.est * dnorm(sd.est * sqrt(2) * XW$x + mu.est, 0, 1) * exp(XW$x^2) * XW$w
        x.star <- sd.est * (sqrt(2) * XW$x) + mu.est
        LIK <- sum(func(f.sd * x.star + f.mean, ...) * w.star)

        # update mu
        mu.est <- sum(x.star * (func(f.sd * x.star + f.mean, ...) * w.star) / LIK)

        # update sd
        var.est <- sum(x.star^2 * (func(f.sd * x.star + f.mean, ...) * w.star) / LIK) - mu.est^2
        sd.est <- sqrt(var.est)

        if (lav_verbose()) {
          cat(
            "i = ", i, "LIK = ", LIK, "mu.est = ", mu.est,
            "sd.est = ", sd.est, "\n"
          )
        }
      }
      out <- LIK

      # Liu and Pierce (1994)
    } else {
      # integrand g(x) = func(x) * dnorm(x; m, s^2)
      log.gz <- function(x, ...) {
        ## FIXME: should we take the log right away?
        log(func(f.sd * x + f.mean, ...) * dnorm(x, mean = 0, sd = 1))
      }
      # find mu hat and sd hat
      mu.est <- optimize(
        f = log.gz, interval = c(-10, 10),
        maximum = TRUE, tol = .Machine$double.eps, ...
      )$maximum
      H <- as.numeric(numDeriv::hessian(func = log.gz, x = mu.est, ...))
      sd.est <- sqrt(1 / -H)

      w.star <- sqrt(2) * sd.est * dnorm(sd.est * (sqrt(2) * XW$x) + mu.est, 0, 1) * exp(XW$x^2) * XW$w
      x.star <- sd.est * (sqrt(2) * XW$x) + mu.est

      out <- sum(func(f.sd * x.star + f.mean, ...) * w.star)
    }
  }

  out
}