vam.R
In rcqls/VAM: Virtual Age Models
# Simulation: sim.vam
sim.vam <- function(formula,data.covariates) {
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(sim.vam,formula=formula(formula),data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
## specify names of variables inside data.frame
if(is.null(model$response)) {
self$response.names <- c("Time","Type")
self$system.name <- "System"
} else {
if(length(model$response)==3) {
self$system.name <- model$response[1]
self$response.names <- tail(model$response,2)
} else if(length(model$response)==2) {
self$system.name <- "System"
self$response.names <- model$response
}
}
self$model<-model#debuggage LD
self$formula <- substitute.vam.formula(model=model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
rcpp <- new(SimVam,model)
rcpp
})
self
}
# TODO: when data provided, complete the data!
simulate.sim.vam <- function(sim, stop.policy, nb.system, cache.size=500,as.list=FALSE,data=NULL) {
# To have a first argument more readable
self <- sim
rcpp <- self$rcpp()
if(missing(stop.policy)) stop.policy <- if(is.null(sim$model.parsed$covariates)) 10 else 1
if(missing(nb.system)) nb.system <- if(is.null(sim$model.parsed$covariates)) 1 else nrow(sim$model.parsed$covariates$data)
self$stop.policy.last <- parse.stop.policy(deparse(substitute(stop.policy)))
stop.policy <- eval(self$stop.policy.last,parent.frame(1))
if(is.numeric(stop.policy)) {
if(stop.policy == as.integer(stop.policy)) {#integer
stop.policy <- EndAt(size=stop.policy)
} else stop.policy <- NULL
} else if(!inherits(stop.policy,"stop.policy")) {
stop.policy <- NULL
}
if(is.null(stop.policy)) {
warning("Argument stop.policy is not a proper one!")
return(invisible(NULL))
}
# default cache size
if(is.null(stop.policy$cache.size)) stop.policy$cache.size <- cache.size
# add stop.policy object
rcpp$add_stop_policy(stop.policy)
# verify the coherence of the data argument
if(!is.null(data)){
if(!all(sim$response.names%in%names(data))) {
warning("The column names of the data argument must contain the response name (if not precise in the R formula those are by default Time and Type)!")
return(invisible(NULL))
}
data<-data[data[,sim$response.names[2]]!=0,]
}
if(nb.system>1) {
# multisystem
if(as.list) df<-list()
for(i in 1:nb.system) {
rcpp$select_current_system(i-1) # i-1 because it is c++
dataTime<-c()
dataType<-c()
if(!is.null(data)){
if(sim$system.name%in%names(data)){
dataInit<-data[data[,sim$system.name]==i,sim$response.names]
dataTime<-dataInit[,sim$response.names[1]]
dataType<-dataInit[,sim$response.names[2]]
} else {
dataTime<-data[,sim$response.names[1]]
dataType<-data[,sim$response.names[2]]
}
}
df2 <- rcpp$simulate(stop.policy$cache.size,as.numeric(dataTime),as.numeric(dataType))[-1,]
names(df2) <- sim$response.names
if(as.list) {
rownames(df2)<-1:nrow(df2)
df[[i]] <- df2 #rbind(data.frame(Time=0,Type=1),df2)
} else {
df2[[sim$system.name]] <- i
df2<-df2[c(3,1:2)]
df <- if(i==1) df2 else rbind(df,df2)
}
}
} else {
rcpp$select_current_system(0) # added for covariates case
df <- rcpp$simulate(stop.policy$cache.size,as.numeric(data[,sim$response.names[1]]),as.numeric(data[,sim$response.names[2]]))[-1,]
names(df) <- sim$response.names
if(as.list) {rownames(df)<-1:nrow(df);df<-list(df)}
}
if(!as.list) {
rownames(df) <- 1:nrow(df)
rcpp$set_data(data.frame.to.list.multi.vam(df,names(df))) #if already list, response not used inside data.frame.to.list.multi.vam
}
else {
names(df) <- paste0(sim$system.name,1:length(df))
rcpp$set_data(unname(df))
}
## put the final result transformed as in mle.vam and model.vam to the model
## return the result as a data.frame
df
}
# Model part
model.vam <- function(formula,data,data.covariates) {
if(missing(data)) data<-NULL
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(model.vam,formula=formula(formula),data=data,data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
self$formula <- substitute.vam.formula(model=model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
if(is.null(self$data)) {## No data
rcpp <- new(ModelVam,model)
rcpp
} else {
response <- model$response
data <- data.frame.to.list.multi.vam(self$data,response)
rcpp <- new(ModelVam,model,data)
rcpp
}
})
self
}
update.model.vam <- function(model,data) {
if(!missing(data)) {
self <- model #to have an argument more readable
response <- parse.vam.formula(self$formula)$response
self$data <- data
data2 <- data.frame.to.list.multi.vam(self$data,response)
self$rcpp()$set_data(data2)
}
}
check.censorship <- function(data) {
lapply(data,function(data.syst) {
ty <- data.syst$Type[-1]
check<-list(type="none",index=0)
if(any(ty==0)) {
if(which(ty==0)[1]==length(ty)) check$type <- 'right' #meaning: right only
else if(!any(ty<0) || (which(ty<0)[1] > which(ty==0)[1])) {
check$type <- 'left' #can have also a right censorship
check$index <- which(ty==0)[1]
}
else check$type <- 'unknown'
}
check
})
}
data.frame.to.list.multi.vam <- function(data,response) {
# return data if it is already only a list!
if(is.list(data) && !is.data.frame(data)) {
data.frame(Time=0,Type=1)->tmp
names(tmp) <- names(data[[1]])
return(lapply(data,function(df) rbind(tmp,df)))
}
if(!("System" %in% response) && ("System" %in% names(data)) ) {
warning(paste0("WARNING: data has variable 'System' when response in formula does not contain this variable!"))
}
# otherwise
if(length(response)==2) {
if(length(intersect(response,names(data))) != 2) stop(paste0("Bad response:",response))
tmp <- data[[response[1]]]
data2 <- list(data.frame(Time=c(0,tmp[order(tmp)]),Type=c(1,data[[response[2]]][order(tmp)])))
} else {
if(length(intersect(response,names(data))) != 3) stop(paste0("Bad response:",response))
syst0 <- unique(syst<-data[[response[1]]])
data2 <- list()
for(i in seq_along(syst0)) {
df <- data[syst==syst0[i],response]
tmp <- df[[response[2]]]
data2[[i]] <- data.frame(Time=c(0,tmp[order(tmp)]),Type=c(1,df[[response[3]]][order(tmp)]))
}
}
data2
}
# TODO: check data
check.data.vam <-function(data,response) {
if(all(data[[response[[-length(response)]]]])) {
}
}
make.censorship <- function(data,rcpp) {
if(any(sapply(censorship <- check.censorship(data),function(e) e$type == 'unknown'))) stop("VAM package does know how to deal with this kind of censorship!")
if(any(sapply(censorship,function(e) e$type == 'left'))) {
leftCensors <- sapply(censorship,function(e) e$index)
rcpp$set_leftCensors(leftCensors)
}
}
mle.vam <- function(formula,data,data.covariates) {
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(mle.vam,formula=formula(formula),data=data,data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
self$formula <- substitute.vam.formula(model=model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
response <- model$response
data <- data.frame.to.list.multi.vam(self$data,response)
rcpp <- new(MLEVam,model,data)
make.censorship(data,rcpp)
rcpp
})
self
}
# to convert in Rcpp
params.model.vam <- params.sim.vam <- params.mle.vam <- params.bayesian.vam <- function(self,param) {
if(missing(param)) {
self$rcpp()$get_params()
} else {
self$rcpp()$set_params(param)
}
}
## Useless since stats:::formula.default do that by default
# formula.model.vam <- formula.sim.vam <- function(self) {
# self$formula
# }
formula.mle.vam <- function(self,origin=FALSE) {
form <- substitute.vam.formula(self$formula,coef(self))
if(origin) list(formula=form,origin=self$formula)
else form
}
update.mle.vam <- function(mle,data) {
if(!missing(data)) {
self <- mle
model <- parse.vam.formula(self$formula)
response <- model$response
self$data <- data
data2 <- data.frame.to.list.multi.vam(self$data,response)
self$rcpp()$set_data(data2)
self$rcpp()$reset_leftCensors()
make.censorship(data2,self$rcpp())
## estimation has to be computed again!
self$mle.coef<-NULL
}
## with
}
#fonction de LD2
contrast.mle.vam <-function(obj,par0,with_value=TRUE,with_gradient=FALSE,with_hessian=FALSE){
type <- c(with_value,with_gradient,with_hessian)
rcpp <- obj$rcpp()
## save the initial param
if(is.null(obj$par0)) obj$par0 <- params(obj)
## parameters stuff!
if(missing(par0)) {
if("par" %in% names(obj)) {
param <- obj$par[-1]
alpha <- obj$par[1]#LD2
} #not the first run
else {
param<-params(obj)[-1] #first run
alpha<-params(obj)[1]#LD2
}
} else if(is.null(par0)) {
param<-obj$par0[-1]
alpha<-obj$par0[1]#LD2
} else {
param<-par0[-1]
alpha<-par0[1]#LD2
}
if(sum(type)==1) {
if(type[1]){
res<-rcpp$contrast(c(alpha,param),FALSE)
} else if(type[2]){
res<-rcpp$gradient(c(alpha,param),FALSE)[-1]
} else if(type[3]){
res<-rcpp$hessian(c(alpha,param),FALSE)[-1,-1]
}
} else {
res<-list()
if (type[1]){
res$contrast<-rcpp$contrast(c(alpha,param),FALSE)
}
if(type[2]){
res$gradient<-rcpp$gradient(c(alpha,param),FALSE)[-1]
}
if(type[3]){
res$hessian<-rcpp$hessian(c(alpha,param),FALSE)[-1,-1]
}
}
res
}
#fonction de LD2
logLik.mle.vam <-function(obj,par0,with_value=TRUE,with_gradient=FALSE,with_hessian=FALSE){
type=c(with_value,with_gradient,with_hessian)
rcpp <- obj$rcpp()
## save the initial param
if(is.null(obj$par0)) obj$par0 <- params(obj)
## parameters stuff!
if(missing(par0)) {
if("par" %in% names(obj)) {
param <- obj$par[-1]
alpha <- obj$par[1]#LD2
} #not the first run
else {
param<-params(obj)[-1] #first run
alpha<-params(obj)[1]#LD2
}
} else if(is.null(par0)) {
param<-obj$par0[-1]
alpha<-obj$par0[1]#LD2
} else {
param<-par0[-1]
alpha<-par0[1]#LD2
}
if(sum(type)==1) {
if(type[1]){
res<-rcpp$contrast(c(alpha,param),TRUE)
} else if(type[2]){
res<-rcpp$gradient(c(alpha,param),TRUE)
} else if(type[3]){
res<-rcpp$hessian(c(alpha,param),TRUE)
}
} else {
res<-list()
if (type[1]){
res$contrast<-rcpp$contrast(c(alpha,param),TRUE)
}
if(type[2]){
res$gradient<-rcpp$gradient(c(alpha,param),TRUE)
}
if(type[3]){
res$hessian<-rcpp$hessian(c(alpha,param),TRUE)
}
}
res
}
# alpha is not considered in the estimation!
run.mle.vam <-function(obj,par0,fixed,method=NULL,verbose=FALSE,...) {
rcpp <- obj$rcpp()
par0.tmp <- init.par0(obj,par0)
param <- par0.tmp$param
alpha <- par0.tmp$alpha
# ## save the initial param
# if(is.null(obj$par0)) obj$par0 <- params(obj)
# ## parameters stuff!
# if(missing(par0)) {
# if("par" %in% names(obj)) {
# param <- obj$par[-1]
# alpha <- obj$par[1] #LD2
# } else {#not the first run
# param<-params(obj)[-1] #first run
# alpha<-params(obj)[1] #LD2
# }
# } else if(is.null(par0)) {
# param<-obj$par0[-1]
# alpha<-obj$par0[1] #LD2
# } else {
# param<-par0[-1]
# alpha<-par0[1] #LD2
# }
## fixed and functions stuff!
fixed.tmp <- init.fixed.param(param,fixed)
fixed <- fixed.tmp$fixed
alpha_fixed <- fixed.tmp$alpha_fixed
fn<-function(par) {
##cat("param->");print(par);print(param[!fixed])
param[!fixed]<-par
#cat("param->");print(param)
## All the commented part allows us to save the param when value is NaN
#res<-
-rcpp$contrast(c(alpha,param),alpha_fixed) #LD2
# if(is.nan(res)) {
# mode_param<-"contrast"
# data_param<- rcpp$get_data(0)
# data_param2 <- obj$data
# save(param,mode_param,data_param,data_param2,file="/Users/remy/tmp/VAM/res.RData")
# }
# res
}
gr <- function(par) {
param[!fixed]<-par
#cat("param2->");print(param)
## All the commented part allows us to save the param when value is NaN
#res <-
(-rcpp$gradient(c(alpha,param),alpha_fixed)[-1])[!fixed] #LD2
# if(any(is.nan(res))) {
# mode_param<-"gradient"
# save(param,mode_param,file="/Users/remy/tmp/VAM/res.RData")
# }
# res
}
## optim stuff!
if(is.null(method) || method=="fast") {
if(length(param[!fixed])>1) param[!fixed]<-(res <- optim(param[!fixed],fn,gr,method="Ne",...))$par
if(is.null(method)) res<-optim(param[!fixed],fn,gr,method="CG",...)
} else {
res<-optim(param[!fixed],fn,gr,method=method,...)
}
#fixed tips
param[!fixed]<-res$par
if(!alpha_fixed){#LD2
## complete the scale parameter
alpha <- rcpp$alpha_est(c(1,param)) #LD2
} #LD2
res$par<-c(alpha,param) #LD2
if(verbose) print(res)
## save stuff
obj$fixed <- c(alpha_fixed,fixed) #LD2
obj$optim<-res
obj$par<-res$par
obj$mle.coef <- res$par
params(obj,obj$mle.coef) #put the result in the c++ part
##obj$mle.coef
obj$optim
}
## Rmk: run.mle.vam is supposed to run many times to get the best estimate!
## Here, par=NULL forces initialisation update but does not ensure that it is the best estimate.
## TODO: try to find a best strategy or many strategies...
## ... added to deal with fixed by example!
coef.mle.vam <- function(obj,par=NULL,method=NULL,verbose=FALSE,...) {
if(is.null(obj$mle.coef) || !is.null(par)) {
res <-run.mle.vam(obj,par,verbose=verbose,method=method,...)
if(verbose && obj$optim$convergence>0) cat("convergence=",obj$optim$convergence,"\n",sep="")
}
obj$mle.coef
}
# bayesian.vam class
bayesian.vam <- function(formula,data,data.covariates) {
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(bayesian.vam,formula=formula,data=data,data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
self$formula <- substitute.vam.formula(model=model)
response <- model$response
data <- data.frame.to.list.multi.vam(self$data,response)
self$priors <- priors.from.vam.formula(model)
##DEBUG: print("priors");print(priors)
##DEBUG: print("modelAV");print(model)
self$prior.params <- sapply(self$priors,mean)
self$mle.formula <- substitute.vam.formula(self$formula,self$prior.params)
## THIS IS LESS CLEVER THAN THE NEXT LINE: print(model<-bayesian.model.to.mle.model(model,priors))
model<-parse.vam.formula(self$mle.formula)
##DEBUG: print("modelAP");print(model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
rcpp <- new(BayesianVam,model,data,self$priors)
rcpp
})
self
}
run.bayesian.vam <- function(obj,par0,fixed,sigma.proposal,nb=100000,burn=10000,profile.alpha=FALSE,method=NULL,verbose=FALSE,history=FALSE,proposal='norm',...) {
rcpp <- obj$rcpp()
obj$history<-history
obj$nb<-nb
obj$burn<-burn
obj$preplots <- NULL
## init via mle: par0 is supposed first to be initialized by mle
if(missing(par0)) {
if(is.null(obj$model.parsed$covariates)) obj$mle <- mle.vam(obj$mle.formula,obj$data) else obj$mle <- mle.vam(obj$mle.formula,obj$data,obj$data.covariates)
obj$mle.init <- TRUE
obj$par0 <- coef(obj$mle,fixed=fixed,method=method,verbose=verbose,...)
} else {
obj$par0<-par0
obj$mle.init <- FALSE
}
fixed.tmp <- init.fixed.param(obj$par0[-1],fixed)
obj$fixed <- c(fixed.tmp$alpha_fixed,fixed.tmp$fixed)
#obj$alpha_fixed <- fixed.tmp$alpha_fixed
obj$profile_alpha<-profile.alpha
if (fixed.tmp$alpha_fixed & profile.alpha){
warning("Parameter alpha can not simultaneously be fixed and optimized: it is fixed !")
obj$profile_alpha<-FALSE
}
##print(obj$mle.init)
if(missing(sigma.proposal)) sigma.proposal <- sapply(obj$priors,sigma)
else {
if(length(sigma.proposal)==1) sigma.proposal <- rep(sigma.proposal,length(obj$priors))
}
if(length(proposal)==1) proposal <- rep(proposal,length(obj$priors))
for(i in (1:length(obj$priors))) {
rcpp$set_sigma(i-1,sigma.proposal[i])
rcpp$set_proposal(i-1,switch(proposal[i],'lnorm'=1,0))
}
obj$sigma_proposal<-sigma.proposal
obj$proposal<-proposal
if(history) {
res <- rcpp$mcmc_history(obj$par0,nb,burn,obj$fixed,obj$profile_alpha)
obj$nb<-res[[3+obj$profile_alpha]]
obj$nb_proposal<-nb
obj$par<-as.data.frame(res[1:(2+obj$profile_alpha)])
names(obj$par) <- c("ind","estimate","alpha")[1:(2+obj$profile_alpha)]
} else {
obj$par <- rcpp$mcmc(obj$par0,nb,burn,obj$fixed,obj$profile_alpha)
}
obj$par
}
coef.bayesian.vam <- function(obj,new.run=FALSE,...) {
if(new.run || is.null(obj$par)) run(obj,...)
if(obj$history){
param<-sapply(obj$profile_alpha:(length(obj$par0)-1),function(j){mean(obj$par$estimate[obj$par$ind==j])})
if(obj$profile_alpha){
## complete the scale parameter
param <- c(mean(obj$par$alpha),param)
}
} else {
param <- sapply(obj$par,mean)
}
param[as.logical(obj$fixed)]<-obj$par0[as.logical(obj$fixed)]
param
}
hist.bayesian.vam <- function(obj,i=1,...) {
if(is.null(obj$par)) run(obj)
if(obj$fixed[i]){
warning("No hist for a fixed parameter!")
} else {
if(obj$history){
if((obj$profile_alpha)&(i==1)){
thetak<-obj$par$alpha
} else {
thetak<-obj$par$estimate[obj$par$ind==(i-1)]
}
hist(thetak,prob=TRUE,...)
abline(v=mean(thetak),col="blue",lwd=2)
} else{
hist(obj$par[[i]],prob=TRUE,...)
abline(v=mean(obj$par[[i]]),col="blue",lwd=2)
}
}
}
summary.bayesian.vam <- function(obj,alpha=0.05,new.run=FALSE,digits=4,...) {
if(new.run || is.null(obj$par)) run(obj,...)
res1<-obj$par0
cat("Initial parameters",if(obj$mle.init) " (by MLE)" else "",": ",paste(signif(res1,digits=digits),collapse=", "),"\n",sep="")
res2<-coef(obj)
cat("(Mean) Bayesian estimates: ", paste(signif(res2,digits=digits),collapse=", "),"\n",sep="")
if(obj$history){
thetak<-lapply(obj$profile_alpha:(length(obj$par0)-1),function(j){obj$par$estimate[obj$par$ind==j]})
if(obj$profile_alpha){thetak<-c(list(obj$par$alpha),thetak)}
} else {
thetak<-obj$par
}
sd<-sapply(thetak,sd)
sd[as.logical(obj$fixed)]<-0
cat("(SD) Bayesian estimates: ",paste(signif(sd,digits=digits),collapse=", "),"\n",sep="")
res4<-sapply(thetak,function(x){quantile(x,probs=alpha/2)})
cat("(", alpha/2,"-Quantile) Bayesian estimates: ",paste(signif(res4,digits=digits),collapse=", "),"\n",sep="")
res5<-sapply(thetak,function(x){quantile(x,probs=1-alpha/2)})
cat("(", 1-alpha/2,"-Quantile) Bayesian estimates: ",paste(signif(res5,digits=digits),collapse=", "),"\n",sep="")
res6<-sapply(thetak,length)
cat("(Number) Bayesian estimates: ",paste(res6,collapse=", "),"\n",sep="")
res7<-sapply(thetak,length)/(obj$nb-obj$burn)
cat("Metropolis-Hasting acceptation rates: ",paste(signif(res7,digits=digits),collapse=", "),"\n",sep="")
res<-data.frame(res1,res2,sd,res4,res5,res6,res7)
names(res)<-c(paste("Init",if(obj$mle.init) "(MLE)" else "",sep=""),"Mean","SD",paste(alpha/2,"-Quantile",sep=""),paste(1-alpha/2,"-Quantile",sep=""),"Number","Accept_Rate")
invisible(res)
}
# for both sim and mle
parse.vam.formula <- function(formula) {
## Needs to have this envir to evaluate params (otherwise, beta was found in baseenv() first before globalenv() for example!)
envir.eval <- parent.frame(4)
eval.vam <- function(e) eval(e,envir.eval)
if(formula[[1]] != as.name("~")) stop("Argument has to be a formula")
if(length(formula) == 2) {
response <- NULL
cm <- formula[[2]]
} else {
tmp <- formula[[2]]
## simplify parenthesis
while(tmp[[1]] == as.name("(")) tmp <- tmp[[2]]
if(tmp[[1]] != as.name("&") && length(tmp) != 3) stop("Left part of formula of the form 'Time & Type'!")
if(length(tmp[[2]])==3 && tmp[[2]][[1]]==as.name("&")) {
response <- c(as.character(tmp[[2]][[2]]),as.character(tmp[[2]][[3]]),as.character(tmp[[3]]))
} else response <- c(as.character(tmp[[2]]),as.character(tmp[[3]]))
cm <- formula[[3]]
}
## simplify parenthesis
while(cm[[1]] == as.name("(")) cm <- cm[[2]]
pms <- list()
policy <- NULL
if(there.is.pm <- (cm[[1]] == as.name("&"))) { # there is a PM part
pm <- cm[[3]]
cm <- cm[[2]]
# deal with PM part
if(pm[[1]] == as.name("(")) {
pm <- pm[[2]]
if(pm[[1]] != as.name("|")) {
## Case: No maintenance policy
#stop("Need a policy to manage Preventive Maintenance")
policy <- NULL
} else {
policy <- pm[[3]]
if(policy[[1]] == as.name("*")) {
## Case: Composition of maintenance policies
# recursive function to detect maintenance policies
run.over.policies<-function(p) {
if(p[[1]] == as.name("*")) {
run.over.policies(p[[2]])
run.over.policies(p[[3]])
} else if(is.name(p[[1]])) {
p[[1]] <- as.name(paste0(as.character(p[[1]]),".maintenance.policy"))
policies <<- c(policies,list(p))
}
}
## init policies and
policies <- list()
run.over.policies(policy)
## print(policies)
policy <- policies ##[[1]]
} else if(is.name(policy[[1]])) {
## Case: One maintenance policy
policy[[1]] <- as.name(paste0(as.character(policy[[1]]),".maintenance.policy"))
}
# PMs
pm <- pm[[2]]
}
# parser for pm
parse.pm <- function(pm) {
if(is.name(pm[[1]])) {
pm[[1]] <- as.name(paste0(as.character(pm[[1]]),".va.model"))
}
pm
}
cpt.pms <- 0
while(pm[[1]] == as.name("+") ) {
if(length(pm) == 3) {
pms[[cpt.pms <- cpt.pms + 1]] <- parse.pm(pm[[3]])
pm <- pm[[2]]
}
}
pms[[cpt.pms <- cpt.pms + 1]] <- parse.pm(pm)
} else stop("Need parenthesis around the Preventive Maintenance terms")
}
# deal with CM PART
cms <- list()
# parser for cm
parse.cm <- function(cm) {
# print(there.is.pm)
# print(cm)
if(there.is.pm) {
if(cm[[1]] == as.name("(")) cm <- cm[[2]]
else stop("CM needs a family!")
}
if(cm[[1]] != as.name("|")) stop("CM needs a family!")
family <- cm[[3]]
if(is.name(family[[1]])) {
family[[1]] <- as.name(paste0(as.character(family[[1]]),".family.cm"))
}
cm <- cm[[2]]
if(is.name(cm[[1]])) {
cm[[1]] <- as.name(paste0(as.character(cm[[1]]),".va.model"))
}
list(model=cm,family=family)
}
cpt.cms <- 0
while( cm[[1]] == as.name("+") ) {
if(length(cm) == 3) {
cms[[cpt.cms <- cpt.cms + 1]] <- parse.cm(cm[[3]])
cm <- cm[[2]]
}
}
cms[[cpt.cms <- cpt.cms + 1]] <- parse.cm(cm)
## Parse covariates
parse.covariates <- function(expr) {
form<-list()
params <- c()
add_term <- function(term,sign) {
##print(term)
if(term[[1]]==as.name("*")) {
form <<- c(as.character(term[[3]]),form)
if(sign == as.name("-")) {
if(!(is.numeric(eval(term[[2]])))) stop("Only + is admitted between covariates definition in Bayes case.") else param_expr <- eval.vam(parse(text=paste0(sign,as.character(eval(term[[2]])))))
} else param_expr <- as.vector(eval.vam(term[[2]]))
params<<- c(param_expr,params)
}
##print(list(form=form,params=params))
}
while(expr[[1]]==as.name("+") || expr[[1]]==as.name("-")) {
add_term(expr[[3]],as.character(expr[[1]]))
expr <- expr[[2]]
}
add_term(expr,"")
list(formula=eval(parse(text=paste0("~",paste(form,collapse="+")))),params=params)
}
convert.family <- function(fam) {
# eval.vam is here to evaluate the value if it is a symbol!
# We want to detect if there is covariates or not. Normally the operator | delimitating covariates has priority, expect possibly in Baysian case.
if(has.covariates <- (length(fam[[length(fam)]])>1 && fam[[length(fam)]][[1]] == as.name("|"))) {
covariates_expr <- fam[[length(fam)]][[3]]
fam[[length(fam)]] <- fam[[length(fam)]][[2]] # first argument of last terms becomes last argument of family
} else {
#In Bayesian case the ~ operator of the description of the prior distribution of the last argument of the family can possibly has priority on the operator | delimitating covariates
if(length(fam[[length(fam)]])>1 && fam[[length(fam)]][[1]] == as.name("~") && length(fam[[length(fam)]][[2]])>1 && fam[[length(fam)]][[2]][[1]] == as.name("|")){
has.covariates <- TRUE
covariates_expr <- fam[[length(fam)]][[2]][[3]]
fam[[length(fam)]][[2]] <- fam[[length(fam)]][[2]][[2]]
}
}
res<-list(
name=as.character(fam[[1]]),
params=sapply(fam[-1],function(e) as.vector(eval.vam(e)))
## instead of : params=sapply(cm$family[-1],as.vector)
## which does not work with negative real since element of tmp[-1] interpreted as call!
)
if(has.covariates) {
res$covariates <- parse.covariates(covariates_expr)
}
return(res)
}
convert.pm <- function(pm) {
n_pip<-c()
if(length(pm)>1){
for(i in 2:length(pm)){
if((length(pm[[i]])==3)&&(pm[[i]][[1]]==as.name("|"))) {
n_pip<-c(n_pip,i)
}
}
}
if(length(n_pip)==0) {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==1) numeric(0) else sapply(pm[2:length(pm)],function(e) as.vector(eval.vam(e))))
)
} else if(length(n_pip)==1) {
if(n_pip<(length(pm)-1)) {
stop("Maximum two arguments after a | in a maintenance effect!")
} else if(n_pip == length(pm)) {
if( typeof(tryCatch( as.double(eval.vam(pm[[length(pm)]][[3]])) ,error=function(e){FALSE},finally=function(e){TRUE}))!="logical"){
if((round(eval.vam(pm[[length(pm)]][[3]])) != eval.vam(pm[[length(pm)]][[3]]))||(round(eval.vam(pm[[length(pm)]][[3]]))<=0)) {
stop("Memory argument of a maintenance model has to be a strictly positive integer!")
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==2) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)]][[2]])))),
m=as.integer(eval.vam(pm[[length(pm)]][[3]]))
)
}
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==2) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)]][[2]])))),
extra=as.character(pm[[length(pm)]][[3]])
)
}
} else {
if( typeof(tryCatch( as.double(eval.vam(pm[[length(pm)-1]][[3]])) ,error=function(e){FALSE},finally=function(e){TRUE}))!="logical"){
if((round(eval.vam(pm[[length(pm)-1]][[3]]))!=eval.vam(pm[[length(pm)-1]][[3]]))||(round(eval.vam(pm[[length(pm)-1]][[3]]))<0)) {
stop("Memory argument of a maintenance model has to be a positive integer!")
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==3) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-2)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)-1]][[2]])))),
m=as.integer(eval.vam(pm[[length(pm)-1]][[3]])),
extra=as.character(pm[[length(pm)]])
)
}
} else {
if( typeof(tryCatch( as.double(eval.vam(pm[[length(pm)]])) ,error=function(e){FALSE},finally=function(e){TRUE}))=="logical"){
stop("At least one of the two argument of maintenance model after a | must be a memory that is to say a non negative positive integer!")
} else {
if((round(eval.vam(pm[[length(pm)]]))!=eval.vam(pm[[length(pm)]]))||(round(eval.vam(pm[[length(pm)]]))<0)) {
stop("Memory argument of a maintenance model has to be a positive integer!")
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==3) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-2)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)-1]][[2]])))),
m=as.integer(eval.vam(pm[[length(pm)]])),
extra=as.character(pm[[length(pm)-1]][[3]])
)
}
}
}
}
} else {
stop("Maximum one | in a maintenance effect!")
}
# if((length(pm)==1)||(pm[[length(pm)]][[1]]!=as.name("|"))) {
# list(
# name=as.character(pm[[1]]),
# params=as.vector(if(length(pm)==1) numeric(0) else sapply(pm[2:length(pm)],function(e) as.vector(eval(e))))
# )
# } else if ( typeof(tryCatch( as.double(eval(pm[[length(pm)]][[3]])) ,error=function(e){FALSE},finally=function(e){TRUE}))!="logical"){
# if((round(eval(pm[[length(pm)]][[3]]))!=eval(pm[[length(pm)]][[3]]))||(round(eval(pm[[length(pm)]][[3]]))<0)) {
# stop("Memory argument of a maintenance model has to be a positive integer!")
# } else {
# list(
# name=as.character(pm[[1]]),
# params=as.vector(if(length(pm)==2) pm[[2]][[2]] else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval(e))),as.vector(eval(pm[[length(pm)]][[2]])))),
# m=as.integer(eval(pm[[length(pm)]][[3]]))
# )
# }
# } else {
# list(
# name=as.character(pm[[1]]),
# params=as.vector(if(length(pm)==2) pm[[2]][[2]] else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval(e))),as.vector(eval(pm[[length(pm)]][[2]])))),
# extra=as.character(pm[[length(pm)]][[3]])
# )
# }
}
convert.mp <- function(mp) {#maintenance policy
if(is.null(mp)) list(name="None")
else if(is.list(mp)) {
list(name="MaintenancePolicyList",policies=lapply(mp,convert.mp))
}
else {
## The function defining the maintenance policy
## (registered in maintenance-policy-register.R or in any other R file)
mp.fct <- eval.vam(mp[[1]])
## params used in the call mp
pars <- as.list(match.call(mp.fct,mp))[-1]
## Default values are then automatically completed using declaration of maintenance policy
pars.default <- (as.list(mp.fct)->tmp)[-length(tmp)]
pars.default <- pars.default[sapply(pars.default,function(e) nchar(as.character(e)))!=0]
for(e in names(pars.default)) if(is.null(pars[[e]])) pars[[e]] <- pars.default[[e]]
##print(list(pars=pars))
## deal with model parameter which has a specific treatment
mod <- NULL
if(!is.null(pars[["model"]])) {
mod <- rcpp(eval.vam(pars[["model"]]))
pars[["model"]] <- NULL
}
res <- list(
name=as.character(mp[[1]]),
params=lapply(pars,eval.vam)
)
res[["with.model"]] <- !is.null(mod)
if(!is.null(mod)) res[["model"]] <- mod
res
}
}
res<-list(
response=response,
models=c(list(convert.pm(cms[[1]]$model)),lapply(pms[rev(seq(pms))],convert.pm)),
family=convert.family(cms[[1]]$family),
pm.policy=convert.mp(policy)
)
## covariates direct acces
res$covariates <- res$family$covariates
res$family$covariates <- NULL
res$max_memory <- max(1,unlist(sapply(res$models,function(e) e$m)),na.rm=TRUE)
res
}
# use substitute coef in vam formula
substitute.vam.formula <- function(formula,coef,model) {
if(missing(model)) model <- parse.vam.formula(formula)
if(missing(coef)) {
coef <- c(model$family$params,sapply(model$models,function(m) m$params))
if(!is.null(model$covariates)) coef <- c(coef,model$covariates$params)
}
if(!is.null(model$covariates)) {
nb_paramsCovariates <- length(model$covariates$params)
} else {
nb_paramsCovariates <- 0
}
coef<-unlist(coef)
nb_paramsFamily <- length(model$family$params)
nb_paramsCM <- length(model$models[[1]]$params)
nb_paramsPM <- sapply(model$models[-1],function(m) length(m$params))
form <- paste0(
paste(model$response,collapse=" & "),
"~ (",
strsplit(model$models[[1]]$name,"\\.")[[1]][1],
"(",
if(nb_paramsCM>0) paste(coef[nb_paramsFamily+(1:nb_paramsCM)],collapse=",") else "",
if(!is.null(model$models[[1]]$m) || !is.null(model$models[[1]]$extra)) {
extra <- c()
if(!is.null(model$models[[1]]$extra)) extra <- c(extra,model$models[[1]]$extra)
if(!is.null(model$models[[1]]$m)) extra <- c(extra,model$models[[1]]$m)
paste0("|",paste(extra,collapse=","))
} else "",
")",
"|",
strsplit(model$family$name,"\\.")[[1]][1],
"(",
paste(coef[1:nb_paramsFamily],collapse=","),
if(!is.null(model$covariates)) {
# unlist(nb_paramsPM) since nb_paramsPM is list() when no PM
tmp<-coef[nb_paramsFamily + nb_paramsCM + sum(unlist(nb_paramsPM)) + (1:nb_paramsCovariates)]
if(is.list(tmp)){##Bayesian case
paste0("| (",
paste(tmp,all.vars(model$covariates$formula),sep=") * ",collapse=" + (")
)
} else {
tmp[tmp<0] <- paste0("(",tmp[tmp<0],")")
paste0("|",
paste(tmp,all.vars(model$covariates$formula),sep=" * ",collapse=" + ")
)
}
} else "",
")",
")"
)
if(length(model$models)>1) {
pms <- model$models[-1]
form <- paste0(form,
" & (",
paste(
sapply(seq(pms),function(i) {
paste0(
strsplit(pms[[i]]$name,"\\.")[[1]][1],
"(",
if(nb_paramsPM[i]>0) paste(coef[nb_paramsFamily+nb_paramsCM+ifelse(i>1,sum(nb_paramsPM[1:(i-1)]),0)+(1:nb_paramsPM[i])],collapse=",") else "",
if(!is.null(pms[[i]]$m) || !is.null(pms[[i]]$extra)) {
extra <- c()
if(!is.null(pms[[i]]$extra)) extra <- c(extra,pms[[i]]$extra)
if(!is.null(pms[[i]]$m)) extra <- c(extra,pms[[i]]$m)
paste0("|",paste(extra,collapse=","))
} else "",
")"
)
}),
collapse=" + "
),
")"
)
}
form <- eval(parse(text=form),envir=globalenv())
form
}
priors.from.vam.formula <- function(model) {
flatten.params <- c(model$family$params,unlist(sapply(model$models,function(e) e$params)))
##if(!is.null(model$family$covariates)) flatten.params <- c(flatten.params,mode$family$covariates$params) ##Strange
if(!is.null(model$covariates)) flatten.params <- c(flatten.params,model$covariates$params)
if(all(sapply(flatten.params,class) == "formula") ) {
prior.families <- c("B","Beta","U","Unif","G","Gamma","Norm","N","NonInform","NInf","LNorm","LogNorm","LN")
## clear "|" expression
flatten.params <- lapply(flatten.params,function(e) if(!as.character(e[[2]][[1]]) %in% prior.families) e[[2]][[2]] else e[[2]])
## transform to list
parse.prior <- function(prior) {
## declare here all the priors
Beta <- B <- Be <- function(a,b) list(name="Beta.prior",params=c(a,b))
Gamma <- G <- function(a,s) list(name="Gamma.prior",params=c(a,s))
Unif <- U <- function(a=0,b=1) list(name="Unif.prior",params=c(a,b))
Norm <- N <- function(m=0,s=1) list(name="Norm.prior",params=c(m,s))
LNorm <- LogNorm <- LN <- function(m=0,s=1) list(name="LNorm.prior",params=c(m,s))
NonInform <- NInf <- NI <- function(init=1,init_sigma=1) list(name="NonInform.prior",params=c(init,init_sigma))
res <- eval(prior) ## TODO or NOT TODO: eval(prior,parent.frame())
class(res) <- res$name #to be accessible as a class in R
res
}
flatten.params <- lapply(flatten.params,parse.prior)
return(flatten.params)
} else {
warning("Not a formula for bayesian.vam object!")
NULL # means not ok!
}
}
### FOUND A BETTER WAY: see bayesian.vam
## substitute prior with mean prior to be used inside mle.vam as initialization of run.bayesian.vam
# bayesian.model.to.mle.model <- function(model,priors) {
# print(priors)
# if(!is.null(priors)) {
# k<-0
# for(i in 1:length(model$family$params)) {
# ##DEBUG: print("prior");print(k+1);print(priors[[k+1]]);print(update(priors[[k+1]]))
# model$family$params[[i]] <- update(priors[[k<-k+1]])
# }
# model$family$params <- unlist(model$family$params)
# for(j in 1:length(model$models)) {
# for(i in 1:length(model$models[[j]]$params)) {
# ##DEBUG: print("prior");print(k+1);print(priors[[k+1]]);print(update(priors[[k+1]]))
# model$models[[j]]$params[[i]] <- update(priors[[k<-k+1]])
# }
# model$models[[j]]$params <- unlist(model$models[[j]]$params)
# }
# return(model)
# }
# return(NULL)
# }
## RMK: these 2 init following functions were only used for mle but are useful now for Bayesian too!
init.par0 <- function(obj,par0) {
## save the initial param
if(is.null(obj$par0)) obj$par0 <- params(obj)
## parameters stuff!
if(missing(par0)) {
if("par" %in% names(obj)) {
param <- obj$par[-1]
alpha <- obj$par[1] #LD2
} else {#not the first run
param<-params(obj)[-1] #first run
alpha<-params(obj)[1] #LD2
}
} else if(is.null(par0)) {
param<-obj$par0[-1]
alpha<-obj$par0[1] #LD2
} else {
param<-par0[-1]
alpha<-par0[1] #LD2
}
list(param=param,alpha=alpha)
}
init.fixed.param <- function(param,fixed) {
## fixed and functions stuff!
if(missing(fixed)) {
fixed<-rep(FALSE,length(param))
alpha_fixed<-FALSE #LD2
} else if(is.numeric(fixed)) {
fixedInd<-fixed
fixed<-rep(FALSE,length(param))
fixed[fixedInd-1]<-TRUE #LD2: ajout du -1
alpha_fixed<-sum(fixedInd==1) #OR ( 1 %in% fixedInd which is more expressive)
} else {#LD2
alpha_fixed<-fixed[1] #LD2
fixed<-fixed[-1] #LD2
}#LD2
list(fixed=fixed,alpha_fixed=alpha_fixed)
}
rcqls/VAM documentation built on Jan. 14, 2024, 9:07 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.
# Simulation: sim.vam
sim.vam <- function(formula,data.covariates) {
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(sim.vam,formula=formula(formula),data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
## specify names of variables inside data.frame
if(is.null(model$response)) {
self$response.names <- c("Time","Type")
self$system.name <- "System"
} else {
if(length(model$response)==3) {
self$system.name <- model$response[1]
self$response.names <- tail(model$response,2)
} else if(length(model$response)==2) {
self$system.name <- "System"
self$response.names <- model$response
}
}
self$model<-model#debuggage LD
self$formula <- substitute.vam.formula(model=model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
rcpp <- new(SimVam,model)
rcpp
})
self
}
# TODO: when data provided, complete the data!
simulate.sim.vam <- function(sim, stop.policy, nb.system, cache.size=500,as.list=FALSE,data=NULL) {
# To have a first argument more readable
self <- sim
rcpp <- self$rcpp()
if(missing(stop.policy)) stop.policy <- if(is.null(sim$model.parsed$covariates)) 10 else 1
if(missing(nb.system)) nb.system <- if(is.null(sim$model.parsed$covariates)) 1 else nrow(sim$model.parsed$covariates$data)
self$stop.policy.last <- parse.stop.policy(deparse(substitute(stop.policy)))
stop.policy <- eval(self$stop.policy.last,parent.frame(1))
if(is.numeric(stop.policy)) {
if(stop.policy == as.integer(stop.policy)) {#integer
stop.policy <- EndAt(size=stop.policy)
} else stop.policy <- NULL
} else if(!inherits(stop.policy,"stop.policy")) {
stop.policy <- NULL
}
if(is.null(stop.policy)) {
warning("Argument stop.policy is not a proper one!")
return(invisible(NULL))
}
# default cache size
if(is.null(stop.policy$cache.size)) stop.policy$cache.size <- cache.size
# add stop.policy object
rcpp$add_stop_policy(stop.policy)
# verify the coherence of the data argument
if(!is.null(data)){
if(!all(sim$response.names%in%names(data))) {
warning("The column names of the data argument must contain the response name (if not precise in the R formula those are by default Time and Type)!")
return(invisible(NULL))
}
data<-data[data[,sim$response.names[2]]!=0,]
}
if(nb.system>1) {
# multisystem
if(as.list) df<-list()
for(i in 1:nb.system) {
rcpp$select_current_system(i-1) # i-1 because it is c++
dataTime<-c()
dataType<-c()
if(!is.null(data)){
if(sim$system.name%in%names(data)){
dataInit<-data[data[,sim$system.name]==i,sim$response.names]
dataTime<-dataInit[,sim$response.names[1]]
dataType<-dataInit[,sim$response.names[2]]
} else {
dataTime<-data[,sim$response.names[1]]
dataType<-data[,sim$response.names[2]]
}
}
df2 <- rcpp$simulate(stop.policy$cache.size,as.numeric(dataTime),as.numeric(dataType))[-1,]
names(df2) <- sim$response.names
if(as.list) {
rownames(df2)<-1:nrow(df2)
df[[i]] <- df2 #rbind(data.frame(Time=0,Type=1),df2)
} else {
df2[[sim$system.name]] <- i
df2<-df2[c(3,1:2)]
df <- if(i==1) df2 else rbind(df,df2)
}
}
} else {
rcpp$select_current_system(0) # added for covariates case
df <- rcpp$simulate(stop.policy$cache.size,as.numeric(data[,sim$response.names[1]]),as.numeric(data[,sim$response.names[2]]))[-1,]
names(df) <- sim$response.names
if(as.list) {rownames(df)<-1:nrow(df);df<-list(df)}
}
if(!as.list) {
rownames(df) <- 1:nrow(df)
rcpp$set_data(data.frame.to.list.multi.vam(df,names(df))) #if already list, response not used inside data.frame.to.list.multi.vam
}
else {
names(df) <- paste0(sim$system.name,1:length(df))
rcpp$set_data(unname(df))
}
## put the final result transformed as in mle.vam and model.vam to the model
## return the result as a data.frame
df
}
# Model part
model.vam <- function(formula,data,data.covariates) {
if(missing(data)) data<-NULL
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(model.vam,formula=formula(formula),data=data,data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
self$formula <- substitute.vam.formula(model=model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
if(is.null(self$data)) {## No data
rcpp <- new(ModelVam,model)
rcpp
} else {
response <- model$response
data <- data.frame.to.list.multi.vam(self$data,response)
rcpp <- new(ModelVam,model,data)
rcpp
}
})
self
}
update.model.vam <- function(model,data) {
if(!missing(data)) {
self <- model #to have an argument more readable
response <- parse.vam.formula(self$formula)$response
self$data <- data
data2 <- data.frame.to.list.multi.vam(self$data,response)
self$rcpp()$set_data(data2)
}
}
check.censorship <- function(data) {
lapply(data,function(data.syst) {
ty <- data.syst$Type[-1]
check<-list(type="none",index=0)
if(any(ty==0)) {
if(which(ty==0)[1]==length(ty)) check$type <- 'right' #meaning: right only
else if(!any(ty<0) || (which(ty<0)[1] > which(ty==0)[1])) {
check$type <- 'left' #can have also a right censorship
check$index <- which(ty==0)[1]
}
else check$type <- 'unknown'
}
check
})
}
data.frame.to.list.multi.vam <- function(data,response) {
# return data if it is already only a list!
if(is.list(data) && !is.data.frame(data)) {
data.frame(Time=0,Type=1)->tmp
names(tmp) <- names(data[[1]])
return(lapply(data,function(df) rbind(tmp,df)))
}
if(!("System" %in% response) && ("System" %in% names(data)) ) {
warning(paste0("WARNING: data has variable 'System' when response in formula does not contain this variable!"))
}
# otherwise
if(length(response)==2) {
if(length(intersect(response,names(data))) != 2) stop(paste0("Bad response:",response))
tmp <- data[[response[1]]]
data2 <- list(data.frame(Time=c(0,tmp[order(tmp)]),Type=c(1,data[[response[2]]][order(tmp)])))
} else {
if(length(intersect(response,names(data))) != 3) stop(paste0("Bad response:",response))
syst0 <- unique(syst<-data[[response[1]]])
data2 <- list()
for(i in seq_along(syst0)) {
df <- data[syst==syst0[i],response]
tmp <- df[[response[2]]]
data2[[i]] <- data.frame(Time=c(0,tmp[order(tmp)]),Type=c(1,df[[response[3]]][order(tmp)]))
}
}
data2
}
# TODO: check data
check.data.vam <-function(data,response) {
if(all(data[[response[[-length(response)]]]])) {
}
}
make.censorship <- function(data,rcpp) {
if(any(sapply(censorship <- check.censorship(data),function(e) e$type == 'unknown'))) stop("VAM package does know how to deal with this kind of censorship!")
if(any(sapply(censorship,function(e) e$type == 'left'))) {
leftCensors <- sapply(censorship,function(e) e$index)
rcpp$set_leftCensors(leftCensors)
}
}
mle.vam <- function(formula,data,data.covariates) {
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(mle.vam,formula=formula(formula),data=data,data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
self$formula <- substitute.vam.formula(model=model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
response <- model$response
data <- data.frame.to.list.multi.vam(self$data,response)
rcpp <- new(MLEVam,model,data)
make.censorship(data,rcpp)
rcpp
})
self
}
# to convert in Rcpp
params.model.vam <- params.sim.vam <- params.mle.vam <- params.bayesian.vam <- function(self,param) {
if(missing(param)) {
self$rcpp()$get_params()
} else {
self$rcpp()$set_params(param)
}
}
## Useless since stats:::formula.default do that by default
# formula.model.vam <- formula.sim.vam <- function(self) {
# self$formula
# }
formula.mle.vam <- function(self,origin=FALSE) {
form <- substitute.vam.formula(self$formula,coef(self))
if(origin) list(formula=form,origin=self$formula)
else form
}
update.mle.vam <- function(mle,data) {
if(!missing(data)) {
self <- mle
model <- parse.vam.formula(self$formula)
response <- model$response
self$data <- data
data2 <- data.frame.to.list.multi.vam(self$data,response)
self$rcpp()$set_data(data2)
self$rcpp()$reset_leftCensors()
make.censorship(data2,self$rcpp())
## estimation has to be computed again!
self$mle.coef<-NULL
}
## with
}
#fonction de LD2
contrast.mle.vam <-function(obj,par0,with_value=TRUE,with_gradient=FALSE,with_hessian=FALSE){
type <- c(with_value,with_gradient,with_hessian)
rcpp <- obj$rcpp()
## save the initial param
if(is.null(obj$par0)) obj$par0 <- params(obj)
## parameters stuff!
if(missing(par0)) {
if("par" %in% names(obj)) {
param <- obj$par[-1]
alpha <- obj$par[1]#LD2
} #not the first run
else {
param<-params(obj)[-1] #first run
alpha<-params(obj)[1]#LD2
}
} else if(is.null(par0)) {
param<-obj$par0[-1]
alpha<-obj$par0[1]#LD2
} else {
param<-par0[-1]
alpha<-par0[1]#LD2
}
if(sum(type)==1) {
if(type[1]){
res<-rcpp$contrast(c(alpha,param),FALSE)
} else if(type[2]){
res<-rcpp$gradient(c(alpha,param),FALSE)[-1]
} else if(type[3]){
res<-rcpp$hessian(c(alpha,param),FALSE)[-1,-1]
}
} else {
res<-list()
if (type[1]){
res$contrast<-rcpp$contrast(c(alpha,param),FALSE)
}
if(type[2]){
res$gradient<-rcpp$gradient(c(alpha,param),FALSE)[-1]
}
if(type[3]){
res$hessian<-rcpp$hessian(c(alpha,param),FALSE)[-1,-1]
}
}
res
}
#fonction de LD2
logLik.mle.vam <-function(obj,par0,with_value=TRUE,with_gradient=FALSE,with_hessian=FALSE){
type=c(with_value,with_gradient,with_hessian)
rcpp <- obj$rcpp()
## save the initial param
if(is.null(obj$par0)) obj$par0 <- params(obj)
## parameters stuff!
if(missing(par0)) {
if("par" %in% names(obj)) {
param <- obj$par[-1]
alpha <- obj$par[1]#LD2
} #not the first run
else {
param<-params(obj)[-1] #first run
alpha<-params(obj)[1]#LD2
}
} else if(is.null(par0)) {
param<-obj$par0[-1]
alpha<-obj$par0[1]#LD2
} else {
param<-par0[-1]
alpha<-par0[1]#LD2
}
if(sum(type)==1) {
if(type[1]){
res<-rcpp$contrast(c(alpha,param),TRUE)
} else if(type[2]){
res<-rcpp$gradient(c(alpha,param),TRUE)
} else if(type[3]){
res<-rcpp$hessian(c(alpha,param),TRUE)
}
} else {
res<-list()
if (type[1]){
res$contrast<-rcpp$contrast(c(alpha,param),TRUE)
}
if(type[2]){
res$gradient<-rcpp$gradient(c(alpha,param),TRUE)
}
if(type[3]){
res$hessian<-rcpp$hessian(c(alpha,param),TRUE)
}
}
res
}
# alpha is not considered in the estimation!
run.mle.vam <-function(obj,par0,fixed,method=NULL,verbose=FALSE,...) {
rcpp <- obj$rcpp()
par0.tmp <- init.par0(obj,par0)
param <- par0.tmp$param
alpha <- par0.tmp$alpha
# ## save the initial param
# if(is.null(obj$par0)) obj$par0 <- params(obj)
# ## parameters stuff!
# if(missing(par0)) {
# if("par" %in% names(obj)) {
# param <- obj$par[-1]
# alpha <- obj$par[1] #LD2
# } else {#not the first run
# param<-params(obj)[-1] #first run
# alpha<-params(obj)[1] #LD2
# }
# } else if(is.null(par0)) {
# param<-obj$par0[-1]
# alpha<-obj$par0[1] #LD2
# } else {
# param<-par0[-1]
# alpha<-par0[1] #LD2
# }
## fixed and functions stuff!
fixed.tmp <- init.fixed.param(param,fixed)
fixed <- fixed.tmp$fixed
alpha_fixed <- fixed.tmp$alpha_fixed
fn<-function(par) {
##cat("param->");print(par);print(param[!fixed])
param[!fixed]<-par
#cat("param->");print(param)
## All the commented part allows us to save the param when value is NaN
#res<-
-rcpp$contrast(c(alpha,param),alpha_fixed) #LD2
# if(is.nan(res)) {
# mode_param<-"contrast"
# data_param<- rcpp$get_data(0)
# data_param2 <- obj$data
# save(param,mode_param,data_param,data_param2,file="/Users/remy/tmp/VAM/res.RData")
# }
# res
}
gr <- function(par) {
param[!fixed]<-par
#cat("param2->");print(param)
## All the commented part allows us to save the param when value is NaN
#res <-
(-rcpp$gradient(c(alpha,param),alpha_fixed)[-1])[!fixed] #LD2
# if(any(is.nan(res))) {
# mode_param<-"gradient"
# save(param,mode_param,file="/Users/remy/tmp/VAM/res.RData")
# }
# res
}
## optim stuff!
if(is.null(method) || method=="fast") {
if(length(param[!fixed])>1) param[!fixed]<-(res <- optim(param[!fixed],fn,gr,method="Ne",...))$par
if(is.null(method)) res<-optim(param[!fixed],fn,gr,method="CG",...)
} else {
res<-optim(param[!fixed],fn,gr,method=method,...)
}
#fixed tips
param[!fixed]<-res$par
if(!alpha_fixed){#LD2
## complete the scale parameter
alpha <- rcpp$alpha_est(c(1,param)) #LD2
} #LD2
res$par<-c(alpha,param) #LD2
if(verbose) print(res)
## save stuff
obj$fixed <- c(alpha_fixed,fixed) #LD2
obj$optim<-res
obj$par<-res$par
obj$mle.coef <- res$par
params(obj,obj$mle.coef) #put the result in the c++ part
##obj$mle.coef
obj$optim
}
## Rmk: run.mle.vam is supposed to run many times to get the best estimate!
## Here, par=NULL forces initialisation update but does not ensure that it is the best estimate.
## TODO: try to find a best strategy or many strategies...
## ... added to deal with fixed by example!
coef.mle.vam <- function(obj,par=NULL,method=NULL,verbose=FALSE,...) {
if(is.null(obj$mle.coef) || !is.null(par)) {
res <-run.mle.vam(obj,par,verbose=verbose,method=method,...)
if(verbose && obj$optim$convergence>0) cat("convergence=",obj$optim$convergence,"\n",sep="")
}
obj$mle.coef
}
# bayesian.vam class
bayesian.vam <- function(formula,data,data.covariates) {
if(missing(data.covariates)) data.covariates<-NULL
self <- newEnv(bayesian.vam,formula=formula,data=data,data.covariates=data.covariates)
PersistentRcppObject(self,new = {
model <- parse.vam.formula(self$formula)
self$formula <- substitute.vam.formula(model=model)
response <- model$response
data <- data.frame.to.list.multi.vam(self$data,response)
self$priors <- priors.from.vam.formula(model)
##DEBUG: print("priors");print(priors)
##DEBUG: print("modelAV");print(model)
self$prior.params <- sapply(self$priors,mean)
self$mle.formula <- substitute.vam.formula(self$formula,self$prior.params)
## THIS IS LESS CLEVER THAN THE NEXT LINE: print(model<-bayesian.model.to.mle.model(model,priors))
model<-parse.vam.formula(self$mle.formula)
##DEBUG: print("modelAP");print(model)
if(!is.null(model$covariates)) {
model$covariates$data <- model.frame(model$covariates$formula,data=self$data.covariates) #ok even if data.covariates is null
}
self$model.parsed <- model #just to keep some info notably on covariates
rcpp <- new(BayesianVam,model,data,self$priors)
rcpp
})
self
}
run.bayesian.vam <- function(obj,par0,fixed,sigma.proposal,nb=100000,burn=10000,profile.alpha=FALSE,method=NULL,verbose=FALSE,history=FALSE,proposal='norm',...) {
rcpp <- obj$rcpp()
obj$history<-history
obj$nb<-nb
obj$burn<-burn
obj$preplots <- NULL
## init via mle: par0 is supposed first to be initialized by mle
if(missing(par0)) {
if(is.null(obj$model.parsed$covariates)) obj$mle <- mle.vam(obj$mle.formula,obj$data) else obj$mle <- mle.vam(obj$mle.formula,obj$data,obj$data.covariates)
obj$mle.init <- TRUE
obj$par0 <- coef(obj$mle,fixed=fixed,method=method,verbose=verbose,...)
} else {
obj$par0<-par0
obj$mle.init <- FALSE
}
fixed.tmp <- init.fixed.param(obj$par0[-1],fixed)
obj$fixed <- c(fixed.tmp$alpha_fixed,fixed.tmp$fixed)
#obj$alpha_fixed <- fixed.tmp$alpha_fixed
obj$profile_alpha<-profile.alpha
if (fixed.tmp$alpha_fixed & profile.alpha){
warning("Parameter alpha can not simultaneously be fixed and optimized: it is fixed !")
obj$profile_alpha<-FALSE
}
##print(obj$mle.init)
if(missing(sigma.proposal)) sigma.proposal <- sapply(obj$priors,sigma)
else {
if(length(sigma.proposal)==1) sigma.proposal <- rep(sigma.proposal,length(obj$priors))
}
if(length(proposal)==1) proposal <- rep(proposal,length(obj$priors))
for(i in (1:length(obj$priors))) {
rcpp$set_sigma(i-1,sigma.proposal[i])
rcpp$set_proposal(i-1,switch(proposal[i],'lnorm'=1,0))
}
obj$sigma_proposal<-sigma.proposal
obj$proposal<-proposal
if(history) {
res <- rcpp$mcmc_history(obj$par0,nb,burn,obj$fixed,obj$profile_alpha)
obj$nb<-res[[3+obj$profile_alpha]]
obj$nb_proposal<-nb
obj$par<-as.data.frame(res[1:(2+obj$profile_alpha)])
names(obj$par) <- c("ind","estimate","alpha")[1:(2+obj$profile_alpha)]
} else {
obj$par <- rcpp$mcmc(obj$par0,nb,burn,obj$fixed,obj$profile_alpha)
}
obj$par
}
coef.bayesian.vam <- function(obj,new.run=FALSE,...) {
if(new.run || is.null(obj$par)) run(obj,...)
if(obj$history){
param<-sapply(obj$profile_alpha:(length(obj$par0)-1),function(j){mean(obj$par$estimate[obj$par$ind==j])})
if(obj$profile_alpha){
## complete the scale parameter
param <- c(mean(obj$par$alpha),param)
}
} else {
param <- sapply(obj$par,mean)
}
param[as.logical(obj$fixed)]<-obj$par0[as.logical(obj$fixed)]
param
}
hist.bayesian.vam <- function(obj,i=1,...) {
if(is.null(obj$par)) run(obj)
if(obj$fixed[i]){
warning("No hist for a fixed parameter!")
} else {
if(obj$history){
if((obj$profile_alpha)&(i==1)){
thetak<-obj$par$alpha
} else {
thetak<-obj$par$estimate[obj$par$ind==(i-1)]
}
hist(thetak,prob=TRUE,...)
abline(v=mean(thetak),col="blue",lwd=2)
} else{
hist(obj$par[[i]],prob=TRUE,...)
abline(v=mean(obj$par[[i]]),col="blue",lwd=2)
}
}
}
summary.bayesian.vam <- function(obj,alpha=0.05,new.run=FALSE,digits=4,...) {
if(new.run || is.null(obj$par)) run(obj,...)
res1<-obj$par0
cat("Initial parameters",if(obj$mle.init) " (by MLE)" else "",": ",paste(signif(res1,digits=digits),collapse=", "),"\n",sep="")
res2<-coef(obj)
cat("(Mean) Bayesian estimates: ", paste(signif(res2,digits=digits),collapse=", "),"\n",sep="")
if(obj$history){
thetak<-lapply(obj$profile_alpha:(length(obj$par0)-1),function(j){obj$par$estimate[obj$par$ind==j]})
if(obj$profile_alpha){thetak<-c(list(obj$par$alpha),thetak)}
} else {
thetak<-obj$par
}
sd<-sapply(thetak,sd)
sd[as.logical(obj$fixed)]<-0
cat("(SD) Bayesian estimates: ",paste(signif(sd,digits=digits),collapse=", "),"\n",sep="")
res4<-sapply(thetak,function(x){quantile(x,probs=alpha/2)})
cat("(", alpha/2,"-Quantile) Bayesian estimates: ",paste(signif(res4,digits=digits),collapse=", "),"\n",sep="")
res5<-sapply(thetak,function(x){quantile(x,probs=1-alpha/2)})
cat("(", 1-alpha/2,"-Quantile) Bayesian estimates: ",paste(signif(res5,digits=digits),collapse=", "),"\n",sep="")
res6<-sapply(thetak,length)
cat("(Number) Bayesian estimates: ",paste(res6,collapse=", "),"\n",sep="")
res7<-sapply(thetak,length)/(obj$nb-obj$burn)
cat("Metropolis-Hasting acceptation rates: ",paste(signif(res7,digits=digits),collapse=", "),"\n",sep="")
res<-data.frame(res1,res2,sd,res4,res5,res6,res7)
names(res)<-c(paste("Init",if(obj$mle.init) "(MLE)" else "",sep=""),"Mean","SD",paste(alpha/2,"-Quantile",sep=""),paste(1-alpha/2,"-Quantile",sep=""),"Number","Accept_Rate")
invisible(res)
}
# for both sim and mle
parse.vam.formula <- function(formula) {
## Needs to have this envir to evaluate params (otherwise, beta was found in baseenv() first before globalenv() for example!)
envir.eval <- parent.frame(4)
eval.vam <- function(e) eval(e,envir.eval)
if(formula[[1]] != as.name("~")) stop("Argument has to be a formula")
if(length(formula) == 2) {
response <- NULL
cm <- formula[[2]]
} else {
tmp <- formula[[2]]
## simplify parenthesis
while(tmp[[1]] == as.name("(")) tmp <- tmp[[2]]
if(tmp[[1]] != as.name("&") && length(tmp) != 3) stop("Left part of formula of the form 'Time & Type'!")
if(length(tmp[[2]])==3 && tmp[[2]][[1]]==as.name("&")) {
response <- c(as.character(tmp[[2]][[2]]),as.character(tmp[[2]][[3]]),as.character(tmp[[3]]))
} else response <- c(as.character(tmp[[2]]),as.character(tmp[[3]]))
cm <- formula[[3]]
}
## simplify parenthesis
while(cm[[1]] == as.name("(")) cm <- cm[[2]]
pms <- list()
policy <- NULL
if(there.is.pm <- (cm[[1]] == as.name("&"))) { # there is a PM part
pm <- cm[[3]]
cm <- cm[[2]]
# deal with PM part
if(pm[[1]] == as.name("(")) {
pm <- pm[[2]]
if(pm[[1]] != as.name("|")) {
## Case: No maintenance policy
#stop("Need a policy to manage Preventive Maintenance")
policy <- NULL
} else {
policy <- pm[[3]]
if(policy[[1]] == as.name("*")) {
## Case: Composition of maintenance policies
# recursive function to detect maintenance policies
run.over.policies<-function(p) {
if(p[[1]] == as.name("*")) {
run.over.policies(p[[2]])
run.over.policies(p[[3]])
} else if(is.name(p[[1]])) {
p[[1]] <- as.name(paste0(as.character(p[[1]]),".maintenance.policy"))
policies <<- c(policies,list(p))
}
}
## init policies and
policies <- list()
run.over.policies(policy)
## print(policies)
policy <- policies ##[[1]]
} else if(is.name(policy[[1]])) {
## Case: One maintenance policy
policy[[1]] <- as.name(paste0(as.character(policy[[1]]),".maintenance.policy"))
}
# PMs
pm <- pm[[2]]
}
# parser for pm
parse.pm <- function(pm) {
if(is.name(pm[[1]])) {
pm[[1]] <- as.name(paste0(as.character(pm[[1]]),".va.model"))
}
pm
}
cpt.pms <- 0
while(pm[[1]] == as.name("+") ) {
if(length(pm) == 3) {
pms[[cpt.pms <- cpt.pms + 1]] <- parse.pm(pm[[3]])
pm <- pm[[2]]
}
}
pms[[cpt.pms <- cpt.pms + 1]] <- parse.pm(pm)
} else stop("Need parenthesis around the Preventive Maintenance terms")
}
# deal with CM PART
cms <- list()
# parser for cm
parse.cm <- function(cm) {
# print(there.is.pm)
# print(cm)
if(there.is.pm) {
if(cm[[1]] == as.name("(")) cm <- cm[[2]]
else stop("CM needs a family!")
}
if(cm[[1]] != as.name("|")) stop("CM needs a family!")
family <- cm[[3]]
if(is.name(family[[1]])) {
family[[1]] <- as.name(paste0(as.character(family[[1]]),".family.cm"))
}
cm <- cm[[2]]
if(is.name(cm[[1]])) {
cm[[1]] <- as.name(paste0(as.character(cm[[1]]),".va.model"))
}
list(model=cm,family=family)
}
cpt.cms <- 0
while( cm[[1]] == as.name("+") ) {
if(length(cm) == 3) {
cms[[cpt.cms <- cpt.cms + 1]] <- parse.cm(cm[[3]])
cm <- cm[[2]]
}
}
cms[[cpt.cms <- cpt.cms + 1]] <- parse.cm(cm)
## Parse covariates
parse.covariates <- function(expr) {
form<-list()
params <- c()
add_term <- function(term,sign) {
##print(term)
if(term[[1]]==as.name("*")) {
form <<- c(as.character(term[[3]]),form)
if(sign == as.name("-")) {
if(!(is.numeric(eval(term[[2]])))) stop("Only + is admitted between covariates definition in Bayes case.") else param_expr <- eval.vam(parse(text=paste0(sign,as.character(eval(term[[2]])))))
} else param_expr <- as.vector(eval.vam(term[[2]]))
params<<- c(param_expr,params)
}
##print(list(form=form,params=params))
}
while(expr[[1]]==as.name("+") || expr[[1]]==as.name("-")) {
add_term(expr[[3]],as.character(expr[[1]]))
expr <- expr[[2]]
}
add_term(expr,"")
list(formula=eval(parse(text=paste0("~",paste(form,collapse="+")))),params=params)
}
convert.family <- function(fam) {
# eval.vam is here to evaluate the value if it is a symbol!
# We want to detect if there is covariates or not. Normally the operator | delimitating covariates has priority, expect possibly in Baysian case.
if(has.covariates <- (length(fam[[length(fam)]])>1 && fam[[length(fam)]][[1]] == as.name("|"))) {
covariates_expr <- fam[[length(fam)]][[3]]
fam[[length(fam)]] <- fam[[length(fam)]][[2]] # first argument of last terms becomes last argument of family
} else {
#In Bayesian case the ~ operator of the description of the prior distribution of the last argument of the family can possibly has priority on the operator | delimitating covariates
if(length(fam[[length(fam)]])>1 && fam[[length(fam)]][[1]] == as.name("~") && length(fam[[length(fam)]][[2]])>1 && fam[[length(fam)]][[2]][[1]] == as.name("|")){
has.covariates <- TRUE
covariates_expr <- fam[[length(fam)]][[2]][[3]]
fam[[length(fam)]][[2]] <- fam[[length(fam)]][[2]][[2]]
}
}
res<-list(
name=as.character(fam[[1]]),
params=sapply(fam[-1],function(e) as.vector(eval.vam(e)))
## instead of : params=sapply(cm$family[-1],as.vector)
## which does not work with negative real since element of tmp[-1] interpreted as call!
)
if(has.covariates) {
res$covariates <- parse.covariates(covariates_expr)
}
return(res)
}
convert.pm <- function(pm) {
n_pip<-c()
if(length(pm)>1){
for(i in 2:length(pm)){
if((length(pm[[i]])==3)&&(pm[[i]][[1]]==as.name("|"))) {
n_pip<-c(n_pip,i)
}
}
}
if(length(n_pip)==0) {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==1) numeric(0) else sapply(pm[2:length(pm)],function(e) as.vector(eval.vam(e))))
)
} else if(length(n_pip)==1) {
if(n_pip<(length(pm)-1)) {
stop("Maximum two arguments after a | in a maintenance effect!")
} else if(n_pip == length(pm)) {
if( typeof(tryCatch( as.double(eval.vam(pm[[length(pm)]][[3]])) ,error=function(e){FALSE},finally=function(e){TRUE}))!="logical"){
if((round(eval.vam(pm[[length(pm)]][[3]])) != eval.vam(pm[[length(pm)]][[3]]))||(round(eval.vam(pm[[length(pm)]][[3]]))<=0)) {
stop("Memory argument of a maintenance model has to be a strictly positive integer!")
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==2) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)]][[2]])))),
m=as.integer(eval.vam(pm[[length(pm)]][[3]]))
)
}
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==2) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)]][[2]])))),
extra=as.character(pm[[length(pm)]][[3]])
)
}
} else {
if( typeof(tryCatch( as.double(eval.vam(pm[[length(pm)-1]][[3]])) ,error=function(e){FALSE},finally=function(e){TRUE}))!="logical"){
if((round(eval.vam(pm[[length(pm)-1]][[3]]))!=eval.vam(pm[[length(pm)-1]][[3]]))||(round(eval.vam(pm[[length(pm)-1]][[3]]))<0)) {
stop("Memory argument of a maintenance model has to be a positive integer!")
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==3) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-2)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)-1]][[2]])))),
m=as.integer(eval.vam(pm[[length(pm)-1]][[3]])),
extra=as.character(pm[[length(pm)]])
)
}
} else {
if( typeof(tryCatch( as.double(eval.vam(pm[[length(pm)]])) ,error=function(e){FALSE},finally=function(e){TRUE}))=="logical"){
stop("At least one of the two argument of maintenance model after a | must be a memory that is to say a non negative positive integer!")
} else {
if((round(eval.vam(pm[[length(pm)]]))!=eval.vam(pm[[length(pm)]]))||(round(eval.vam(pm[[length(pm)]]))<0)) {
stop("Memory argument of a maintenance model has to be a positive integer!")
} else {
list(
name=as.character(pm[[1]]),
params=as.vector(if(length(pm)==3) as.vector(eval.vam(pm[[2]][[2]])) else c(sapply(pm[2:(length(pm)-2)],function(e) as.vector(eval.vam(e))),as.vector(eval.vam(pm[[length(pm)-1]][[2]])))),
m=as.integer(eval.vam(pm[[length(pm)]])),
extra=as.character(pm[[length(pm)-1]][[3]])
)
}
}
}
}
} else {
stop("Maximum one | in a maintenance effect!")
}
# if((length(pm)==1)||(pm[[length(pm)]][[1]]!=as.name("|"))) {
# list(
# name=as.character(pm[[1]]),
# params=as.vector(if(length(pm)==1) numeric(0) else sapply(pm[2:length(pm)],function(e) as.vector(eval(e))))
# )
# } else if ( typeof(tryCatch( as.double(eval(pm[[length(pm)]][[3]])) ,error=function(e){FALSE},finally=function(e){TRUE}))!="logical"){
# if((round(eval(pm[[length(pm)]][[3]]))!=eval(pm[[length(pm)]][[3]]))||(round(eval(pm[[length(pm)]][[3]]))<0)) {
# stop("Memory argument of a maintenance model has to be a positive integer!")
# } else {
# list(
# name=as.character(pm[[1]]),
# params=as.vector(if(length(pm)==2) pm[[2]][[2]] else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval(e))),as.vector(eval(pm[[length(pm)]][[2]])))),
# m=as.integer(eval(pm[[length(pm)]][[3]]))
# )
# }
# } else {
# list(
# name=as.character(pm[[1]]),
# params=as.vector(if(length(pm)==2) pm[[2]][[2]] else c(sapply(pm[2:(length(pm)-1)],function(e) as.vector(eval(e))),as.vector(eval(pm[[length(pm)]][[2]])))),
# extra=as.character(pm[[length(pm)]][[3]])
# )
# }
}
convert.mp <- function(mp) {#maintenance policy
if(is.null(mp)) list(name="None")
else if(is.list(mp)) {
list(name="MaintenancePolicyList",policies=lapply(mp,convert.mp))
}
else {
## The function defining the maintenance policy
## (registered in maintenance-policy-register.R or in any other R file)
mp.fct <- eval.vam(mp[[1]])
## params used in the call mp
pars <- as.list(match.call(mp.fct,mp))[-1]
## Default values are then automatically completed using declaration of maintenance policy
pars.default <- (as.list(mp.fct)->tmp)[-length(tmp)]
pars.default <- pars.default[sapply(pars.default,function(e) nchar(as.character(e)))!=0]
for(e in names(pars.default)) if(is.null(pars[[e]])) pars[[e]] <- pars.default[[e]]
##print(list(pars=pars))
## deal with model parameter which has a specific treatment
mod <- NULL
if(!is.null(pars[["model"]])) {
mod <- rcpp(eval.vam(pars[["model"]]))
pars[["model"]] <- NULL
}
res <- list(
name=as.character(mp[[1]]),
params=lapply(pars,eval.vam)
)
res[["with.model"]] <- !is.null(mod)
if(!is.null(mod)) res[["model"]] <- mod
res
}
}
res<-list(
response=response,
models=c(list(convert.pm(cms[[1]]$model)),lapply(pms[rev(seq(pms))],convert.pm)),
family=convert.family(cms[[1]]$family),
pm.policy=convert.mp(policy)
)
## covariates direct acces
res$covariates <- res$family$covariates
res$family$covariates <- NULL
res$max_memory <- max(1,unlist(sapply(res$models,function(e) e$m)),na.rm=TRUE)
res
}
# use substitute coef in vam formula
substitute.vam.formula <- function(formula,coef,model) {
if(missing(model)) model <- parse.vam.formula(formula)
if(missing(coef)) {
coef <- c(model$family$params,sapply(model$models,function(m) m$params))
if(!is.null(model$covariates)) coef <- c(coef,model$covariates$params)
}
if(!is.null(model$covariates)) {
nb_paramsCovariates <- length(model$covariates$params)
} else {
nb_paramsCovariates <- 0
}
coef<-unlist(coef)
nb_paramsFamily <- length(model$family$params)
nb_paramsCM <- length(model$models[[1]]$params)
nb_paramsPM <- sapply(model$models[-1],function(m) length(m$params))
form <- paste0(
paste(model$response,collapse=" & "),
"~ (",
strsplit(model$models[[1]]$name,"\\.")[[1]][1],
"(",
if(nb_paramsCM>0) paste(coef[nb_paramsFamily+(1:nb_paramsCM)],collapse=",") else "",
if(!is.null(model$models[[1]]$m) || !is.null(model$models[[1]]$extra)) {
extra <- c()
if(!is.null(model$models[[1]]$extra)) extra <- c(extra,model$models[[1]]$extra)
if(!is.null(model$models[[1]]$m)) extra <- c(extra,model$models[[1]]$m)
paste0("|",paste(extra,collapse=","))
} else "",
")",
"|",
strsplit(model$family$name,"\\.")[[1]][1],
"(",
paste(coef[1:nb_paramsFamily],collapse=","),
if(!is.null(model$covariates)) {
# unlist(nb_paramsPM) since nb_paramsPM is list() when no PM
tmp<-coef[nb_paramsFamily + nb_paramsCM + sum(unlist(nb_paramsPM)) + (1:nb_paramsCovariates)]
if(is.list(tmp)){##Bayesian case
paste0("| (",
paste(tmp,all.vars(model$covariates$formula),sep=") * ",collapse=" + (")
)
} else {
tmp[tmp<0] <- paste0("(",tmp[tmp<0],")")
paste0("|",
paste(tmp,all.vars(model$covariates$formula),sep=" * ",collapse=" + ")
)
}
} else "",
")",
")"
)
if(length(model$models)>1) {
pms <- model$models[-1]
form <- paste0(form,
" & (",
paste(
sapply(seq(pms),function(i) {
paste0(
strsplit(pms[[i]]$name,"\\.")[[1]][1],
"(",
if(nb_paramsPM[i]>0) paste(coef[nb_paramsFamily+nb_paramsCM+ifelse(i>1,sum(nb_paramsPM[1:(i-1)]),0)+(1:nb_paramsPM[i])],collapse=",") else "",
if(!is.null(pms[[i]]$m) || !is.null(pms[[i]]$extra)) {
extra <- c()
if(!is.null(pms[[i]]$extra)) extra <- c(extra,pms[[i]]$extra)
if(!is.null(pms[[i]]$m)) extra <- c(extra,pms[[i]]$m)
paste0("|",paste(extra,collapse=","))
} else "",
")"
)
}),
collapse=" + "
),
")"
)
}
form <- eval(parse(text=form),envir=globalenv())
form
}
priors.from.vam.formula <- function(model) {
flatten.params <- c(model$family$params,unlist(sapply(model$models,function(e) e$params)))
##if(!is.null(model$family$covariates)) flatten.params <- c(flatten.params,mode$family$covariates$params) ##Strange
if(!is.null(model$covariates)) flatten.params <- c(flatten.params,model$covariates$params)
if(all(sapply(flatten.params,class) == "formula") ) {
prior.families <- c("B","Beta","U","Unif","G","Gamma","Norm","N","NonInform","NInf","LNorm","LogNorm","LN")
## clear "|" expression
flatten.params <- lapply(flatten.params,function(e) if(!as.character(e[[2]][[1]]) %in% prior.families) e[[2]][[2]] else e[[2]])
## transform to list
parse.prior <- function(prior) {
## declare here all the priors
Beta <- B <- Be <- function(a,b) list(name="Beta.prior",params=c(a,b))
Gamma <- G <- function(a,s) list(name="Gamma.prior",params=c(a,s))
Unif <- U <- function(a=0,b=1) list(name="Unif.prior",params=c(a,b))
Norm <- N <- function(m=0,s=1) list(name="Norm.prior",params=c(m,s))
LNorm <- LogNorm <- LN <- function(m=0,s=1) list(name="LNorm.prior",params=c(m,s))
NonInform <- NInf <- NI <- function(init=1,init_sigma=1) list(name="NonInform.prior",params=c(init,init_sigma))
res <- eval(prior) ## TODO or NOT TODO: eval(prior,parent.frame())
class(res) <- res$name #to be accessible as a class in R
res
}
flatten.params <- lapply(flatten.params,parse.prior)
return(flatten.params)
} else {
warning("Not a formula for bayesian.vam object!")
NULL # means not ok!
}
}
### FOUND A BETTER WAY: see bayesian.vam
## substitute prior with mean prior to be used inside mle.vam as initialization of run.bayesian.vam
# bayesian.model.to.mle.model <- function(model,priors) {
# print(priors)
# if(!is.null(priors)) {
# k<-0
# for(i in 1:length(model$family$params)) {
# ##DEBUG: print("prior");print(k+1);print(priors[[k+1]]);print(update(priors[[k+1]]))
# model$family$params[[i]] <- update(priors[[k<-k+1]])
# }
# model$family$params <- unlist(model$family$params)
# for(j in 1:length(model$models)) {
# for(i in 1:length(model$models[[j]]$params)) {
# ##DEBUG: print("prior");print(k+1);print(priors[[k+1]]);print(update(priors[[k+1]]))
# model$models[[j]]$params[[i]] <- update(priors[[k<-k+1]])
# }
# model$models[[j]]$params <- unlist(model$models[[j]]$params)
# }
# return(model)
# }
# return(NULL)
# }
## RMK: these 2 init following functions were only used for mle but are useful now for Bayesian too!
init.par0 <- function(obj,par0) {
## save the initial param
if(is.null(obj$par0)) obj$par0 <- params(obj)
## parameters stuff!
if(missing(par0)) {
if("par" %in% names(obj)) {
param <- obj$par[-1]
alpha <- obj$par[1] #LD2
} else {#not the first run
param<-params(obj)[-1] #first run
alpha<-params(obj)[1] #LD2
}
} else if(is.null(par0)) {
param<-obj$par0[-1]
alpha<-obj$par0[1] #LD2
} else {
param<-par0[-1]
alpha<-par0[1] #LD2
}
list(param=param,alpha=alpha)
}
init.fixed.param <- function(param,fixed) {
## fixed and functions stuff!
if(missing(fixed)) {
fixed<-rep(FALSE,length(param))
alpha_fixed<-FALSE #LD2
} else if(is.numeric(fixed)) {
fixedInd<-fixed
fixed<-rep(FALSE,length(param))
fixed[fixedInd-1]<-TRUE #LD2: ajout du -1
alpha_fixed<-sum(fixedInd==1) #OR ( 1 %in% fixedInd which is more expressive)
} else {#LD2
alpha_fixed<-fixed[1] #LD2
fixed<-fixed[-1] #LD2
}#LD2
list(fixed=fixed,alpha_fixed=alpha_fixed)
}
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.