test-estimate-calculus.R
In lava: Latent Variable Models

# estimate calculus

set.seed(1)
# Generate mean-zero random IC matrix (for internal testing)
center_ic <- function(n, p = 1) {
  x <- matrix(rnorm(n * p), n, p)
  scale(x, center = TRUE, scale = FALSE)
}
a1 <- estimate(coef=1, IC=center_ic(10), id=1:10, labels="a1")
a2 <- estimate(coef=2, IC=center_ic(10), id=1:10, labels="a2")
a <- merge(a1, a2)

b1 <- estimate(coef=0.5, IC=center_ic(10), id=1:10, labels="b1")
b2 <- estimate(coef=0.9, IC=center_ic(10), id=1:10, labels="b2")
b <- merge(b1, b2)

test_that("+.estimate", {
  ## check addition: estimate + numeric with dim(estimate)=1
  e1 <- a1 + c(2,4,6)
  e2 <- estimate(a1, \(p) p + c(2,4,6))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## check addition: estimate + numeric with dim(estimate)>1
  e1 <- a + 2
  e2 <- estimate(a, \(p) p + 2)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a + c(2,4)
  e2 <- estimate(a, \(p) p + c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  # commutative
  e1 <- c(2,4) + a # 2.dim. estimate
  e2 <- estimate(a, \(p) p + c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- c(2,4) + a1 # 1.dim
  e2 <- estimate(a1, \(p) p + c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # addition with two estimate objects
  e1 <- a1 + a2
  e2 <- estimate(a, cbind(1, 1))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a + b # 2d. estimates
  e2 <- estimate(merge(a, b), rbind(c(1, 0, 1, 0), c(0, 1, 0, 1)))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

test_that("-.estimate", {
  ## check minus operator
  e1 <- -a
  e2 <- estimate(a, function(p) -p)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## estimate + numeric with dim(estimate)=1
  e1 <- a1 - c(2,4,6)
  e2 <- estimate(a1, \(p) p - c(2,4,6))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## estimate + numeric with dim(estimate)>1
  e1 <- a - 2
  e2 <- estimate(a, \(p) p - 2)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a - c(2,4)
  e2 <- estimate(a, \(p) p - c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # reverse order of numeric and estimate
  e1 <- c(2,4) - a # 2.dim. estimate
  e2 <- estimate(a, \(p) c(2,4) - p)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- c(2,4) - a1 # 1.dim
  e2 <- estimate(a1, \(p) c(2,4) - p)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # subtraction of two estimate objects
  e1 <- a1 - a2
  e2 <- estimate(a, cbind(1, -1))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a - b # 2d. estimates
  e2 <- estimate(merge(a, b), rbind(c(1, 0, -1, 0), c(0, 1, 0, -1)))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

test_that("*.estimate", {
  ## check multiplication: estimate + numeric with dim(estimate)=1
  e1 <- a1 * c(2,4,6)
  e2 <- estimate(a1, \(p) p * c(2,4,6))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## estimate + numeric with dim(estimate)>1
  e1 <- a * 2
  e2 <- estimate(a, \(p) p * 2)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a * c(2,4)
  e2 <- estimate(a, \(p) p * c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  # change order
  e1 <- c(2,4) * a # 2.dim. estimate
  e2 <- estimate(a, \(p) p * c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- c(2,4) * a1 # 1.dim
  e2 <- estimate(a1, \(p) p * c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # with two estimate objects
  e1 <- a1 * a2
  e2 <- estimate(a, prod)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a * b # 2d. estimates
  e2 <- estimate(merge(a, b), function(p) c(p[1]*p[3], p[2]*p[4]))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

test_that("/.estimate", {
  ## check division: estimate + numeric with dim(estimate)=1
  e1 <- a1 / c(2,4,6)
  e2 <- estimate(a1, \(p) p / c(2,4,6))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## estimate + numeric with dim(estimate)>1
  e1 <- a / 2
  e2 <- estimate(a, \(p) p / 2)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a / c(2,4)
  e2 <- estimate(a, \(p) p / c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  # numeric first
  e1 <- c(2,4) / a # 2.dim. estimate
  e2 <- estimate(a, \(p) c(2,4) / p)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- c(2,4) / a1 # 1.dim
  e2 <- estimate(a1, \(p) c(2,4) / p)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # with two estimate objects
  e1 <- a1 / a2
  e2 <- estimate(a, function(p) p[1] / p[2])
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a / b # 2d. estimates
  e2 <- estimate(merge(a, b), function(p) c(p[1]/p[3], p[2]/p[4]))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

test_that("^.estimate", {
  ## check power-func: estimate + numeric with dim(estimate)=1
  e1 <- a1^0.25
  e2 <- estimate(a1, function(p) p**0.25)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a1 / c(2,4,6)
  e2 <- estimate(a1, \(p) p / c(2,4,6))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## estimate + numeric with dim(estimate)>1
  e1 <- a ^ 3
  e2 <- estimate(a, \(p) p ^ 3)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a ^ c(2,4)
  e2 <- estimate(a, \(p) p ^ c(2,4))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  # numeric first
  e1 <- c(2,4) ^ a # 2.dim. estimate
  e2 <- estimate(a, \(p) c(2,4) ^ p)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- c(2,4) ^ a1 # 1.dim
  e2 <- estimate(a1, \(p) c(2,4) ^ p)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # addition with two estimate objects
  e1 <- a1 ^ a2
  e2 <- estimate(a, function(p) p[1] ^ p[2])
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a ^ b # 2d. estimates
  e2 <- estimate(merge(a, b), function(p) c(p[1]^p[3], p[2]^p[4]))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

test_that("math functions", {
  e0 <- c(a, b) # merge
  e <- merge(a, b)
  expect_equal(e, e0, ignore_attr = TRUE)

  # subset and merge
  ## e1 <- c(a=a["a1"], b=b["b1"])
  e1 <- merge(a["a1"], b["b1"], labels=c("a", "b"))
  expect_true(length(coef(e1)) == 2L)
  expect_equal(names(coef(e1)), c("a", "b"), ignore_attr = TRUE)

  # product
  e1 <- prod(e, labels="prod")
  e2 <- estimate(e, function(p) prod(p))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # sum
  e1 <- sum(e, labels="sum")
  e2 <- estimate(e, function(p) sum(p))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # inner-prod
  e1 <- a %*% b
  e2 <- estimate(e, function(p) sum(p[1:2]*p[3:4]))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- a %*% c(2, 3)
  e2 <- estimate(a, function(p) sum(p[1]*2 + p[2]*3))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # log
  e1 <- log(e["a1"])
  e2 <- estimate(e, function(p) log(p[1]))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- log(e)
  e2 <- estimate(e, function(p) log(p))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  e1 <- log1p(b)
  e2 <- estimate(b, function(p) log1p(p))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # exp
  e1 <- exp(e["a1"])
  e2 <- estimate(e, function(p) exp(p[1]))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  e1 <- exp(e)
  e2 <- estimate(e, function(p) exp(p))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  e1 <- expm1(b)
  e2 <- estimate(b, function(p) expm1(p))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # pow, sqrt
  e1 <- sqrt(e)
  e2 <- estimate(e, function(p) sqrt(p))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  # trigonometric
  for (f in c(cos, sin, tan, acos, asin, atan)) {
      e1 <- f(b)
      e2 <- estimate(b, function(p) f(p))
      testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
      testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  }
  e1 <- sin(asin(b))
  testthat::expect_equal(IC(e1), IC(b), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(b), ignore_attr = TRUE)

  # hyberbolic
  for (f in c(cosh, sinh, tanh)) {
      e1 <- f(b)
      e2 <- estimate(b, function(p) f(p))
      testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
      testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  }
  e1 <- acosh(cosh(b))
  testthat::expect_equal(IC(e1), IC(b), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(b), ignore_attr = TRUE)
  e1 <- asinh(sinh(b))
  testthat::expect_equal(IC(e1), IC(b), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(b), ignore_attr = TRUE)
  e1 <- atanh(tanh(b))
  testthat::expect_equal(IC(e1), IC(b), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(b), ignore_attr = TRUE)

})

test_that("custom functions", {
  e1 <- logit(expit(a))
  testthat::expect_equal(IC(e1), IC(a), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(a), ignore_attr = TRUE)

  e1 <- 2 * log(a["a2"]) * sqrt(b["b2"])
  e2 <- estimate(merge(a,b),
                 function(p) 2 * log(p["a2"]) * sqrt(p["b2"]))
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  e1 <- 1 / cos(a) - sqrt(a^2 + 1) + 1
  e2 <- estimate(a,
                 function(p) 1 / cos(p) - sqrt(p^2 + 1) + 1)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  e1 <- odds(b)
  e2 <- estimate(b, odds)
  testthat::expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  testthat::expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

## ---------------------------------------------------------------------------
## operator_grad shape correctness vs numerical jacobian (Issue C)
## ---------------------------------------------------------------------------
make_est <- function(coefs, seed) {
  set.seed(seed)
  k <- length(coefs)
  ic <- center_ic(20, k)
  nm <- if (is.null(names(coefs))) paste0("p", seq_along(coefs)) else names(coefs)
  names(coefs) <- nm
  estimate(coef = coefs, IC = ic, id = 1:20)
}

ops <- list(
  "+" = `+`,
  "-" = `-`,
  "*" = `*`,
  "/" = `/`,
  "^" = `^`
)

shapes <- list(
  c(1, 1),
  c(1, 3),
  c(3, 1),
  c(3, 3)
)

for (opname in names(ops)) {
  for (sh in shapes) {
    nx <- sh[1]; ny <- sh[2]
    label <- sprintf("operator_grad: %s with shape (%d,%d)", opname, nx, ny)
    test_that(label, {
      ## Use strictly positive coefs to keep ^ and / well-defined.
      a <- make_est(setNames(seq_len(nx) + 1.0, paste0("a", seq_len(nx))),
                    seed = 100 + nx)
      b <- make_est(setNames(seq_len(ny) + 0.5, paste0("b", seq_len(ny))),
                    seed = 200 + ny)
      op_est <- ops[[opname]]
      op_num <- switch(opname,
                       "+" = function(p) p[seq_len(nx)] + p[nx + seq_len(ny)],
                       "-" = function(p) p[seq_len(nx)] - p[nx + seq_len(ny)],
                       "*" = function(p) p[seq_len(nx)] * p[nx + seq_len(ny)],
                       "/" = function(p) p[seq_len(nx)] / p[nx + seq_len(ny)],
                       "^" = function(p) p[seq_len(nx)] ^ p[nx + seq_len(ny)])
      ## Analytic via the operator
      res <- op_est(a, b)
      res_coef <- coef(res)
      ## Numerical jacobian
      pp <- c(coef(a), coef(b))
      Jnum <- numDeriv::jacobian(op_num, pp)
      ## Reconstruct analytic D from IC: IC_res = IC_ab %*% t(D)
      ## We instead validate by comparing variance: an incorrect D shape
      ## would produce wrong vcov dims or wrong values.
      ## Build merged IC manually:
      ic_ab <- cbind(IC(a), IC(b))
      V_expected <- Jnum %*% var_ic(ic_ab) %*% t(Jnum)
      V_actual   <- vcov(res)
      expect_equal(dim(V_actual), dim(V_expected))
      expect_equal(V_actual, V_expected, tolerance = 1e-5, ignore_attr = TRUE)
      expect_equal(length(res_coef), max(nx, ny))
    })
  }
}

## ---------------------------------------------------------------------------
## Normal distribution: normcdf / normpdf / norminv
## ---------------------------------------------------------------------------

test_that("normcdf/normpdf/norminv default methods", {
  x <- c(-1, 0, 1)
  expect_equal(normcdf(x, 0.5, 2), pnorm(x, 0.5, 2))
  expect_equal(normcdf(x, 0.5, 2, lower.tail = FALSE, log.p = TRUE),
               pnorm(x, 0.5, 2, lower.tail = FALSE, log.p = TRUE))
  expect_equal(normpdf(x, 0.5, 2), dnorm(x, 0.5, 2))
  expect_equal(normpdf(x, 0.5, 2, log = TRUE), dnorm(x, 0.5, 2, log = TRUE))
  pr <- c(0.1, 0.5, 0.9)
  expect_equal(norminv(pr, 0.5, 2), qnorm(pr, 0.5, 2))
  expect_equal(norminv(log(pr), 0.5, 2, lower.tail = FALSE, log.p = TRUE),
               qnorm(log(pr), 0.5, 2, lower.tail = FALSE, log.p = TRUE))
})

test_that("normcdf.estimate", {
  e1 <- normcdf(a)
  e2 <- estimate(a, function(p) pnorm(p))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## lower.tail / log.p
  e1 <- normcdf(a, lower.tail = FALSE, log.p = TRUE)
  e2 <- estimate(a, function(p) pnorm(p, lower.tail = FALSE, log.p = TRUE))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## constant mean/sd
  e1 <- normcdf(a, mean = 2, sd = 3)
  e2 <- estimate(a, function(p) pnorm(p, 2, 3))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## mean and sd given as estimate objects
  s <- estimate(coef = c(s = 1.5), IC = center_ic(10), id = 1:10)
  e1 <- normcdf(a, mean = b, sd = s)
  e2 <- estimate(merge(a, b, s),
                 function(p) pnorm(p[1:2], p[3:4], p[5]))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## numeric first argument, estimate mean/sd (dispatch via default method)
  e1 <- normcdf(1.5, mean = b, sd = s)
  e2 <- estimate(merge(b, s), function(p) pnorm(1.5, p[1:2], p[3]))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## scalar estimate recycled against longer numeric argument
  e1 <- normcdf(c(1, 2, 3), mean = s)
  e2 <- estimate(s, function(p) pnorm(c(1, 2, 3), p))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  ## labels inherited
  expect_equal(names(coef(normcdf(a))), names(coef(a)), ignore_attr = TRUE)
})

test_that("normpdf.estimate", {
  s <- estimate(coef = c(s = 1.5), IC = center_ic(10), id = 1:10)
  for (lg in c(FALSE, TRUE)) {
    e1 <- normpdf(a, log = lg)
    e2 <- estimate(a, function(p) dnorm(p, log = lg))
    expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
    expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

    e1 <- normpdf(a, mean = b, sd = s, log = lg)
    e2 <- estimate(merge(a, b, s),
                   function(p) dnorm(p[1:2], p[3:4], p[5], log = lg))
    expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
    expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
  }
  ## exp(log-density) == density
  e1 <- exp(normpdf(a, log = TRUE))
  e2 <- normpdf(a)
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

test_that("norminv.estimate", {
  s <- estimate(coef = c(s = 1.5), IC = center_ic(10), id = 1:10)
  pr <- estimate(coef = c(p1 = 0.3, p2 = 0.75),
                 IC = center_ic(10, 2) / 20, id = 1:10)
  e1 <- norminv(pr)
  e2 <- estimate(pr, function(p) qnorm(p))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  e1 <- norminv(pr, lower.tail = FALSE)
  e2 <- estimate(pr, function(p) qnorm(p, lower.tail = FALSE))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  lp <- log(pr)
  e1 <- norminv(lp, lower.tail = FALSE, log.p = TRUE)
  e2 <- estimate(lp, function(p) qnorm(p, lower.tail = FALSE, log.p = TRUE))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  e1 <- norminv(pr, mean = b, sd = s)
  e2 <- estimate(merge(pr, b, s), function(p) qnorm(p[1:2], p[3:4], p[5]))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)

  e1 <- norminv(0.975, mean = b, sd = s)
  e2 <- estimate(merge(b, s), function(p) qnorm(0.975, p[1:2], p[3]))
  expect_equal(IC(e1), IC(e2), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(e2), ignore_attr = TRUE)
})

test_that("normcdf/norminv are inverse transformations", {
  pr <- estimate(coef = c(p1 = 0.3, p2 = 0.75),
                 IC = center_ic(10, 2) / 20, id = 1:10)
  e1 <- norminv(normcdf(a))
  expect_equal(IC(e1), IC(a), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(a), ignore_attr = TRUE)
  e1 <- normcdf(norminv(pr))
  expect_equal(IC(e1), IC(pr), ignore_attr = TRUE)
  expect_equal(coef(e1), coef(pr), ignore_attr = TRUE)
})

test_that("normcdf/normpdf/norminv work without influence function", {
  e <- estimate(coef = c(a = 1, b = 2), vcov = diag(2) * 0.1)
  res <- normcdf(e)
  expect_equal(coef(res), pnorm(coef(e)), ignore_attr = TRUE)
  d <- diag(dnorm(coef(e)))
  expect_equal(vcov(res), d %*% (diag(2) * 0.1) %*% d, ignore_attr = TRUE)
})

test_that("normcdf.estimate incompatible lengths", {
  expect_error(normcdf(a, mean = c(1, 2, 3)))
})

## ---------------------------------------------------------------------------
## diff
## ---------------------------------------------------------------------------

test_that("diff.estimate", {
  e <- c(a, b)
  f <- diff(e)
  expect_true(length(coef(f))==3L)
  expect_equal(unname(coef(f)), unname(diff(coef(e))))
  lbl <- gsub(" ", "", names(coef(f)))
  pname <- paste0("[",names(coef(e)),"]")
  fname <- character(length(pname)-1)
  for (i in seq_len(length(pname)-1)) {
    fname[i] <- paste0("-",pname[i],"+",pname[i+1])
  }
  expect_equal(lbl, fname)
  expect_equal(IC(f)[,1], -IC(a)[,1]+IC(a)[,2])
  expect_equal(IC(f)[,2], -IC(a)[,2]+IC(b)[,1])
  expect_equal(IC(f)[,3], -IC(b)[,1]+IC(b)[,2])
})

test_that("pairwise_diff", {
  e <- c(a, b)
  f <- pairwise_diff(e)
  expect_equal(length(coef(f)), choose(length(coef(e)), 2))
})

## ---------------------------------------------------------------------------
## gamma
## ---------------------------------------------------------------------------

test_that("gamma.estimate", {
  a1 <- a + c(2,8)
  e1 <- gamma(a1)
  e2 <- exp(lgamma(a1))
  expect_equal(e1, e2)
  expect_equal(coef(e1), gamma(coef(a1)))
})

## ---------------------------------------------------------------------------
## subset methods
## ---------------------------------------------------------------------------

test_that("subset estimate methods", {
  p <- 5
  cc <- seq_len(p)
  e <- estimate(
    coef = cc,
    IC = matrix(c(cc, -cc), nrow=2, ncol=p, byrow=TRUE),
    labels = rev(letters[seq_len(p)])
  )

  f <- head(e, 2)
  expect_equal(coef(e)[1:2], coef(f))
  expect_equal(f, subset(e, 1:2))
  expect_equal(f, subset(e, names(coef(e))[1:2]))
  expect_equal(f, e[1:2])

  f <- tail(e, 2)
  idx <- seq(p-1, p)
  expect_equal(coef(e)[idx], coef(f))
  expect_equal(f, subset(e, idx))
  expect_equal(f, e[idx])

  f <- rev(e)
  idx <- seq(p, 1)
  expect_equal(f, subset(e, idx))

  f <- sort(e)
  expect_equal(names(coef(f)), sort(names(coef(e))))

  f <- sort(e, decreasing=TRUE)
  expect_equal(sort(coef(e)), coef(f))

})