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.

> 
>  ##### Bernoulli, lower limit #####
> 
>  theta <- seq(-3, 3, 0.1)
> 
>  cumfun <- function(theta) rep(0, length(theta))
>  moofun <- cumfun
>  voofun <- cumfun
>  thoofun <- cumfun
> 
>  zeroth <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = -1)$zeroth, theta)
>  all.equal(zeroth, cumfun(theta))
[1] TRUE
> 
>  first <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = -1, deriv = 1)$first, theta)
>  all.equal(first, moofun(theta))
[1] TRUE
> 
>  second <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = - 1, deriv = 2)$second, theta)
>  all.equal(second, voofun(theta))
[1] TRUE
> 
>  third <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = - 1, deriv = 3)$third, theta)
>  all.equal(third, thoofun(theta))
[1] TRUE
> 
>  foo <- link(0, fam.bernoulli(), delta = -1, deriv = 1)
>  # does not matter what zeroth is, theta is not identifiable for
>  # degenerate model
>  is.finite(foo$zeroth)
[1] TRUE
>  all.equal(0, foo$first)
[1] TRUE
> 
>  ##### Bernoulli, upper limit #####
> 
>  cumfun <- function(theta) theta
>  moofun <- function(theta) rep(1, length(theta))
>  voofun <- function(theta) rep(0, length(theta))
>  thoofun <- voofun
> 
>  zeroth <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = 1)$zeroth, theta)
>  all.equal(zeroth, cumfun(theta))
[1] TRUE
> 
>  first <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = 1, deriv = 1)$first, theta)
>  all.equal(first, moofun(theta))
[1] TRUE
> 
>  second <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = 1, deriv = 2)$second, theta)
>  all.equal(second, voofun(theta))
[1] TRUE
> 
>  third <- mapply(function(theta) cumulant(theta, fam.bernoulli(),
+      delta = 1, deriv = 3)$third, theta)
>  all.equal(third, thoofun(theta))
[1] TRUE
> 
>  foo <- link(1, fam.bernoulli(), delta = 1, deriv = 1)
>  # does not matter what zeroth is, theta is not identifiable for
>  # degenerate model
>  is.finite(foo$zeroth)
[1] TRUE
>  all.equal(0, foo$first)
[1] TRUE
> 
>  ##### Poisson, lower limit #####
> 
>  theta <- seq(-3, 3, 0.1)
> 
>  cumfun <- function(theta) rep(0, length(theta))
>  moofun <- cumfun
>  voofun <- cumfun
>  thoofun <- cumfun
> 
>  zeroth <- mapply(function(theta) cumulant(theta, fam.poisson(),
+      delta = -1)$zeroth, theta)
>  all.equal(zeroth, cumfun(theta))
[1] TRUE
> 
>  first <- mapply(function(theta) cumulant(theta, fam.poisson(),
+      delta = -1, deriv = 1)$first, theta)
>  all.equal(first, moofun(theta))
[1] TRUE
> 
>  second <- mapply(function(theta) cumulant(theta, fam.poisson(),
+      delta = - 1, deriv = 2)$second, theta)
>  all.equal(second, voofun(theta))
[1] TRUE
> 
>  third <- mapply(function(theta) cumulant(theta, fam.poisson(),
+      delta = - 1, deriv = 3)$third, theta)
>  all.equal(third, thoofun(theta))
[1] TRUE
> 
>  foo <- link(0, fam.poisson(), delta = -1, deriv = 1)
>  # does not matter what zeroth is, theta is not identifiable for
>  # degenerate model
>  is.finite(foo$zeroth)
[1] TRUE
>  all.equal(0, foo$first)
[1] TRUE
> 
>  ##### zero-truncated Poisson, lower limit #####
> 
>  cumfun <- function(theta) theta
>  moofun <- function(theta) rep(1, length(theta))
>  voofun <- function(theta) rep(0, length(theta))
>  thoofun <- voofun
> 
>  zeroth <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
+      delta = -1)$zeroth, theta)
>  all.equal(zeroth, cumfun(theta))
[1] TRUE
> 
>  first <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
+      delta = -1, deriv = 1)$first, theta)
>  all.equal(first, moofun(theta))
[1] TRUE
> 
>  second <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
+      delta = -1, deriv = 2)$second, theta)
>  all.equal(second, voofun(theta))
[1] TRUE
> 
>  third <- mapply(function(theta) cumulant(theta, fam.zero.truncated.poisson(),
+      delta = -1, deriv = 3)$third, theta)
>  all.equal(third, thoofun(theta))
[1] TRUE
> 
>  foo <- link(1, fam.zero.truncated.poisson(), delta = -1, deriv = 1)
>  # does not matter what zeroth is, theta is not identifiable for
>  # degenerate model
>  is.finite(foo$zeroth)
[1] TRUE
>  all.equal(0, foo$first)
[1] TRUE
> 
>  ##### multinomial #####
> 
>  set.seed(42)
> 
>  d <- 4
> 
>  theta <- matrix(rnorm(d * 25), ncol = d)
> 
>  cumfun <- function(theta, delta) {
+      stopifnot(is.numeric(theta))
+      stopifnot(is.finite(theta))
+      stopifnot(length(theta) == d)
+      stopifnot(is.numeric(delta))
+      stopifnot(is.finite(delta))
+      stopifnot(length(delta) == d)
+      stopifnot(delta <= 0)
+      stopifnot(any(delta == 0))
+      inies <- delta == 0
+      d.too <- sum(inies)
+      if (d.too == 1) return(theta[inies])
+      theta.too <- theta[inies]
+      return(cumulant(theta.too, fam.multinomial(d.too))$zeroth)
+  }
> 
>  moofun <- function(theta, delta) {
+      stopifnot(is.numeric(theta))
+      stopifnot(is.finite(theta))
+      stopifnot(length(theta) == d)
+      stopifnot(is.numeric(delta))
+      stopifnot(is.finite(delta))
+      stopifnot(length(delta) == d)
+      stopifnot(delta <= 0)
+      stopifnot(any(delta == 0))
+      inies <- delta == 0
+      d.too <- sum(inies)
+      if (d.too == 1) as.numeric(inies)
+      theta.too <- theta[inies]
+      foo <- cumulant(theta.too, fam.multinomial(d.too), deriv = 1)$first
+      bar <- rep(0, d)
+      bar[inies] <- foo
+      return(bar)
+  }
> 
>  voofun <- function(theta, delta) {
+      stopifnot(is.numeric(theta))
+      stopifnot(is.finite(theta))
+      stopifnot(length(theta) == d)
+      stopifnot(is.numeric(delta))
+      stopifnot(is.finite(delta))
+      stopifnot(length(delta) == d)
+      stopifnot(delta <= 0)
+      stopifnot(any(delta == 0))
+      inies <- delta == 0
+      d.too <- sum(inies)
+      if (d.too == 1) return(matrix(0, d, d))
+      theta.too <- theta[inies]
+      foo <- cumulant(theta.too, fam.multinomial(d.too), deriv = 2)$second
+      bar <- matrix(0, d, d)
+      baz <- matrix(0, d, d.too)
+      baz[inies, ] <- foo
+      bar[ , inies] <- baz
+      return(bar)
+  }
> 
>  thoofun <- function(theta, delta) {
+      stopifnot(is.numeric(theta))
+      stopifnot(is.finite(theta))
+      stopifnot(length(theta) == d)
+      stopifnot(is.numeric(delta))
+      stopifnot(is.finite(delta))
+      stopifnot(length(delta) == d)
+      stopifnot(delta <= 0)
+      stopifnot(any(delta == 0))
+      inies <- delta == 0
+      d.too <- sum(inies)
+      if (d.too == 1) return(array(0, rep(d, 3)))
+      theta.too <- theta[inies]
+      foo <- cumulant(theta.too, fam.multinomial(d.too), deriv = 3)$third
+      bar <- array(0, rep(d, 3))
+      baz <- array(0, c(d, d, d.too))
+      qux <- array(0, c(d, d.too, d.too))
+      qux[inies, , ] <- foo
+      baz[ , inies, ] <- qux
+      bar[ , , inies] <- baz
+      return(bar)
+  }
> 
>  ##### one extra degeneracy #####
> 
>  delta <- c(0, 0, 0, -1)
> 
>  zeroth <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta)$zeroth)
>  my.zeroth <- apply(theta, 1, cumfun, delta = delta)
>  all.equal(zeroth, my.zeroth)
[1] TRUE
> 
>  first <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  my.first <- apply(theta, 1, moofun, delta = delta)
>  all.equal(first, my.first)
[1] TRUE
> 
>  second <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 2)$second)
>  my.second <- apply(theta, 1, voofun, delta = delta)
>  all.equal(second, my.second)
[1] TRUE
> 
>  third <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 3)$third)
>  my.third <- apply(theta, 1, thoofun, delta = delta)
>  all.equal(third, my.third)
[1] TRUE
> 
>  link(first[ , 1], fam.multinomial(d), delta = delta, deriv = 1)
$zeroth
[1]  0.000000 -1.801428 -1.049033  0.000000

