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)
> 
>  data(test1)
> 
>  fred <- asterdata(test1,
+      vars = c("m1", "m2", "m3", "n1", "n2", "b1", "p1", "z1"),
+      pred = c(0, 0, 0, 1, 1, 2, 3, 6), group = c(0, 1, 2, 0, 4, 0, 0, 0),
+      code = c(1, 1, 1, 2, 2, 3, 4, 5),
+      families = list(fam.multinomial(3), "normal.location.scale",
+      "bernoulli", "poisson", "zero.truncated.poisson"))
>  is.validasterdata(fred)
[1] TRUE
>  names(fred)
 [1] "redata"        "repred"        "regroup"       "recode"       
 [5] "families"      "redelta"       "initial"       "response.name"
 [9] "varb.name"     "pred"          "group"         "code"         
>  names(fred$redata)
[1] "varb" "resp" "id"  
> 
>  theta <- rnorm(length(fred))
>  try(validtheta(fred, theta))
Error in validtheta(fred, theta) : 
  theta[2] not negative for normal location-scale
>  try(validtheta(fred, theta, mod = "cond"))
Error in validtheta(fred, theta, mod = "cond") : 
  theta[2] not negative for normal location-scale
> 
>  varb <- as.character(fred$redata[[fred$varb.name]]) 
>  pred <- fred$repred
>  resp <- fred$redata[[fred$response.name]]
>  init <- fred$initial
>  ypred <- ifelse(pred == 0, init, c(NA, resp)[pred + 1])
> 
>  theta.cond <- theta
>  theta.unco <- theta
>  theta.unco[varb == "n2"] <- (- abs(theta.unco[varb == "n2"]))
>  theta.cond[varb == "n2" & ypred > 0] <-
+      (- abs(theta.unco[varb == "n2" & ypred > 0]))
> 
>  is.validtheta(fred, theta.unco)
[1] TRUE
>  is.validtheta(fred, theta.cond, mod = "cond")
[1] TRUE
>  
>  xi <- rnorm(length(fred))
>  try(validxi(fred, xi))
Error in validxi(fred, xi) : 
  zero-truncated Poisson xi not strictly greater than 1
>  try(validxi(fred, xi, mod = "cond"))
Error in validxi(fred, xi, mod = "cond") : 
  zero-truncated Poisson xi not strictly greater than 1
> 
>  fixup <- function(xi, varb, ypred) {
+      stopifnot(length(xi) == length(varb))
+      stopifnot(is.numeric(xi))
+      stopifnot(is.finite(xi))
+      stopifnot(is.character(varb))
+      if (missing(ypred)) {
+          ypred <- rep(1, length(xi))
+      } else {
+          stopifnot(is.numeric(ypred))
+          stopifnot(length(ypred) == length(xi))
+      }
+      todo <- ypred != 0
+      xi[todo & varb == "z1"] <- 1 + abs(xi[todo & varb == "z1"])
+      xi[todo & varb == "p1"] <- abs(xi[todo & varb == "p1"])
+      xi[todo & varb == "b1"] <- 1 / (1 + abs(xi[todo & varb == "b1"]))
+      xi[todo & varb == "n2"] <- 1 + xi[todo & varb == "n1"]^2
+      my.xi.m1 <- exp(xi[todo & varb == "m1"])
+      my.xi.m2 <- exp(xi[todo & varb == "m2"])
+      my.xi.m3 <- exp(xi[todo & varb == "m3"])
+      xi[todo & varb == "m1"] <- my.xi.m1 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
+      xi[todo & varb == "m2"] <- my.xi.m2 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
+      xi[todo & varb == "m3"] <- my.xi.m3 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
+      return(xi)
+  }
> 
>  xi.unco <- fixup(xi, varb)
>  xi.cond <- fixup(xi, varb, ypred)
> 
>  is.validxi(fred, xi.unco)
[1] TRUE
>  is.validxi(fred, xi.cond, mod = "cond")
[1] TRUE
> 
>  ########## now with direction of recession ##########
> 
>  delta <- fred$redelta
>  all(delta == 0)
[1] TRUE
>  delta[grepl("[mbp]", varb) & resp == 0 & runif(length(fred)) < 1 / 5] <- -1
>  fred$redelta <- delta
>  is.validasterdata(fred)
[1] TRUE
>  sum(delta != 0)
[1] 65
>  
>  try(validtheta(fred, theta))
Error in validtheta(fred, theta) : 
  theta[2] not negative for normal location-scale
