##' @title Maximum likelihood optimization
##' @param data The angular data to be used for inference
##' @param model A list made of \describe{
##'\item{likelihood}{The likelihood function, see \code{\link{dpairbeta}}
##' for a template}.
##' \item{npar}{The length of the parameter vector}
##' }
##' @param init NULL or a real vector of size \code{model$npar} giving the initial values for \code{link{par}}.
##' @param maxit maximum number of iterations to be performed by
##' function \code{optim}
##' @param method The method to be used by \code{optim}
##' @param hess logical: should an approximation of the hessian be performed ?
##' @param link the link function from the natural marginal parameter spaces to the real line.
##' @param unlink the inverse link function. If \code{x} is any real number, then \code{unlink(x)} should be in the admissible range for the likelihood function and the prior function.
##' @return The list returned by \code{optim} and the AIC and BIC criteria
##' @export
maxLikelihood <-
function(data, model,## = list(likelihood, npar),
init = NULL, maxit = 500, method="L-BFGS-B", hess = T, link, unlink)
## @param prior The prior function (see \code{\link{prior.pb}} for a template) for generating initial parameters in case the initial value results in non finite log-likelihood.
## @param Hpar the prior hyper parameters
{
p = dim(data)[2]
ndat=dim(data)[1]
npar <- model$npar
lhood <- function(vec)
{
-model$likelihood(x=data, par=unlink(vec),
log=TRUE, vectorial=FALSE)
}
if(is.null(init))
init <- rep(0,model$npar)
count <- 0
converged <- FALSE
while(!converged & count<20 )
{
count <- count+1
opt.test <- tryCatch(optim(init, lhood, method = method,
control =list(maxit = maxit,trace=0 ),
hessian = hess),
error=function(e){
# print(e)
return(list(convergence=100))
})
if(as.integer(opt.test$convergence)>0){
init <- rnorm(model$npar,mean=0,sd=1)
}
else{
converged <- TRUE
opt <- opt.test
}
}
rtn <- list()
if(as.integer(opt.test$convergence)>0){
rtn$message <- "optimisation failed"
rtn$counts <- count
rtn$convergence <- 100
rtn$linkedpar <- rep(0,npar)
rtn$par <- unlink(rep(0,npar))
rtn$value <- -Inf
rtn$aic <- rtn$aicc <- rtn$bic <- Inf
rtn$linkedHessian <- 1
return(rtn)
}
rtn$message <- opt$message
rtn$counts <- opt$counts
rtn$convergence <- opt$convergence
rtn$linkedpar <- opt$par
rtn$par <- unlink(opt$par)
rtn$value <- opt$value
rtn$aic <- 2* (opt$value + npar)
rtn$aicc <- 2* (opt$value + npar) +2*npar*(npar+1)/(ndat-npar-1 )
rtn$bic <- 2*opt$value + npar*log(ndat)
if(!hess)
{
rtn$linkedHessian <- 1
}
else
{
rtn$linkedHessian <- opt$hessian
tryCatch(expr={
asympt.variance <- chol2inv(chol(opt$hessian))
rtn$asympt.variance <- asympt.variance
rtn$linked.esterr <- sqrt(diag(asympt.variance))},
error=function(e){## rtn$asympt.variance <- NULL
## rtn$linked.esterr <- NULL
return(rtn)} )
}
return(rtn)
}
##?optim
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.