$first
          [,1]     [,2]     [,3] [,4]
[1,]  0.000000 0.000000 0.000000    0
[2,] -1.515339 9.180365 0.000000    0
[3,] -1.515339 0.000000 4.326126    0
[4,]  0.000000 0.000000 0.000000    0

> 
>  ##### two extra degeneracies #####
> 
>  delta <- c(0, 0, -1, -2)
> 
>  zeroth <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta)$zeroth)
>  my.zeroth <- apply(theta, 1, cumfun, delta = delta)
>  all.equal(zeroth, my.zeroth)
[1] TRUE
> 
>  first <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  my.first <- apply(theta, 1, moofun, delta = delta)
>  all.equal(first, my.first)
[1] TRUE
> 
>  second <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 2)$second)
>  my.second <- apply(theta, 1, voofun, delta = delta)
>  all.equal(second, my.second)
[1] TRUE
> 
>  third <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 3)$third)
>  my.third <- apply(theta, 1, thoofun, delta = delta)
>  all.equal(third, my.third)
[1] TRUE
> 
>  link(first[ , 1], fam.multinomial(d), delta = delta, deriv = 1)
$zeroth
[1]  0.000000 -1.801428  0.000000  0.000000

$first
          [,1]    [,2] [,3] [,4]
[1,]  0.000000 0.00000    0    0
[2,] -1.165063 7.05829    0    0
[3,]  0.000000 0.00000    0    0
[4,]  0.000000 0.00000    0    0