>  try(validtheta(fred, theta, mod = "cond"))
Error in validtheta(fred, theta, mod = "cond") : 
  theta[2] not negative for normal location-scale
> 
>  is.validtheta(fred, theta.unco)
[1] TRUE
>  is.validtheta(fred, theta.cond, mod = "cond")
[1] TRUE
>  
>  my.sally <- 2 * as.numeric(delta < 0)
>  # when loop done
>  # my.sally == 2 means delta < 0 and predecessor NOT a. s. zero
>  # my.sally == 1 means predecessor a. s. zero (regardless of delta)
>  repeat {
+      save.my.sally <- my.sally
+      foom <- c(0, my.sally)[pred + 1]
+      my.sally[foom != 0] <- 1
+      if (identical(my.sally, save.my.sally)) break
+  }
>  my.sally <- my.sally == 1
> 
>  try(validxi(fred, xi))
Error in validxi(fred, xi) : 
  zero-truncated Poisson xi not strictly greater than 1
>  try(validxi(fred, xi, mod = "cond"))
Error in validxi(fred, xi, mod = "cond") : 
  zero-truncated Poisson xi not strictly greater than 1
> 
>  fixup <- function(xi, varb, ypred, delta) {
+      stopifnot(is.character(varb))
+      stopifnot(length(xi) == length(varb))
+      stopifnot(is.numeric(xi))
+      stopifnot(is.finite(xi))
+      stopifnot(is.numeric(ypred))
+      stopifnot(is.finite(ypred))
+      stopifnot(length(ypred) == length(varb))
+      stopifnot(is.numeric(delta))
+      stopifnot(is.finite(delta))
+      stopifnot(length(delta) == length(varb))
+      todo <- ypred != 0
+      lower <- delta < 0
+      xi[todo & varb == "z1"] <- 1 + abs(xi[todo & varb == "z1"])
+      xi[todo & varb == "z1" & lower] <- 1
+      xi[todo & varb == "p1"] <- abs(xi[todo & varb == "p1"])
+      xi[todo & varb == "p1" & lower] <- 0
+      xi[todo & varb == "b1"] <- 1 / (1 + abs(xi[todo & varb == "b1"]))
+      xi[todo & varb == "b1" & lower] <- 0
+      xi[todo & varb == "n2"] <- 1 + xi[todo & varb == "n1"]^2
+      my.xi <- ifelse(lower, 0, exp(xi))
+      my.xi.m1 <- my.xi[todo & varb == "m1"]
+      my.xi.m2 <- my.xi[todo & varb == "m2"]
+      my.xi.m3 <- my.xi[todo & varb == "m3"]
+      xi[todo & varb == "m1"] <- my.xi.m1 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
+      xi[todo & varb == "m2"] <- my.xi.m2 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
+      xi[todo & varb == "m3"] <- my.xi.m3 / (my.xi.m1 + my.xi.m2 + my.xi.m3)
+      return(xi)
+  }
> 
>  xi.unco <- fixup(xi, varb, as.numeric(! my.sally), delta)
>  xi.cond <- fixup(xi, varb, ypred, delta)
> 
>  is.validxi(fred, xi.unco)
[1] TRUE
>  is.validxi(fred, xi.cond, mod = "cond")
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   0.92    0.06    0.96