R Under development (unstable) (2024-06-17 r86768 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(epigrowthfit)
> options(warn = 2L, error = if (interactive()) recover, egf.cores = 2L)
> 
> 
> ## exponential #########################################################
> 
> r <- log(2) / 20
> c0 <- 100
> s <- 0.2
> 
> mu <- log(c(r, c0))
> Sigma <- diag(rep.int(s^2, length(mu)))
> 
> zz <- simulate(egf_model(curve = "exponential", family = "pois"),
+                nsim = 20L,
+                seed = 775494L,
+                mu = mu,
+                Sigma = Sigma,
+                cstart = 10)
> mm <- egf(zz,
+           formula_priors = list(Sigma ~ LKJ(eta = 2)))
> 
> p1 <- as.list(coef(zz))
> p2 <- as.list(coef(mm))
> 
> stopifnot(exprs = {
+ 	max(abs(mm[["gradient"]])) < 5e-05
+ 	all.equal(p1[["beta"]], p2[["beta"]], tolerance = 5e-02)
+ 	all.equal(theta2cov(p1[["theta"]]), theta2cov(p2[["theta"]]), tolerance = 5e-02)
+ })
> 
> 
> ## subexponential ######################################################
> 
> alpha <- log(2) / 20
> c0 <- 100
> p <- 0.95
> s <- 0.2
> 
> mu <- c(log(alpha), log(c0), qlogis(p))
> Sigma <- diag(rep.int(s^2, length(mu)))
> 
> zz <- simulate(egf_model(curve = "subexponential", family = "pois"),
+                nsim = 20L,
+                seed = 653927L,
+                mu = mu,
+                Sigma = Sigma,
+                cstart = 10)
> mm <- egf(zz,
+           formula_priors = list(beta[3L] ~ Normal(mu = qlogis(p), sigma = 0.05),
+                                 theta[3L] ~ Normal(mu = log(s), sigma = 0.25),
+                                 Sigma ~ LKJ(eta = 2)))
> 
> p1 <- as.list(coef(zz))
> p2 <- as.list(coef(mm))
> 
> stopifnot(exprs = {
+ 	max(abs(mm[["gradient"]])) < 5e-04
+ 	all.equal(p1[["beta"]], p2[["beta"]], tolerance = 5e-02)
+ 	all.equal(theta2cov(p1[["theta"]]), theta2cov(p2[["theta"]]), tolerance = 2e-02)
+ })
> 
> 
> ## gompertz ############################################################
> 
> alpha <- log(2) / 20
> tinfl <- 100
> K <- 25000
> s <- 0.2
> 
> mu <- log(c(alpha, tinfl, K))
> Sigma <- diag(rep.int(s^2, length(mu)))
> 
> zz <- simulate(egf_model(curve = "gompertz", family = "pois"),
+                nsim = 20L,
+                seed = 685399L,
+                mu = mu,
+                Sigma = Sigma,
+                cstart = 10)
> oo <- options(warn = 1L) # FIXME: diagnose NA/NaN function evaluation
> mm <- egf(zz,
+           formula_priors = list(Sigma ~ LKJ(eta = 2)))
Warning in nlminb(start = par, objective = fn, gradient = gr, control = control,  :
  NA/NaN function evaluation
> options(oo)
> 
> p1 <- as.list(coef(zz))
> p2 <- as.list(coef(mm))
> 
> stopifnot(exprs = {
+ 	max(abs(mm[["gradient"]])) < 5e-04
+ 	all.equal(p1[["beta"]], p2[["beta"]], tolerance = 5e-02)
+ 	all.equal(theta2cov(p1[["theta"]]), theta2cov(p2[["theta"]]), tolerance = 2e-02)
+ })
> 
> 
> ## logistic ############################################################
> 
> r <- log(2) / 20
> tinfl <- 100
> K <- 25000
> s <- 0.2
> 
> mu <- log(c(r, tinfl, K))
> Sigma <- diag(rep.int(s^2, length(mu)))
> 
> zz <- simulate(egf_model(curve = "logistic", family = "pois"),
+                nsim = 20L,
+                seed = 397981L,
+                mu = mu,
+                Sigma = Sigma,
+                cstart = 10)
> mm <- egf(zz,
+           formula_priors = list(Sigma ~ LKJ(eta = 2)))
> 
> p1 <- as.list(coef(zz))
> p2 <- as.list(coef(mm))
> 
> stopifnot(exprs = {
+ 	max(abs(mm[["gradient"]])) < 1e-02
+ 	all.equal(p1[["beta"]], p2[["beta"]], tolerance = 1e-02)
+ 	all.equal(theta2cov(p1[["theta"]]), theta2cov(p2[["theta"]]), tolerance = 2e-02)
+ })
> 
> 
> ## richards ############################################################
> 
> r <- log(2) / 20
> tinfl <- 100
> K <- 25000
> a <- 1.005
> s <- 0.2
> 
> mu <- log(c(r, tinfl, K, a))
> Sigma <- diag(rep.int(s^2, length(mu)))
> 
> zz <- simulate(egf_model(curve = "richards", family = "pois"),
+                nsim = 20L,
+                seed = 949642L,
+                mu = mu,
+                Sigma = Sigma,
+                cstart = 10)
> mm <- egf(zz,
+           formula_priors = list(beta[4L] ~ Normal(mu = log(a), sigma = 0.005),
+                                 theta[4L] ~ Normal(mu = log(s), sigma = 0.25),
+                                 Sigma ~ LKJ(eta = 2)))
> 
> p1 <- as.list(coef(zz))
> p2 <- as.list(coef(mm))
> 
> stopifnot(exprs = {
+ 	max(abs(mm[["gradient"]])) < 2e-02
+ 	all.equal(p1[["beta"]], p2[["beta"]], tolerance = 5e-03)
+ 	all.equal(theta2cov(p1[["theta"]]), theta2cov(p2[["theta"]]), tolerance = 2e-02)
+ })
> 
> proc.time()
   user  system elapsed 
 106.26    0.68   61.46