Permalink
Browse files

Exact match of random effect names and related changes

    - Exact match of random effect names new default
    - Old behaviour possible by new argument regexp=TRUE
    - Examples modified to new default
    - Remove mle2 dependency
    - Optimize deprecated in randomregression example
  • Loading branch information...
1 parent fed240b commit bae04f2faf329516c29e6d9376519d1aaa5ff69a @kaskr committed Feb 25, 2014
View
@@ -39,7 +39,8 @@ updateCholesky <- function(L,H,t=0){
##' More advanced parameter mapping, such as collecting parameters between different vectors etc., must be implemented from the template.
##'
##' Random effects are specified via the argument \code{random}: A component of the parameter list is marked as random if its name is matched
-##' by any of the regular expressions of the character vector \code{random}. If some parameters are specified as random effects, these will
+##' by any of the characters of the vector \code{random} (Regular expression match is performed if \code{regexp=TRUE}).
+##' If some parameters are specified as random effects, these will
##' be integrated out of the objective function via the Laplace approximation. In this situation all of the functions \code{fn} and \code{gr}
##' automatically performs an optimization of random effects for each function evaluation. This is referred to as
##' the 'inner optimization'. Strategies for choosing initial values of the inner optimization can be controlled
@@ -82,7 +83,7 @@ updateCholesky <- function(L,H,t=0){
##' @param parameters List of all parameter objects required by the user template (both random and fixed effects).
##' @param map List defining how to optionally collect and fix parameters - see details.
##' @param type Character vector defining which operation stacks are generated from the users template - see details.
-##' @param random Character vector of regular expression defining the random effect parameters.
+##' @param random Character vector defining the random effect parameters. See also \code{regexp}.
##' @param random.start Expression defining the strategy for choosing random effect initial values as function of previous function evaluations - see details.
##' @param hessian Calculate Hessian at optimum?
##' @param method Outer optimization method.
@@ -94,6 +95,7 @@ updateCholesky <- function(L,H,t=0){
##' @param LaplaceNonZeroGradient Allow taylor expansion around non-stationary point?
##' @param DLL Name of shared object file compiled by user.
##' @param checkParameterOrder Optional check for correct parameter order.
+##' @param regexp Match random effects by regular expressions?
##' @param ...
##' @return List with components (fn,gr, etc) suitable for an optim call.
MakeADFun <- function(data,parameters,map=list(),
@@ -109,6 +111,7 @@ MakeADFun <- function(data,parameters,map=list(),
LaplaceNonZeroGradient=FALSE, ## Experimental feature: Allow expansion around non-stationary point
DLL=getUserDLL(),
checkParameterOrder=TRUE, ## Optional check
+ regexp=FALSE,
...){
env <- environment() ## This environment
if(!is.list(data))
@@ -259,10 +262,27 @@ MakeADFun <- function(data,parameters,map=list(),
.Call("EvalDoubleFunObject",Fun$ptr,unlist(parameters),control=list(order=as.integer(0)),PACKAGE=DLL)
}
if(is.character(random)){
- random <<- grepRandomParameters(parameters,random)
- if(length(random)==0){
- cat("Selected random effects did not match any model parameters.\n")
- random <<- NULL
+ if(!regexp){ ## Default: do exact match
+ if(!all(random %in% names(parameters))){
+ cat("Some 'random' effect names does not match 'parameter' list:\n")
+ print(setdiff(random,names(parameters)))
+ cat("(Note that regular expression match is disabled by default)\n")
+ stop()
+ }
+ if(any(duplicated(random))){
+ cat("Duplicates in 'random' - will be removed\n")
+ random <- unique(random)
+ }
+ tmp <- lapply(parameters,function(x)x*0)
+ tmp[random] <- lapply(tmp[random],function(x)x*0+1)
+ random <<- which(as.logical(unlist(tmp)))
+ }
+ if(regexp){ ## Original regular expression match
+ random <<- grepRandomParameters(parameters,random)
+ if(length(random)==0){
+ cat("Selected random effects did not match any model parameters.\n")
+ random <<- NULL
+ }
}
par <<- unlist(parameters)
}
@@ -24,16 +24,13 @@ parameters <- list(
)
obj <- MakeADFun(data=data,
parameters=parameters,
- random=c("^a","^b"),
+ random=c("a","b"),
type = c("ADFun", "Fun")
)
obj$fn()
obj$gr()
-## Optimize
-system.time(optimize(obj))
-
## Estimate
obj$control <- list(trace=1)
obj$hessian <- TRUE
View
@@ -48,7 +48,7 @@ parameters <- list(
logsdObs=sdObs*0
)
newtonOption(smartsearch=FALSE)
-obj <- MakeADFun(data,parameters,random="^u")
+obj <- MakeADFun(data,parameters,random="u")
obj$hessian <- TRUE
obj$fn()
@@ -21,12 +21,13 @@ x <- as.numeric(A%*%beta+B%*%u+eps)
map <- NULL
obj <- MakeADFun(data=list(x=x,B=B,A=A),
parameters=list(u=u*0+.1,beta=beta*0+.1,logsdu=.1,logsd0=0.1),
- random="^u"
+ random="u"
)
newtonOption(smartsearch=FALSE)
obj$fn(obj$par)
obj$gr(obj$par)
obj$control <- list(parscale=obj$par*0+1e-1,trace=10)
obj$hessian <- TRUE
opt <- do.call("optim",obj)
-summary(as.mle2(opt))
+##summary(as.mle2(opt))
+summary(as.mle(opt))
@@ -3,7 +3,7 @@ dyn.load(dynlib("sumtest"))
parameters=list(x=numeric(5)+0.1 )
df <- list()
-obj <- MakeADFun(data=df,parameters=parameters,random=c("^x"))
+obj <- MakeADFun(data=df,parameters=parameters,random=c("x"))
r <- obj$env$random
par <- obj$env$par
h <- obj$env$spHess(par)[r,r]
View
@@ -43,7 +43,7 @@ parameters <- list(
newtonOption(smartsearch=FALSE)
##compile("rw_parallel.cpp")
dyn.load("rw_parallel.so")
-obj <- MakeADFun(data,parameters,random="^u",DLL="rw_parallel")
+obj <- MakeADFun(data,parameters,random="u",DLL="rw_parallel")
ben <- benchmarkParallel(obj,n=1000,cores=1:8)
ben
pdf("../slides/results/scalability.pdf")
View
@@ -40,7 +40,7 @@ parameters <- list(
logsdObs=sdObs*0
)
newtonOption(smartsearch=FALSE)
-obj <- MakeADFun(data,parameters,random="^u",DLL="mvrw")
+obj <- MakeADFun(data,parameters,random="u",DLL="mvrw")
obj$fn()
obj$gr()
@@ -40,7 +40,7 @@ parameters <- list(
logsdObs=sdObs*0
)
newtonOption(smartsearch=FALSE)
-obj <- MakeADFun(data,parameters,random="^u",DLL="mvrw_sparse")
+obj <- MakeADFun(data,parameters,random="u",DLL="mvrw_sparse")
obj$fn()
obj$gr()
View
@@ -22,7 +22,7 @@ obj <- MakeADFun(data,parameters,map=map,DLL="nmix")
opt <- nlminb(obj$par,obj$fn,obj$gr)
pl <- obj$env$parList(opt$par) ## Parameter estimate after phase 1
## Phase 2
-obj <- MakeADFun(data,pl,random="^u",DLL="nmix")
+obj <- MakeADFun(data,pl,random="u",DLL="nmix")
system.time( opt <- nlminb(obj$par,obj$fn,obj$gr) )
rep <- sdreport(obj)
rep
@@ -15,7 +15,7 @@ obj <- MakeADFun(data=data_orange,
log_sigma=1,
log_sigma_u=2,
u = rep(0,Mmultiply)),
- random=c("^u")
+ random=c("u")
)
obj$control <- list(trace=1,parscale=c(1,1,1,1,1)*1e-2,REPORT=1,reltol=1e-12,maxit=100)
obj$hessian <- F
@@ -71,7 +71,7 @@ parameters <- list(
# The inner problem is quadratic
newtonOption(smartsearch=FALSE)
-obj <- MakeADFun(data,parameters,random=c("^X"),DLL="sde_linear")
+obj <- MakeADFun(data,parameters,random=c("X"),DLL="sde_linear")
# Estimate parameters and latent variables
system.time(opt <- nlminb(obj$par,obj$fn,obj$gr))
View
@@ -15,7 +15,7 @@ obj <- MakeADFun(data=
off_diag_x=rep(0.0,p),
h=t(matrix(0.0,nrow=n,ncol=p))
),
- random=c("^h")
+ random=c("h")
# map=list(off_diag_x=factor(rep(NA,3)))
)
#obj$control <- list(trace=1,parscale=rep(1,12)*1e-2,REPORT=1,reltol=1e-12,maxit=100)
@@ -15,7 +15,7 @@ obj <- MakeADFun(data=
off_diag_x=rep(0.0,p),
h=t(matrix(0.0,nrow=n,ncol=p))
),
- random=c("^h")
+ random=c("h")
# map=list(off_diag_x=factor(rep(NA,3)))
)
#obj$control <- list(trace=1,parscale=rep(1,12)*1e-2,REPORT=1,reltol=1e-12,maxit=100)
View
@@ -22,7 +22,7 @@ parameters <- list(
logsigma=1,
tmpk=rep(0,d$S-1)
)
-obj <- MakeADFun(data,parameters,random="^u",DLL="socatt")
+obj <- MakeADFun(data,parameters,random="u",DLL="socatt")
obj$fn()
obj$gr()
system.time(opt <- do.call("optim",obj))
View
@@ -13,7 +13,7 @@ parameters <- list(
logR=0
)
newtonOption(smartsearch=FALSE)
-obj <- MakeADFun(data,parameters,random="^X",DLL="thetalog")
+obj <- MakeADFun(data,parameters,random="X",DLL="thetalog")
obj$hessian <- TRUE
#obj$control$reltol<-1e-12

0 comments on commit bae04f2

Please sign in to comment.