https://github.com/cran/aster
Tip revision: 7016e6b97f24943bdab11323884baf030f38260b authored by Charles J. Geyer on 06 July 2018, 07:20:08 UTC
version 1.0-2
version 1.0-2
Tip revision: 7016e6b
pickle.Rout.save
R Under development (unstable) (2018-02-14 r74250) -- "Unsuffered Consequences"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
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(aster)
Loading required package: trust
>
> options(digits=4) # avoid rounding differences
>
> data(radish)
>
> pred <- c(0,1,2)
> fam <- c(1,3,2)
>
> ### need object of type aster to supply to penmlogl and pickle
>
> aout <- aster(resp ~ varb + fit : (Site * Region + Block + Pop),
+ pred, fam, varb, id, root, data = radish)
>
> ### model matrices for fixed and random effects
>
> modmat.fix <- model.matrix(resp ~ varb + fit : (Site * Region),
+ data = radish)
> modmat.blk <- model.matrix(resp ~ 0 + fit:Block, data = radish)
> modmat.pop <- model.matrix(resp ~ 0 + fit:Pop, data = radish)
>
> rownames(modmat.fix) <- NULL
> rownames(modmat.blk) <- NULL
> rownames(modmat.pop) <- NULL
>
> idrop <- match(aout$dropped, colnames(modmat.fix))
> idrop <- idrop[! is.na(idrop)]
> modmat.fix <- modmat.fix[ , - idrop]
>
> nfix <- ncol(modmat.fix)
> nblk <- ncol(modmat.blk)
> npop <- ncol(modmat.pop)
>
> ### try penmlogl
>
> sigma.start <- c(1, 1)
>
> alpha.start <- aout$coefficients[match(colnames(modmat.fix),
+ names(aout$coefficients))]
> parm.start <- c(alpha.start, rep(0, nblk + npop))
>
> tout <- trust(objfun = penmlogl, parm.start, rinit = 1, rmax = 10,
+ sigma = sigma.start, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout)
>
> eff <- tout$argument * tout$scale
> eff.blk <- eff[seq(nfix + 1, nfix + nblk)]
> eff.pop <- eff[seq(nfix + nblk + 1, nfix + nblk + npop)]
> sigma.crude <- sqrt(c(var(eff.blk), var(eff.pop)))
>
> pout <- pickle(sigma.crude, tout$argument, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout)
>
> # try optim with method = "Nelder-Mead" and pickle
>
> cache <- new.env(parent = emptyenv())
> oout <- optim(sigma.crude, pickle, parm = tout$argument,
+ fixed = modmat.fix, random = list(modmat.blk, modmat.pop),
+ obj = aout, cache = cache)
> oout$convergence == 0
[1] TRUE
>
> sigma.mle <- oout$par
> tout <- trust(objfun = penmlogl, cache$parm, rinit = 1, rmax = 10,
+ sigma = sigma.mle, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout)
> stopifnot(tout$converged)
> parm.mle <- tout$argument
>
> # try pickle1
>
> zwz <- makezwz(sigma.mle, parm.mle, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout)
>
> pout <- pickle1(sigma.mle, parm.mle, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout, zwz = zwz, deriv = 1)
>
> foo <- function(sigma) {
+ pout <- pickle1(sigma.mle, parm.mle, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout, zwz = zwz,
+ deriv = 1)
+ result <- pout$value
+ attr(result, "gradient") <- pout$gradient
+ return(result)
+ }
>
> nout <- nlm(foo, sigma.mle)
> nout$code
[1] 3
> nout$iterations
[1] 1
>
> nrand <- c(nblk, npop)
> qux <- function(parm, sigma, zwz) {
+ foo <- function(alphaceesigma) pickle3(alphaceesigma, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout, zwz = zwz, deriv = 2)
+ sigma <- as.vector(sigma)
+ parm <- as.vector(parm)
+ zwz <- as.matrix(zwz)
+ iter <- NULL
+ repeat {
+ tout <- trust(foo, c(parm, sigma), rinit = 1, rmax = 10)
+ stopifnot(tout$converged)
+ iter <- c(iter, tout$iterations)
+ sigma.old <- sigma
+ sigma <- tout$argument[nfix + sum(nrand) + seq(along = nrand)]
+ parm <- tout$argument[seq(1, nfix + sum(nrand))]
+ zwz <- makezwz(sigma, parm, fixed = modmat.fix,
+ random = list(modmat.blk, modmat.pop), obj = aout)
+ # cat("iteration", length(iter), ":",
+ # all.equal(sigma, sigma.old, check.attributes = FALSE), "\n")
+ if (isTRUE(all.equal(sigma, sigma.old))) break
+ }
+ return(list(sigma = sigma, parm = parm, zwz = zwz, iterations = iter))
+ }
>
> qout <- qux(parm.mle, sigma.mle, zwz)
>
> sigma.mle <- qout$sigma
> parm.mle <- qout$parm
> zwz.mle <- qout$zwz
> alpha.mle <- parm.mle[1:nfix]
> c.mle <- parm.mle[-(1:nfix)]
> a.mle <- rep(sigma.mle, times = nrand)
> b.mle <- a.mle * c.mle
>
> # use optim to get hessian of q(alpha, sigma)
> # note: nice analytic formula in inst/doc/re.pdf doesn't work
> # with inexact computer arithmetic due to catastrophic cancellation
>
> alphasigma.mle <- c(alpha.mle, sigma.mle)
>
> objfun <- function(alphasigma) pickle2(alphasigma, parm = c.mle,
+ fixed = modmat.fix, random = list(modmat.blk, modmat.pop),
+ obj = aout, zwz = zwz.mle)$value
> gradfun <- function(alphasigma) pickle2(alphasigma, parm = c.mle,
+ fixed = modmat.fix, random = list(modmat.blk, modmat.pop),
+ obj = aout, zwz = zwz.mle, deriv = 1)$gradient
>
> oout2 <- optim(alphasigma.mle, objfun, gradfun, method = "BFGS", hessian = TRUE)
>
> foo <- new.env(parent = emptyenv())
> bar <- suppressWarnings(try(load("pickle.rda", foo), silent = TRUE))
> if (inherits(bar, "try-error")) {
+ save(oout, qout, oout2, file = "pickle.rda")
+ } else {
+ print(all.equal(oout, foo$oout))
+ qout$iterations <- NULL
+ foo$qout$iterations <- NULL
+ ## changes needed for alternative BLASes
+ qout$counts <- foo$qout$counts <- NULL
+ oout2$counts <- foo$oout2$counts <- NULL
+ print(all.equal(qout, foo$qout, tolerance = 1e-4))
+ print(all.equal(oout2, foo$oout2, tolerance = 1e-4))
+ }
[1] TRUE
[1] TRUE
[1] TRUE
>
> ########## ?????????? LOOKS LIKE UNFINISHED ??????????
>
>
> proc.time()
user system elapsed
8.507 0.060 8.598