R Under development (unstable) (2024-09-15 r87152 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
>  library(aster2)
Loading required package: Matrix
This is beta software.
Unless you need to do aster models with dependence groups,
    use package "aster" instead.
See help(aster2-package) for differences from package "aster"
    and examples.

> 
>  set.seed(42)
> 
>  epsilon <- 1e-6
> 
>  ##### Bernoulli #####
> 
>  theta <- seq(-3, 3, 0.1)
> 
>  cumfun <- function(theta)
+      ifelse(theta > 0, theta + log1p(exp(- theta)), log1p(exp(theta)))
> 
>  moofun <- function(theta)
+      ifelse(theta > 0, 1 / (1 + exp(- theta)), exp(theta) / (1 + exp(theta)))
> 
>  voofun <- function(theta) {
+      pp <- moofun(theta)
+      qq <- moofun(- theta)
+      pp * qq
+  }
> 
>  thoofun <- function(theta) {
+      pp <- moofun(theta)
+      qq <- moofun(- theta)
+      - pp * qq * tanh(theta / 2)
+  }
> 
>  zeroth <- mapply(function(theta) cumulant(theta, fam.bernoulli())$zeroth,
+      theta)
> 
>  all.equal(zeroth, cumfun(theta))
[1] TRUE
> 
>  first <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      deriv = 1)$first, theta)
> 
>  all.equal(first, moofun(theta))
[1] TRUE
> 
>  second <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      deriv = 2)$second, theta)
> 
>  all.equal(second, voofun(theta))
[1] TRUE
> 
>  third <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      deriv = 3)$third, theta)
> 
>  all.equal(third, thoofun(theta))
[1] TRUE
> 
>  grad <- moofun(theta)
>  fdgrad <- (cumfun(theta + epsilon) - cumfun(theta - epsilon)) / (2 * epsilon)
>  all.equal(grad, fdgrad)
[1] TRUE
> 
>  hess <- voofun(theta)
>  fdhess <- (moofun(theta + epsilon) - moofun(theta - epsilon)) / (2 * epsilon)
>  all.equal(hess, fdhess)
[1] TRUE
> 
>  thss <- thoofun(theta)
>  fdth <- (voofun(theta + epsilon) - voofun(theta - epsilon)) / (2 * epsilon)
>  all.equal(thss, fdth)
[1] TRUE
> 
>  foo <- mapply(function(xi) link(xi, fam.bernoulli(), deriv = 0)$zeroth, first)
> 
>  all.equal(theta, foo)
[1] TRUE
> 
>  foo <- mapply(function(xi) link(xi, fam.bernoulli(), deriv = 1)$first, first)
> 
>  all.equal(second, 1 / foo)
[1] TRUE
> 
>  ##### Poisson #####
> 
>  zeroth <- mapply(function(theta) cumulant(theta, fam.poisson())$zeroth, theta)
> 
>  all.equal(zeroth, exp(theta))
[1] TRUE
> 
>  first <- mapply(function(theta) cumulant(theta, fam.poisson(),
+      deriv = 1)$first, theta)
> 
>  all.equal(first, exp(theta))
[1] TRUE
> 
>  second <- mapply(function(theta) cumulant(theta, fam.poisson(),
+      deriv = 2)$second, theta)
> 
>  all.equal(second, exp(theta))
[1] TRUE
> 
>  third <- mapply(function(theta) cumulant(theta, fam.poisson(),
+      deriv = 3)$third, theta)
> 
>  all.equal(third, exp(theta))
[1] TRUE
> 
>  foo <- mapply(function(xi) link(xi, fam.poisson(), deriv = 0)$zeroth, first)
> 
>  all.equal(theta, foo)
[1] TRUE
> 
>  foo <- mapply(function(xi) link(xi, fam.poisson(), deriv = 1)$first, first)
> 
>  all.equal(second, 1 / foo)
[1] TRUE
> 
>  ##### zero-truncated Poisson #####
> 
>  cumfun <- function(theta) {
+      m <- exp(theta)
+      result <- m + log1p(- exp(- m))
+      bar <- m / 2 * (1 + m / 3 * (1 + m / 4 * (1 + m / 5 * (1 + m / 6 *
+          (1 + m / 7 * (1 + m / 8))))))
+      inies <- theta < (- 4)
+      result[inies] <- (theta + log1p(bar))[inies]
+      return(result)
+  }
> 
>  moofun <- function(theta) {
+      m <- exp(theta)
+      result <- m / (1 - exp(- m))
+      bar <- m / 2 * (1 + m / 3 * (1 + m / 4 * (1 + m / 5 * (1 + m / 6 *
+          (1 + m / 7 * (1 + m / 8))))))
+      inies <- theta < (- 4)
+      result[inies] <- (m + 1 / (1 + bar))[inies]
+      return(result)
+  }
> 
>  voofun <- function(theta) {
+      m <- exp(theta)
+      mu <- moofun(theta)
+      return(mu * (1 + m - mu))
+  }
> 
>  thoofun <- function(theta) {
+      m <- exp(theta)
+      mu <- moofun(theta)
+      return(mu * ((1 + m - mu) * (1 + m - 2 * mu) + m))
+  }
> 
>  theta <- seq(-12, 2, 0.1)
> 
>  zeroth <- mapply(function(theta)
+      cumulant(theta, fam.zero.truncated.poisson())$zeroth, theta)
> 
>  all.equal(zeroth, cumfun(theta))
[1] TRUE
> 
>  first <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
+      deriv = 1)$first, theta)
> 
>  all.equal(first, moofun(theta))
[1] TRUE
> 
>  second <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
+      deriv = 2)$second, theta)
> 
>  all.equal(second, voofun(theta))
[1] TRUE
> 
>  third <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
+      deriv = 3)$third, theta)
> 
>  all.equal(third, thoofun(theta))
[1] TRUE
> 
>  grad <- moofun(theta)
>  fdgrad <- (cumfun(theta + epsilon) - cumfun(theta - epsilon)) / (2 * epsilon)
>  all.equal(grad, fdgrad)
[1] TRUE
> 
>  hess <- voofun(theta)
>  fdhess <- (moofun(theta + epsilon) - moofun(theta - epsilon)) / (2 * epsilon)
>  all.equal(hess, fdhess)
[1] TRUE
> 
>  thss <- thoofun(theta)
>  fdth <- (voofun(theta + epsilon) - voofun(theta - epsilon)) / (2 * epsilon)
>  all.equal(thss, fdth)
[1] TRUE
> 
>  foo <- mapply(function(xi) link(xi, fam.zero.truncated.poisson(),
+      deriv = 0)$zeroth, first)
> 
>  all.equal(theta, foo)
[1] TRUE
> 
>  foo <- mapply(function(xi) link(xi, fam.zero.truncated.poisson(),
+      deriv = 1)$first, first)
> 
>  all.equal(second, 1 / foo)
[1] TRUE
> 
>  ##### Multinomial #####
> 
>  d <- 4
> 
>  cumfun <- function(theta) {
+      stopifnot(length(theta) == d)
+      theta.max = max(theta)
+      j <- seq(along = theta)[theta == theta.max]
+      j <- j[1]
+      return(theta.max + log1p(sum(exp(theta - theta.max)[-j])))
+  }
> 
>  moofun <- function(theta) {
+      stopifnot(length(theta) == d)
+      theta.max = max(theta)
+      return(exp(theta - theta.max) / sum(exp(theta - theta.max)))
+  }
> 
>  voofun <- function(theta) {
+      mu <- moofun(theta)
+      diag(mu) - outer(mu, mu)
+  }
> 
>  thoofun <- function(theta) {
+      mu <- moofun(theta)
+      result <- 2 * outer(mu, outer(mu, mu))
+      for (i in 1:d) {
+          for (j in 1:d) {
+              if (i != j) {
+                  qux <- mu[i] * mu[j]
+                  result[i, i, j] <- result[i, i, j] - qux
+                  result[i, j, i] <- result[i, j, i] - qux
+                  result[j, i, i] <- result[j, i, i] - qux
+              }
+          }
+          qux <- mu[i]
+          result[i, i, i] <- qux * (1 - qux) * (1 - 2 * qux)
+      }
+      return(result)
+  }
> 
>  theta <- matrix(rnorm(d * 25), ncol = d)
> 
>  itself <- apply(theta, 1, cumfun)
> 
>  grad <- t(apply(theta, 1, moofun))
>  all.equal(rowSums(grad), rep(1, nrow(grad)))
[1] TRUE
> 
>  fred <- function(theta) {
+      stopifnot(length(theta) == d)
+      result <- double(d)
+      for (i in 1:d) {
+          theta.plus <- theta
+          theta.plus[i] <- theta.plus[i] + epsilon
+          theta.minus <- theta
+          theta.minus[i] <- theta.minus[i] - epsilon
+          result[i] <- (cumfun(theta.plus) - cumfun(theta.minus)) / (2 * epsilon)
+      }
+      return(result)
+  }
>  fdgrad <- t(apply(theta, 1, fred))
>  all.equal(grad, fdgrad)
[1] TRUE
> 
>  hess <- apply(theta, 1, voofun)
>  hess <- array(as.vector(hess), dim = c(d, d, ncol(hess)))
>  dim(hess)
[1]  4  4 25
> 
>  fred <- function(theta) {
+      stopifnot(length(theta) == d)
+      result <- matrix(NA, d, d)
+      for (i in 1:d) {
+          theta.plus <- theta
+          theta.plus[i] <- theta.plus[i] + epsilon
+          theta.minus <- theta
+          theta.minus[i] <- theta.minus[i] - epsilon
+          result[i, ] <- (moofun(theta.plus) - moofun(theta.minus)) /
+              (2 * epsilon)
+      }
+      return(result)
+  }
> 
>  fdhess <- apply(theta, 1, fred)
>  fdhess <- array(as.vector(fdhess), dim = c(d, d, ncol(fdhess)))
>  all.equal(hess, fdhess)
[1] TRUE
> 
>  thss <- apply(theta, 1, thoofun)
>  thss <- array(as.vector(thss), dim = c(d, d, d, ncol(thss)))
>  dim(thss)
[1]  4  4  4 25
> 
>  fred <- function(theta) {
+      stopifnot(length(theta) == d)
+      result <- array(NA, dim = c(d, d, d))
+      for (i in 1:d) {
+          theta.plus <- theta
+          theta.plus[i] <- theta.plus[i] + epsilon
+          theta.minus <- theta
+          theta.minus[i] <- theta.minus[i] - epsilon
+          result[i, , ] <- (voofun(theta.plus) - voofun(theta.minus)) /
+              (2 * epsilon)
+      }
+      return(result)
+  }
>  fdthss <- apply(theta, 1, fred)
>  fdthss <- array(as.vector(fdthss), dim = c(d, d, d, ncol(fdthss)))
>  all.equal(thss, fdthss)
[1] TRUE
> 
>  cumfun <- function(theta) cumulant(theta, fam.multinomial(d))$zeroth
>  moofun <- function(theta) cumulant(theta, fam.multinomial(d), deriv = 1)$first
>  voofun <- function(theta) cumulant(theta, fam.multinomial(d), deriv = 2)$second
>  thoofun <- function(theta) cumulant(theta, fam.multinomial(d), deriv = 3)$third
> 
>  all.equal(itself, apply(theta, 1, cumfun))
[1] TRUE
>  all.equal(grad, t(apply(theta, 1, moofun)))
[1] TRUE
>  chess <- apply(theta, 1, voofun)
>  chess <- array(as.vector(chess), dim = c(d, d, ncol(chess)))
>  all.equal(hess, chess)
[1] TRUE
>  cthss <- apply(theta, 1, thoofun)
>  cthss <- array(as.vector(cthss), dim = c(d, d, d, ncol(cthss)))
>  all.equal(cthss, thss)
[1] TRUE
> 
>  foo <- apply(grad, 1, function(xi) link(xi, fam.multinomial(d),
+      deriv = 0)$zeroth)
>  bar <- apply(t(foo), 1, function(theta) cumulant(theta, fam.multinomial(d),
+      deriv = 1)$first)
>  all.equal(t(bar), grad)
[1] TRUE
> 
>  baz <- apply(grad, 1, function(xi) link(xi, fam.multinomial(d),
+      deriv = 1)$first)
>  baz <- array(as.vector(baz), dim = c(d, d, nrow(theta)))
>  fdjackfun <- function(xi) {
+      result <- matrix(NA, d, d)
+      for (i in 1:d) {
+          xi.plus <- xi
+          xi.plus[i] <- xi[i] + epsilon
+          xi.minus <- xi
+          xi.minus[i] <- xi[i] - epsilon
+          result[ , i] <- (link(xi.plus, fam.multinomial(d))$zeroth -
+              link(xi.minus, fam.multinomial(d))$zeroth) / (2 * epsilon)
+      }
+      return(result)
+  }
>  qux <- apply(grad, 1, fdjackfun)
>  qux <- array(as.vector(qux), dim = c(d, d, nrow(theta)))
>  all.equal(baz, qux)
[1] TRUE
> 
>  ##### Normal location-scale #####
> 
>  d <- 2
> 
>  cumfun <- function(theta) {
+      stopifnot(length(theta) == d)
+      stopifnot(theta[2] < 0)
+      - theta[1]^2 / (4 * theta[2]) + (1 / 2) * log(- 1 / (2 * theta[2]))
+  }
> 
>  moofun <- function(theta) {
+      stopifnot(length(theta) == d)
+      stopifnot(theta[2] < 0)
+      r1 <- (- theta[1] / (2 * theta[2]))
+      r2 <- theta[1]^2 / (4 * theta[2]^2) - 1 / (2 * theta[2])
+      return(c(r1, r2))
+  }
> 
>  voofun <- function(theta) {
+      stopifnot(length(theta) == d)
+      stopifnot(theta[2] < 0)
+      r11 <- (- 1 / (2 * theta[2]))
+      r12 <- theta[1] / (2 * theta[2]^2)
+      r22 <- (- theta[1]^2 / (2 * theta[2]^3)) + 1 / (2 * theta[2]^2)
+      return(matrix(c(r11, r12, r12, r22), 2, 2))
+  }
> 
>  thoofun <- function(theta) {
+      stopifnot(length(theta) == d)
+      stopifnot(theta[2] < 0)
+      r111 <- 0
+      r112 <- 1 / (2 * theta[2]^2)
+      r122 <- (- theta[1] / theta[2]^3)
+      r222 <- 3 * theta[1]^2 / (2 * theta[2]^4) - 1 / theta[2]^3
+      return(array(c(r111, r112, r112, r122, r112, r122, r122, r222), dim = c(2, 2, 2)))
+  }
> 
>  theta <- cbind(rnorm(25), - rexp(25))
> 
>  itself <- apply(theta, 1, cumfun)
> 
>  grad <- t(apply(theta, 1, moofun))
> 
>  fred <- function(theta) {
+      stopifnot(length(theta) == d)
+      result <- double(d)
+      for (i in 1:d) {
+          theta.plus <- theta
+          theta.plus[i] <- theta.plus[i] + epsilon
+          theta.minus <- theta
+          theta.minus[i] <- theta.minus[i] - epsilon
+          result[i] <- (cumfun(theta.plus) - cumfun(theta.minus)) / (2 * epsilon)
+      }
+      return(result)
+  }
>  fdgrad <- t(apply(theta, 1, fred))
>  all.equal(grad, fdgrad)
[1] TRUE
> 
>  hess <- apply(theta, 1, voofun)
>  hess <- array(as.vector(hess), dim = c(d, d, ncol(hess)))
>  dim(hess)
[1]  2  2 25
> 
>  fred <- function(theta) {
+      stopifnot(length(theta) == d)
+      result <- matrix(NA, d, d)
+      for (i in 1:d) {
+          theta.plus <- theta
+          theta.plus[i] <- theta.plus[i] + epsilon
+          theta.minus <- theta
+          theta.minus[i] <- theta.minus[i] - epsilon
+          result[i, ] <- (moofun(theta.plus) - moofun(theta.minus)) /
+              (2 * epsilon)
+      }
+      return(result)
+  }
> 
>  fdhess <- apply(theta, 1, fred)
>  fdhess <- array(as.vector(fdhess), dim = c(d, d, ncol(fdhess)))
>  all.equal(hess, fdhess)
[1] TRUE
> 
>  thss <- apply(theta, 1, thoofun)
>  thss <- array(as.vector(thss), dim = c(d, d, d, ncol(thss)))
>  dim(thss)
[1]  2  2  2 25
> 
>  fred <- function(theta) {
+      stopifnot(length(theta) == d)
+      result <- array(NA, dim = c(d, d, d))
+      for (i in 1:d) {
+          theta.plus <- theta
+          theta.plus[i] <- theta.plus[i] + epsilon
+          theta.minus <- theta
+          theta.minus[i] <- theta.minus[i] - epsilon
+          result[i, , ] <- (voofun(theta.plus) - voofun(theta.minus)) /
+              (2 * epsilon)
+      }
+      return(result)
+  }
> 
>  fdthss <- apply(theta, 1, fred)
>  fdthss <- array(as.vector(fdthss), dim = c(d, d, d, ncol(fdthss)))
>  all.equal(thss, fdthss)
[1] TRUE
> 
>  cumfun <- function(theta) cumulant(theta, fam.normal.location.scale())$zeroth
>  moofun <- function(theta) cumulant(theta, fam.normal.location.scale(),
+      deriv = 1)$first
>  voofun <- function(theta) cumulant(theta, fam.normal.location.scale(),
+      deriv = 2)$second
>  thoofun <- function(theta) cumulant(theta, fam.normal.location.scale(),
+      deriv = 3)$third
> 
>  all.equal(itself, apply(theta, 1, cumfun))
[1] TRUE
>  all.equal(grad, t(apply(theta, 1, moofun)))
[1] TRUE
>  chess <- apply(theta, 1, voofun)
>  chess <- array(as.vector(chess), dim = c(d, d, ncol(chess)))
>  all.equal(hess, chess)
[1] TRUE
>  cthss <- apply(theta, 1, thoofun)
>  cthss <- array(as.vector(cthss), dim = c(d, d, d, ncol(cthss)))
>  all.equal(cthss, thss)
[1] TRUE
> 
>  foo <- apply(grad, 1, function(xi) link(xi, fam.normal.location.scale(),
+      deriv = 0)$zeroth)
>  foo <- t(foo)
>  all.equal(theta, foo)
[1] TRUE
> 
>  foo <- apply(grad, 1, function(xi) link(xi, fam.normal.location.scale(),
+      deriv = 1)$first)
>  foo <- array(as.vector(foo), dim = c(2, 2, nrow(theta)))
>  is.ok <- as.logical(rep(0, nrow(theta)))
>  for (i in 1:nrow(theta))
+      is.ok[i] <- all.equal(foo[ , , i] %*% hess[ , , i], diag(2))
>  all(is.ok)
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   1.87    0.20    2.06