Nothing
#'Continuous multi-dimensional optimization
#'@references Yashin, A.I. et al (2007). Stochastic model for analysis of longitudinal data on aging
#'and mortality. Mathematical Biosciences, 208(2), 538-551.<DOI:10.1016/j.mbs.2006.11.006>.
#'@param dat A data table.
#'@param a A starting value of the rate of adaptive response to any deviation of Y from f1(t).
#'@param f1 A starting value of the average age trajectories of the variables which process is forced to follow.
#'@param Q Starting values of the quadratic hazard term.
#'@param f A starting value of the "optimal" value of variable which corresponds to the minimum of hazard rate at a respective time.
#'@param b A starting value of a diffusion coefficient representing a strength of the random disturbance from Wiener Process.
#'@param mu0 A starting value of the baseline hazard.
#'@param theta A starting value of the parameter theta (axe displacement of Gompertz function).
#'@param verbose An indicator of verbosing output.
#'@param stopifbound Estimation stops if at least one parameter achieves lower or upper boundaries.
#'#'Check the NLopt website for a description of
#'the algorithms. Default: NLOPT_LN_NELDERMEAD
#'@param lb Lower bound of parameters under estimation.
#'@param ub Upper bound of parameters under estimation.
#'The program stops when the number of function evaluations exceeds maxeval. Default: 500.
#'@param pinv.tol A tolerance value for pseudo-inverse of matrix gamma (see Yashin, A.I. et al (2007). Stochastic model for analysis of longitudinal data on aging
#'and mortality. Mathematical Biosciences, 208(2), 538-551.<DOI:10.1016/j.mbs.2006.11.006>.)
#'@param gomp A flag (FALSE by default). When it is set, then time-dependent exponential form of mu0 is used:
#' mu0 = mu0*exp(theta*t).
#'@param opts A list of options for \code{nloptr}.
#'Default value: \code{opt=list(algorithm="NLOPT_LN_NELDERMEAD",
#'maxeval=100, ftol_rel=1e-8)}.
#'Please see \code{nloptr} documentation for more information.
#'@param logmu0 Natural logarith of baseline mortality.
#'Default: \code{FALSE}.
#'@return A set of estimated parameters a, f1, Q, f, b, mu0, theta and
#'additional variable \code{limit} which indicates if any parameter
#'achieved lower or upper boundary conditions (FALSE by default).
#'@return status Optimization status (see documentation for nloptr package).
#'@return LogLik A logarithm likelihood.
#'@return objective A value of objective function (given by nloptr).
#'@return message A message given by nloptr optimization function (see documentation for nloptr package).
#'@details \code{spm_continuous} runs much slower that discrete but more precise and can handle time
#'intervals with different lengths.
#'@export
#'@examples
#'library(stpm)
#'set.seed(123)
#'#Reading the data:
#'data <- simdata_cont(N=2)
#'head(data)
#'#Parameters estimation:
#'pars <- spm_continuous(dat=data,a=-0.05, f1=80,
#' Q=2e-8, f=80, b=5, mu0=2e-5)
#'pars
#'
spm_continuous <- function(dat,
a=-0.05,
f1=80,
Q=2e-8,
f=80,
b=5,
mu0=2e-5,
theta=0.08,
stopifbound=FALSE,
lb=NULL, ub=NULL,
verbose=FALSE,
pinv.tol=0.01,
gomp=FALSE,
opts=list(algorithm="NLOPT_LN_NELDERMEAD",
maxeval=100, ftol_rel=1e-8),
logmu0=FALSE) {
setlb <- function(k, params) {
# This function sets lower and upper boundaries for optim.
# - k - number of dimensions
# - params - initial parameters, a vector
#
# Lower boundaries:
lower_bound <- c()
# Setting boundaries for coefficients:
# aH
start=1; end=k^2
lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) })))
# f1H
start=end+1; end=start+k-1
lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) })))
# QH
start=end+1; end=start+k^2-1
lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ ifelse(params[n] > 0, 0, params[n]+0.1*params[n]) })))
# fH
start=end+1; end=start+k-1
lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) })) )
# bH
start=end+1; end=start+k-1
lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) })) )
# mu0
start=end+1; end=start
lower_bound <- c( lower_bound, params[start:end] + ifelse(params[start:end] > 0, -0.1*params[start:end], 0.1*params[start:end]))
# theta
start=end+1; end=start
lower_bound <- c( lower_bound, params[start:end] - 0.1*params[start:end])
lower_bound
}
setub <- function(k, params) {
# This function sets lower and upper boundaries for optim.
# - k - number of dimensions
# - params - initial parameters, a vector
#
#Upper boundaries:
upper_bound <- c()
# Setting boundaries for coefficients:
# aH
start=1; end=k^2
upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n] + 0.1*params[n], 0*params[n]) })))
# f1H
start=end+1; end=start+k-1
upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n] + 0.1*params[n], params[n]- 0.1*params[n]) })))
# QH
start=end+1; end=start+k^2-1
upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n) {ifelse(params[n] > 0, params[n] + 0.1*params[n], params[n]-0.1*params[n] )})))
# fH
start=end+1; end=start+k-1
upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n]+0.1*params[n], params[n]-0.1*params[n]) })))
# bH
start=end+1; end=start+k-1
upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n]+0.1*params[n], params[n]-0.1*params[n]) })))
# mu0
start=end+1; end=start
upper_bound <- c( upper_bound, ifelse(params[start:end] > 0, params[start:end]+0.1*params[start:end], params[start:end]-0.1*params[start:end]))
# theta
start=end+1; end=start
upper_bound <- c( upper_bound, params[start:end] + 0.1*params[start:end])
upper_bound
}
###=======For DEBUG========###
#dat = dat
#
#a=matrix(c(-0.05, 0.001, 0.001, -0.05), nrow = 2, ncol = 2, byrow = T)
#f1=t(matrix(c(100, 200), nrow = 2, ncol = 1, byrow = F))
#Q=matrix(c(1e-06, 1e-7, 1e-7, 1e-06), nrow = 2, ncol = 2, byrow = T)
#f=t(matrix(c(100, 200), nrow = 2, ncol = 1, byrow = F))
#b=matrix(c(2, 5), nrow = 2, ncol = 1, byrow = F)
#mu0=1e-4
#theta=0.08
#k=2
#
#stopifbound=FALSE
#algorithm="NLOPT_LN_NELDERMEAD"
#lb=NULL
#ub=NULL
#verbose=FALSE
###=========================###
#avail_algorithms <- c("NLOPT_GN_DIRECT", "NLOPT_GN_DIRECT_L",
# "NLOPT_GN_DIRECT_L_RAND", "NLOPT_GN_DIRECT_NOSCAL",
# "NLOPT_GN_DIRECT_L_NOSCAL",
# "NLOPT_GN_DIRECT_L_RAND_NOSCAL",
# "NLOPT_GN_ORIG_DIRECT", "NLOPT_GN_ORIG_DIRECT_L",
# "NLOPT_GD_STOGO", "NLOPT_GD_STOGO_RAND",
# "NLOPT_LD_SLSQP", "NLOPT_LD_LBFGS_NOCEDAL",
# "NLOPT_LD_LBFGS", "NLOPT_LN_PRAXIS", "NLOPT_LD_VAR1",
# "NLOPT_LD_VAR2", "NLOPT_LD_TNEWTON",
# "NLOPT_LD_TNEWTON_RESTART",
# "NLOPT_LD_TNEWTON_PRECOND",
# "NLOPT_LD_TNEWTON_PRECOND_RESTART",
# "NLOPT_GN_CRS2_LM", "NLOPT_GN_MLSL", "NLOPT_GD_MLSL",
# "NLOPT_GN_MLSL_LDS", "NLOPT_GD_MLSL_LDS",
# "NLOPT_LD_MMA", "NLOPT_LN_COBYLA", "NLOPT_LN_NEWUOA",
# "NLOPT_LN_NEWUOA_BOUND", "NLOPT_LN_NELDERMEAD",
# "NLOPT_LN_SBPLX", "NLOPT_LN_AUGLAG", "NLOPT_LD_AUGLAG",
# "NLOPT_LN_AUGLAG_EQ", "NLOPT_LD_AUGLAG_EQ",
# "NLOPT_LN_BOBYQA", "NLOPT_GN_ISRES")
#
#if(!(algorithm %in% avail_algorithms)) {
# stop(cat("Provided algorithm", algorithm, "not in the list of available optimization methods."))
#}
dat <- as.matrix(dat[, 2:dim(dat)[2]])
k <- dim(as.matrix(a))[1]
final_res <- list()
if(mu0 < 0) {mu0 <- 0}
parameters <- c(t(a), f1, t(Q), f, b, mu0, theta)
# Current results:
#results <- list(a=NULL, f1=NULL, Q=NULL, f=NULL, b=NULL, mu0=NULL, theta=NULL)
results_tmp <- list(a=NULL, f1=NULL, Q=NULL, f=NULL, b=NULL, mu0=NULL, theta=NULL)
e <- new.env()
iteration <- 0
bounds <- list()
if(is.null(lb)) {
bounds$lower_bound <- setlb(k, parameters)
} else {
bounds$lower_bound <- lb
}
if(is.null(ub)) {
bounds$upper_bound <- setub(k, parameters)
} else {
bounds$upper_bound <- ub
}
# Reading parameters:
#start=1; end=k^2
#a <- matrix(parameters[start:end],ncol=k, nrow=k, byrow=T)
#results$a <- a
##print(results$a)
#start=end+1; end=start+k-1
#f1 <- matrix(parameters[start:end],ncol=1, nrow=k, byrow=T)
#results$f1 <- f1
##print(results$f1)
#start=end+1; end=start+k^2-1
#Q <- matrix(parameters[start:end],ncol=k, nrow=k, byrow=T)
#results$Q <- Q
##print(results$Q)
#start=end+1; end=start+k-1
#f <- matrix(parameters[start:end],ncol=1, nrow=k, byrow=F)
#results$f <- f
##print(results$f)
#start=end+1; end=start+k-1
#b <- matrix(parameters[start:end],nrow=k, ncol=1, byrow=F)
#results$b <- b
##print(results$b)
#start=end+1; end=start
#mu0 <- parameters[start:end]
#results$mu0 <- mu0
##print(results$mu0)
#start=end+1; end=start
#theta <- parameters[start:end]
#results$theta <- theta
##print(results$theta)
## End of reading parameters
L.prev <- NA
maxlik <- function(par) {
stopflag <- F
# Reading parameters:
start=1; end=k^2
a <- matrix(par[start:end],ncol=k, , nrow=k, byrow=TRUE)
results_tmp$a <<- a
start=end+1; end=start+k-1
f1 <- matrix(par[start:end],ncol=1, nrow=k, byrow=FALSE)
results_tmp$f1 <<- f1
start=end+1; end=start+k^2-1
Q <- matrix(par[start:end],ncol=k, nrow=k, byrow=TRUE)
results_tmp$Q <<- Q
start=end+1; end=start+k-1
f <- matrix(par[start:end],ncol=1, nrow=k, byrow=FALSE)
results_tmp$f <<- f
start=end+1; end=start+k-1
b <- matrix(par[start:end],nrow=k, ncol=1, byrow=FALSE)
results_tmp$b <<- b
start=end+1; end=start
mu0 <- par[start:end]
results_tmp$mu0 <<- mu0
start=end+1; end=start
theta <- par[start:end]
results_tmp$theta <<- theta
# End reading parameters
if(stopifbound) {
for(i in 1:length(results_tmp)) {
if(length(intersect(results_tmp[[i]],c(bounds$lower_bound[i], bounds$upper_bound[i]))) >= 1) {
cat("Parameter", names(results_tmp)[i], "achieved lower/upper bound. Process stopped.\n")
cat(results_tmp[[i]],"\n")
stopflag <- T
break
}
}
}
if(stopflag == FALSE) {
dims <- dim(dat)
res <- .Call("complikMD", dat, dims[1], dims[2], a, f1, Q, b, f, mu0, theta, k, pinv.tol, gomp, logmu0)
assign("results", results_tmp, envir=e)
iteration <<- iteration + 1
L.prev <<- res
if(verbose) {
cat("L = ", res,"\n")
cat("Iteration: ", iteration, "\nResults:\n")
print(results_tmp)
}
}
return(as.numeric(-1*res))
}
# Optimization:
if(verbose) {
cat("Lower bound:\n")
print(bounds$lower_bound)
cat("Upper bound:\n")
print(bounds$upper_bound)
}
#tryCatch({ans <- nloptr(x0 = parameters,
# eval_f = maxlik, opts = list("algorithm"=algorithm,
# "xtol_rel"=1.0e-1, "maxeval"=maxeval),
# lb = bounds$lower_bound, ub = bounds$upper_bound)
#
# },
# error=function(e) {if(verbose == TRUE) {print(e)}},
# finally=NA)
#print(parameters)
#print(bounds$lower_bound)
#print(bounds$upper_bound)
ans <- nloptr(x0 = parameters,
eval_f = maxlik,
opts = opts,
lb = bounds$lower_bound, ub = bounds$upper_bound)
final_results <- get("results",envir=e)
#final_results <- results_tmp
final_results[["status"]] <- ans$status
final_results[["LogLik"]] <- L.prev
final_results[["objective"]] <- ans$objective
final_results[["message"]] <- ans$message
# Check if any parameter achieved upper/lower limit and report it:
limit <- FALSE
for(i in 1:length(final_results)) {
if(length(intersect(final_results[[i]],c(bounds$lower_bound[i], bounds$upper_bound[i]))) >= 1) {
cat("Parameter", names(final_results)[i], "achieved lower/upper bound.\n")
cat(final_results[[i]],"\n")
limit <- TRUE
}
}
if(logmu0) {
final_results[["mu0"]] <- exp(final_results[["mu0"]])
}
final_results$limit <- limit
#assign("results", final_results, envir=baseenv())
class(final_results) <- "spm.continuous"
invisible(final_results)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.