..Rcheck/aghq-Ex.R
In awstringer1/aghq: Adaptive Gauss Hermite Quadrature for Bayesian Inference
pkgname <- "aghq"
source(file.path(R.home("share"), "R", "examples-header.R"))
options(warn = 1)
base::assign(".ExTimings", "aghq-Ex.timings", pos = 'CheckExEnv')
base::cat("name\tuser\tsystem\telapsed\n", file=base::get(".ExTimings", pos = 'CheckExEnv'))
base::assign(".format_ptime",
function(x) {
if(!is.na(x[4L])) x[1L] <- x[1L] + x[4L]
if(!is.na(x[5L])) x[2L] <- x[2L] + x[5L]
options(OutDec = '.')
format(x[1L:3L], digits = 7L)
},
pos = 'CheckExEnv')
### * </HEADER>
library('aghq')
base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
cleanEx()
nameEx("aghq")
### * aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: aghq
### Title: Adaptive Gauss-Hermite Quadrature
### Aliases: aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("compute_moment")
### * compute_moment
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: compute_moment
### Title: Compute moments
### Aliases: compute_moment compute_moment.default compute_moment.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
norm_sparse_2d_7 <- normalize_logpost(opt_sparsetrust_2d,7,1)
# ff = function(x) 1 should return 1,
# the normalizing constant of the (already normalized) posterior:
compute_moment(norm_sparse_2d_7)
# Compute the mean of theta1 and theta2
compute_moment(norm_sparse_2d_7,ff = function(x) x)
# Compute the mean of lambda1 = exp(theta1) and lambda2 = exp(theta2)
lambdameans <- compute_moment(norm_sparse_2d_7,ff = function(x) exp(x))
lambdameans
# Compare them to the truth:
(sum(y1) + 1)/(length(y1) + 1)
(sum(y2) + 1)/(length(y2) + 1)
# Compute the standard deviation of lambda1
lambda1sd <- sqrt(compute_moment(norm_sparse_2d_7,ff = function(x) (exp(x) - lambdameans[1])^2))[1]
# ...and so on.
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("compute_moment", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("compute_pdf_and_cdf")
### * compute_pdf_and_cdf
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: compute_pdf_and_cdf
### Title: Density and Cumulative Distribution Function
### Aliases: compute_pdf_and_cdf compute_pdf_and_cdf.default
### compute_pdf_and_cdf.list compute_pdf_and_cdf.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
margpost <- marginal_posterior(opt_sparsetrust_2d,3,1) # margpost for theta1
thepdfandcdf <- compute_pdf_and_cdf(margpost)
with(thepdfandcdf,{
plot(pdf~theta,type='l')
plot(cdf~theta,type='l')
})
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("compute_pdf_and_cdf", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("compute_quantiles")
### * compute_quantiles
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: compute_quantiles
### Title: Quantiles
### Aliases: compute_quantiles compute_quantiles.default
### compute_quantiles.list compute_quantiles.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
margpost <- marginal_posterior(opt_sparsetrust_2d,3,1) # margpost for theta1
etaquant <- compute_quantiles(margpost)
etaquant
# lambda = exp(eta)
exp(etaquant)
# Compare to truth
qgamma(.025,1+sum(y1),1+n1)
qgamma(.975,1+sum(y1),1+n1)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("compute_quantiles", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("default_control")
### * default_control
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: default_control
### Title: Default control arguments for 'aghq::aghq()'.
### Aliases: default_control
### ** Examples
default_control()
default_control(method = "trust")
default_control(negate = TRUE)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("default_control", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("default_control_marglaplace")
### * default_control_marglaplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: default_control_marglaplace
### Title: Default control arguments for 'aghq::marginal_laplace()'.
### Aliases: default_control_marglaplace
### ** Examples
default_control_marglaplace()
default_control_marglaplace(method = "trust")
default_control_marglaplace(method = "trust",inner_method = "trust")
default_control_marglaplace(negate = TRUE)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("default_control_marglaplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("default_control_tmb")
### * default_control_tmb
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: default_control_tmb
### Title: Default control arguments for 'aghq::marginal_laplace_tmb()'.
### Aliases: default_control_tmb
### ** Examples
default_control_marglaplace()
default_control_marglaplace(method = "trust")
default_control_marglaplace(method = "trust",inner_method = "trust")
default_control_marglaplace(negate = TRUE)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("default_control_tmb", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("laplace_approximation")
### * laplace_approximation
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: laplace_approximation
### Title: Laplace Approximation
### Aliases: laplace_approximation
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("laplace_approximation", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("marginal_laplace")
### * marginal_laplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: marginal_laplace
### Title: Marginal Laplace approximation
### Aliases: marginal_laplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("marginal_laplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("marginal_posterior")
### * marginal_posterior
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: marginal_posterior
### Title: Marginal Posteriors
### Aliases: marginal_posterior
### ** Examples
## A 2d example ##
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
# Now actually do the marginal posteriors
marginal_posterior(opt_sparsetrust_2d,3,1)
marginal_posterior(opt_sparsetrust_2d,3,2)
marginal_posterior(opt_sparsetrust_2d,7,2)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("marginal_posterior", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("normalize_logpost")
### * normalize_logpost
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: normalize_logpost
### Title: Normalize the joint posterior using AGHQ
### Aliases: normalize_logpost
### ** Examples
# Same setup as optimize_theta
logfteta <- function(eta,y) {
sum(y) * eta - (length(y) + 1) * exp(eta) - sum(lgamma(y+1)) + eta
}
set.seed(84343124)
y <- rpois(10,5) # Mode should be sum(y) / (10 + 1)
truemode <- log((sum(y) + 1)/(length(y) + 1))
objfunc <- function(x) logfteta(x,y)
funlist <- list(
fn = objfunc,
gr = function(x) numDeriv::grad(objfunc,x),
he = function(x) numDeriv::hessian(objfunc,x)
)
opt_sparsetrust <- optimize_theta(funlist,1.5)
opt_trust <- optimize_theta(funlist,1.5,control = default_control(method = "trust"))
opt_bfgs <- optimize_theta(funlist,1.5,control = default_control(method = "BFGS"))
# Quadrature with 3, 5, and 7 points using sparse trust region optimization:
norm_sparse_3 <- normalize_logpost(opt_sparsetrust,3,1)
norm_sparse_5 <- normalize_logpost(opt_sparsetrust,5,1)
norm_sparse_7 <- normalize_logpost(opt_sparsetrust,7,1)
# Quadrature with 3, 5, and 7 points using dense trust region optimization:
norm_trust_3 <- normalize_logpost(opt_trust,3,1)
norm_trust_5 <- normalize_logpost(opt_trust,5,1)
norm_trust_7 <- normalize_logpost(opt_trust,7,1)
# Quadrature with 3, 5, and 7 points using BFGS optimization:
norm_bfgs_3 <- normalize_logpost(opt_bfgs,3,1)
norm_bfgs_5 <- normalize_logpost(opt_bfgs,5,1)
norm_bfgs_7 <- normalize_logpost(opt_bfgs,7,1)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("normalize_logpost", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("optimize_theta")
### * optimize_theta
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: optimize_theta
### Title: Obtain function information necessary for performing quadrature
### Aliases: optimize_theta
### ** Examples
# Poisson/Exponential example
logfteta <- function(eta,y) {
sum(y) * eta - (length(y) + 1) * exp(eta) - sum(lgamma(y+1)) + eta
}
y <- rpois(10,5) # Mode should be (sum(y) + 1) / (length(y) + 1)
objfunc <- function(x) logfteta(x,y)
funlist <- list(
fn = objfunc,
gr = function(x) numDeriv::grad(objfunc,x),
he = function(x) numDeriv::hessian(objfunc,x)
)
optimize_theta(funlist,1.5)
optimize_theta(funlist,1.5,control = default_control(method = "trust"))
optimize_theta(funlist,1.5,control = default_control(method = "BFGS"))
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("optimize_theta", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("plot.aghq")
### * plot.aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: plot.aghq
### Title: Plot method for AGHQ objects
### Aliases: plot.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
plot(thequadrature)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("plot.aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.aghq")
### * print.aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.aghq
### Title: Print method for AGHQ objects
### Aliases: print.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
thequadrature
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.aghqsummary")
### * print.aghqsummary
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.aghqsummary
### Title: Print method for AGHQ summary objects
### Aliases: print.aghqsummary
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thequadrature)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.aghqsummary", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.laplace")
### * print.laplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.laplace
### Title: Print method for AGHQ objects
### Aliases: print.laplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
thequadrature
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.laplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.laplacesummary")
### * print.laplacesummary
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.laplacesummary
### Title: Print method for laplacesummary objects
### Aliases: print.laplacesummary
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thelaplace <- laplace_approximation(funlist2d,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thelaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.laplacesummary", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.marginallaplacesummary")
### * print.marginallaplacesummary
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.marginallaplacesummary
### Title: Summary statistics for models using marginal Laplace
### approximations
### Aliases: print.marginallaplacesummary
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
themarginallaplace <- aghq::marginal_laplace(funlist2dmarg,3,list(W = 0,theta = 0))
summary(themarginallaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.marginallaplacesummary", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("sample_marginal")
### * sample_marginal
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: sample_marginal
### Title: Exact independent samples from an approximate posterior
### distribution
### Aliases: sample_marginal sample_marginal.aghq
### sample_marginal.marginallaplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("sample_marginal", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("summary.aghq")
### * summary.aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: summary.aghq
### Title: Summary statistics computed using AGHQ
### Aliases: summary.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thequadrature)
# or, compute the summary and save for further processing:
ss <- summary(thequadrature)
str(ss)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("summary.aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("summary.laplace")
### * summary.laplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: summary.laplace
### Title: Summary method for Laplace Approximation objects
### Aliases: summary.laplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thelaplace <- laplace_approximation(funlist2d,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thelaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("summary.laplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("summary.marginallaplace")
### * summary.marginallaplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: summary.marginallaplace
### Title: Summary statistics for models using marginal Laplace
### approximations
### Aliases: summary.marginallaplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
themarginallaplace <- aghq::marginal_laplace(funlist2dmarg,3,list(W = 0,theta = 0))
summary(themarginallaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("summary.marginallaplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("validate_control")
### * validate_control
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: validate_control
### Title: Validate a control list
### Aliases: validate_control
### ** Examples
validate_control(default_control())
validate_control(default_control_marglaplace(),type = "marglaplace")
validate_control(default_control_tmb(),type = "tmb")
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("validate_control", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
### * <FOOTER>
###
cleanEx()
options(digits = 7L)
base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
grDevices::dev.off()
###
### Local variables: ***
### mode: outline-minor ***
### outline-regexp: "\\(> \\)?### [*]+" ***
### End: ***
quit('no')
awstringer1/aghq documentation built on March 4, 2024, 4:48 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
pkgname <- "aghq"
source(file.path(R.home("share"), "R", "examples-header.R"))
options(warn = 1)
base::assign(".ExTimings", "aghq-Ex.timings", pos = 'CheckExEnv')
base::cat("name\tuser\tsystem\telapsed\n", file=base::get(".ExTimings", pos = 'CheckExEnv'))
base::assign(".format_ptime",
function(x) {
if(!is.na(x[4L])) x[1L] <- x[1L] + x[4L]
if(!is.na(x[5L])) x[2L] <- x[2L] + x[5L]
options(OutDec = '.')
format(x[1L:3L], digits = 7L)
},
pos = 'CheckExEnv')
### * </HEADER>
library('aghq')
base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
cleanEx()
nameEx("aghq")
### * aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: aghq
### Title: Adaptive Gauss-Hermite Quadrature
### Aliases: aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("compute_moment")
### * compute_moment
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: compute_moment
### Title: Compute moments
### Aliases: compute_moment compute_moment.default compute_moment.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
norm_sparse_2d_7 <- normalize_logpost(opt_sparsetrust_2d,7,1)
# ff = function(x) 1 should return 1,
# the normalizing constant of the (already normalized) posterior:
compute_moment(norm_sparse_2d_7)
# Compute the mean of theta1 and theta2
compute_moment(norm_sparse_2d_7,ff = function(x) x)
# Compute the mean of lambda1 = exp(theta1) and lambda2 = exp(theta2)
lambdameans <- compute_moment(norm_sparse_2d_7,ff = function(x) exp(x))
lambdameans
# Compare them to the truth:
(sum(y1) + 1)/(length(y1) + 1)
(sum(y2) + 1)/(length(y2) + 1)
# Compute the standard deviation of lambda1
lambda1sd <- sqrt(compute_moment(norm_sparse_2d_7,ff = function(x) (exp(x) - lambdameans[1])^2))[1]
# ...and so on.
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("compute_moment", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("compute_pdf_and_cdf")
### * compute_pdf_and_cdf
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: compute_pdf_and_cdf
### Title: Density and Cumulative Distribution Function
### Aliases: compute_pdf_and_cdf compute_pdf_and_cdf.default
### compute_pdf_and_cdf.list compute_pdf_and_cdf.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
margpost <- marginal_posterior(opt_sparsetrust_2d,3,1) # margpost for theta1
thepdfandcdf <- compute_pdf_and_cdf(margpost)
with(thepdfandcdf,{
plot(pdf~theta,type='l')
plot(cdf~theta,type='l')
})
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("compute_pdf_and_cdf", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("compute_quantiles")
### * compute_quantiles
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: compute_quantiles
### Title: Quantiles
### Aliases: compute_quantiles compute_quantiles.default
### compute_quantiles.list compute_quantiles.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
margpost <- marginal_posterior(opt_sparsetrust_2d,3,1) # margpost for theta1
etaquant <- compute_quantiles(margpost)
etaquant
# lambda = exp(eta)
exp(etaquant)
# Compare to truth
qgamma(.025,1+sum(y1),1+n1)
qgamma(.975,1+sum(y1),1+n1)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("compute_quantiles", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("default_control")
### * default_control
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: default_control
### Title: Default control arguments for 'aghq::aghq()'.
### Aliases: default_control
### ** Examples
default_control()
default_control(method = "trust")
default_control(negate = TRUE)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("default_control", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("default_control_marglaplace")
### * default_control_marglaplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: default_control_marglaplace
### Title: Default control arguments for 'aghq::marginal_laplace()'.
### Aliases: default_control_marglaplace
### ** Examples
default_control_marglaplace()
default_control_marglaplace(method = "trust")
default_control_marglaplace(method = "trust",inner_method = "trust")
default_control_marglaplace(negate = TRUE)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("default_control_marglaplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("default_control_tmb")
### * default_control_tmb
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: default_control_tmb
### Title: Default control arguments for 'aghq::marginal_laplace_tmb()'.
### Aliases: default_control_tmb
### ** Examples
default_control_marglaplace()
default_control_marglaplace(method = "trust")
default_control_marglaplace(method = "trust",inner_method = "trust")
default_control_marglaplace(negate = TRUE)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("default_control_tmb", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("laplace_approximation")
### * laplace_approximation
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: laplace_approximation
### Title: Laplace Approximation
### Aliases: laplace_approximation
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("laplace_approximation", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("marginal_laplace")
### * marginal_laplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: marginal_laplace
### Title: Marginal Laplace approximation
### Aliases: marginal_laplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("marginal_laplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("marginal_posterior")
### * marginal_posterior
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: marginal_posterior
### Title: Marginal Posteriors
### Aliases: marginal_posterior
### ** Examples
## A 2d example ##
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
# Now actually do the marginal posteriors
marginal_posterior(opt_sparsetrust_2d,3,1)
marginal_posterior(opt_sparsetrust_2d,3,2)
marginal_posterior(opt_sparsetrust_2d,7,2)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("marginal_posterior", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("normalize_logpost")
### * normalize_logpost
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: normalize_logpost
### Title: Normalize the joint posterior using AGHQ
### Aliases: normalize_logpost
### ** Examples
# Same setup as optimize_theta
logfteta <- function(eta,y) {
sum(y) * eta - (length(y) + 1) * exp(eta) - sum(lgamma(y+1)) + eta
}
set.seed(84343124)
y <- rpois(10,5) # Mode should be sum(y) / (10 + 1)
truemode <- log((sum(y) + 1)/(length(y) + 1))
objfunc <- function(x) logfteta(x,y)
funlist <- list(
fn = objfunc,
gr = function(x) numDeriv::grad(objfunc,x),
he = function(x) numDeriv::hessian(objfunc,x)
)
opt_sparsetrust <- optimize_theta(funlist,1.5)
opt_trust <- optimize_theta(funlist,1.5,control = default_control(method = "trust"))
opt_bfgs <- optimize_theta(funlist,1.5,control = default_control(method = "BFGS"))
# Quadrature with 3, 5, and 7 points using sparse trust region optimization:
norm_sparse_3 <- normalize_logpost(opt_sparsetrust,3,1)
norm_sparse_5 <- normalize_logpost(opt_sparsetrust,5,1)
norm_sparse_7 <- normalize_logpost(opt_sparsetrust,7,1)
# Quadrature with 3, 5, and 7 points using dense trust region optimization:
norm_trust_3 <- normalize_logpost(opt_trust,3,1)
norm_trust_5 <- normalize_logpost(opt_trust,5,1)
norm_trust_7 <- normalize_logpost(opt_trust,7,1)
# Quadrature with 3, 5, and 7 points using BFGS optimization:
norm_bfgs_3 <- normalize_logpost(opt_bfgs,3,1)
norm_bfgs_5 <- normalize_logpost(opt_bfgs,5,1)
norm_bfgs_7 <- normalize_logpost(opt_bfgs,7,1)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("normalize_logpost", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("optimize_theta")
### * optimize_theta
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: optimize_theta
### Title: Obtain function information necessary for performing quadrature
### Aliases: optimize_theta
### ** Examples
# Poisson/Exponential example
logfteta <- function(eta,y) {
sum(y) * eta - (length(y) + 1) * exp(eta) - sum(lgamma(y+1)) + eta
}
y <- rpois(10,5) # Mode should be (sum(y) + 1) / (length(y) + 1)
objfunc <- function(x) logfteta(x,y)
funlist <- list(
fn = objfunc,
gr = function(x) numDeriv::grad(objfunc,x),
he = function(x) numDeriv::hessian(objfunc,x)
)
optimize_theta(funlist,1.5)
optimize_theta(funlist,1.5,control = default_control(method = "trust"))
optimize_theta(funlist,1.5,control = default_control(method = "BFGS"))
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("optimize_theta", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("plot.aghq")
### * plot.aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: plot.aghq
### Title: Plot method for AGHQ objects
### Aliases: plot.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
plot(thequadrature)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("plot.aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.aghq")
### * print.aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.aghq
### Title: Print method for AGHQ objects
### Aliases: print.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
thequadrature
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.aghqsummary")
### * print.aghqsummary
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.aghqsummary
### Title: Print method for AGHQ summary objects
### Aliases: print.aghqsummary
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thequadrature)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.aghqsummary", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.laplace")
### * print.laplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.laplace
### Title: Print method for AGHQ objects
### Aliases: print.laplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
thequadrature
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.laplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.laplacesummary")
### * print.laplacesummary
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.laplacesummary
### Title: Print method for laplacesummary objects
### Aliases: print.laplacesummary
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thelaplace <- laplace_approximation(funlist2d,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thelaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.laplacesummary", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("print.marginallaplacesummary")
### * print.marginallaplacesummary
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: print.marginallaplacesummary
### Title: Summary statistics for models using marginal Laplace
### approximations
### Aliases: print.marginallaplacesummary
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
themarginallaplace <- aghq::marginal_laplace(funlist2dmarg,3,list(W = 0,theta = 0))
summary(themarginallaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("print.marginallaplacesummary", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("sample_marginal")
### * sample_marginal
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: sample_marginal
### Title: Exact independent samples from an approximate posterior
### distribution
### Aliases: sample_marginal sample_marginal.aghq
### sample_marginal.marginallaplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("sample_marginal", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("summary.aghq")
### * summary.aghq
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: summary.aghq
### Title: Summary statistics computed using AGHQ
### Aliases: summary.aghq
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thequadrature)
# or, compute the summary and save for further processing:
ss <- summary(thequadrature)
str(ss)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("summary.aghq", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("summary.laplace")
### * summary.laplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: summary.laplace
### Title: Summary method for Laplace Approximation objects
### Aliases: summary.laplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
fn = objfunc2d,
gr = function(x) numDeriv::grad(objfunc2d,x),
he = function(x) numDeriv::hessian(objfunc2d,x)
)
thelaplace <- laplace_approximation(funlist2d,c(0,0))
# Summarize and automatically call its print() method when called interactively:
summary(thelaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("summary.laplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("summary.marginallaplace")
### * summary.marginallaplace
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: summary.marginallaplace
### Title: Summary statistics for models using marginal Laplace
### approximations
### Aliases: summary.marginallaplace
### ** Examples
logfteta2d <- function(eta,y) {
# eta is now (eta1,eta2)
# y is now (y1,y2)
n <- length(y)
n1 <- ceiling(n/2)
n2 <- floor(n/2)
y1 <- y[1:n1]
y2 <- y[(n1+1):(n1+n2)]
eta1 <- eta[1]
eta2 <- eta[2]
sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
fn <- function(W) objfunc2dmarg(W,theta)
numDeriv::hessian(fn,W)
}
funlist2dmarg <- list(
fn = objfunc2dmarg,
gr = objfunc2dmarggr,
he = objfunc2dmarghe
)
themarginallaplace <- aghq::marginal_laplace(funlist2dmarg,3,list(W = 0,theta = 0))
summary(themarginallaplace)
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("summary.marginallaplace", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
cleanEx()
nameEx("validate_control")
### * validate_control
flush(stderr()); flush(stdout())
base::assign(".ptime", proc.time(), pos = "CheckExEnv")
### Name: validate_control
### Title: Validate a control list
### Aliases: validate_control
### ** Examples
validate_control(default_control())
validate_control(default_control_marglaplace(),type = "marglaplace")
validate_control(default_control_tmb(),type = "tmb")
base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
base::cat("validate_control", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
### * <FOOTER>
###
cleanEx()
options(digits = 7L)
base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
grDevices::dev.off()
###
### Local variables: ***
### mode: outline-minor ***
### outline-regexp: "\\(> \\)?### [*]+" ***
### End: ***
quit('no')
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.