> 
>  ##### fully degenerate #####
> 
>  delta <- c(0, -0.5, -1, -2)
> 
>  zeroth <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta)$zeroth)
>  my.zeroth <- apply(theta, 1, cumfun, delta = delta)
>  all.equal(zeroth, my.zeroth)
[1] TRUE
> 
>  first <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  my.first <- apply(theta, 1, moofun, delta = delta)
>  all.equal(first, my.first)
[1] TRUE
> 
>  second <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 2)$second)
>  my.second <- apply(theta, 1, voofun, delta = delta)
>  all.equal(second, my.second)
[1] TRUE
> 
>  third <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 3)$third)
>  my.third <- apply(theta, 1, thoofun, delta = delta)
>  all.equal(third, my.third)
[1] TRUE
> 
>  ##### now link #####
> 
>  linkfun <- function(xi, delta) {
+      stopifnot(is.numeric(xi))
+      stopifnot(is.finite(xi))
+      stopifnot(length(xi) == d)
+      stopifnot(is.numeric(delta))
+      stopifnot(is.finite(delta))
+      stopifnot(length(delta) == d)
+      stopifnot(delta <= 0)
+      stopifnot(any(delta == 0))
+      inies <- delta == 0
+      d.too <- sum(inies)
+      if (d.too == 1) return(rep(0, d))
+      xi.too <- xi[inies]
+      foo <- link(xi.too, fam.multinomial(d.too))$zeroth
+      bar <- rep(0, d)
+      bar[inies] <- foo
+      return(bar)
+  }
> 
>  dlinkfun <- function(xi, delta) {
+      stopifnot(is.numeric(xi))
+      stopifnot(is.finite(xi))
+      stopifnot(length(xi) == d)
+      stopifnot(is.numeric(delta))
+      stopifnot(is.finite(delta))
+      stopifnot(length(delta) == d)
+      stopifnot(delta <= 0)
+      stopifnot(any(delta == 0))
+      inies <- delta == 0
+      d.too <- sum(inies)
+      if (d.too == 1) return(matrix(0, d, d))
+      xi.too <- xi[inies]
+      foo <- link(xi.too, fam.multinomial(d.too), deriv = 1)$first
+      bar <- matrix(0, d, d)
+      baz <- matrix(0, d, d.too)
+      baz[inies, ] <- foo
+      bar[ , inies] <- baz
+      return(bar)
+      bar[inies] <- foo
+      return(bar)
+  }
> 
>  ##### one extra degeneracy #####
> 
>  delta <- c(0, 0, 0, -1)
> 
>  xi <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  xi <- t(xi)
> 
>  zeroth <- apply(xi, 1, function(xi) link(xi, fam.multinomial(d),
+      delta = delta)$zeroth)
>  my.zeroth <- apply(xi, 1, linkfun, delta = delta)
>  all.equal(zeroth, my.zeroth)
[1] TRUE
> 
>  first <- apply(xi, 1, function(xi) link(xi, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  my.first <- apply(xi, 1, dlinkfun, delta = delta)
>  all.equal(first, my.first)
[1] TRUE
> 
>  ##### two extra degeneracies #####
> 
>  delta <- c(0, 0, -1, -2)
> 
>  xi <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  xi <- t(xi)
> 
>  zeroth <- apply(xi, 1, function(xi) link(xi, fam.multinomial(d),
+      delta = delta)$zeroth)
>  my.zeroth <- apply(xi, 1, linkfun, delta = delta)
>  all.equal(zeroth, my.zeroth)
[1] TRUE
> 
>  first <- apply(xi, 1, function(xi) link(xi, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  my.first <- apply(xi, 1, dlinkfun, delta = delta)
>  all.equal(first, my.first)
[1] TRUE
> 
>  ##### fully degenerate #####
> 
>  delta <- c(0, -0.5, -1, -2)
> 
>  xi <- apply(theta, 1, function(theta) cumulant(theta, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  xi <- t(xi)
> 
>  zeroth <- apply(xi, 1, function(xi) link(xi, fam.multinomial(d),
+      delta = delta)$zeroth)
>  my.zeroth <- apply(xi, 1, linkfun, delta = delta)
>  all.equal(zeroth, my.zeroth)
[1] TRUE
> 
>  first <- apply(xi, 1, function(xi) link(xi, fam.multinomial(d),
+      delta = delta, deriv = 1)$first)
>  my.first <- apply(xi, 1, dlinkfun, delta = delta)
>  all.equal(first, my.first)
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   1.81    0.14    1.93