Nothing
###################################################
### chunk number 1:
###################################################
algo.hmm <- function(disProgObj, control = list(range=range, Mtilde=-1, noStates=2, trend=TRUE, noHarmonics=1,covEffectEqual=FALSE, saveHMMs = FALSE, extraMSMargs=list() )){
# check if the msm package is available
if (!requireNamespace("msm")) {
stop("the HMM method requires package ", sQuote("msm"))
}
# Set the default values if not yet set
if(is.null(control$Mtilde)){ control$Mtilde <- -1 }
if(is.null(control$noStates)){ control$noStates <- 2 }
if(is.null(control$trend)){ control$trend <- TRUE }
if(is.null(control$noHarmonics)){ control$noHarmonics <- 1 }
if(is.null(control$covEffectEqual)){ control$covEffectEqual <- FALSE }
if(is.null(control$saveHMMs)){ control$saveHMMs <- FALSE }
if(is.null(control$extraMSMargs)){ control$extraMSMargs <- list() }
#Stop if not enough for estimation
if(min(control$range) < 2) {
stop("Error: Too few values as reference values")
}
# initialize the necessary vectors
alarm <- matrix(data = 0, nrow = length(control$range), ncol = 1)
upperbound <- matrix(data = 0, nrow = length(control$range), ncol = 1)
control$hmms <- list()
##############################################
#Repeat for each time point to monitor on-line
##############################################
for (i in 1:length(control$range)) {
#Function is so slow some sort of performance indicator is usually necessary
cat(paste("i=",i," (out of ",length(control$range),")\n",sep=""))
#Initialize observations for each round -- can be done sequentally
first <- ifelse(control$Mtilde== -1, 1, max(control$range[i]-control$Mtilde+1,1))
t <- first:control$range[i]
observed <- disProgObj$observed[t]
#Init data
counts <- data.frame(observed, t)
names(counts) <- c("observed","t")
#Initialize formula
formulaStr <- ifelse(control$trend, "~ 1 + t ", "~ 1 ")
#Create formula and add harmonics as covariates -- this could be done recursively?
for (j in seq_len(control$noHarmonics)) {
counts[,paste("cos",j,"t",sep="")] <- cos(2*j*pi*(t-1)/disProgObj$freq)
counts[,paste("sin",j,"t",sep="")] <- sin(2*j*pi*(t-1)/disProgObj$freq)
formulaStr <- paste(formulaStr,"+ cos",j,"t + sin",j,"t ",sep="")
}
#Obtain crude inits
q <- quantile(observed,seq(0,1,length.out=control$noStates+1))
lvl <- cut(observed,breaks=q,include.lowest=TRUE)
crudeMean <- as.numeric(tapply(observed, lvl, mean))
hcovariates <- list()
hmodel <- list()
for (j in seq_len(control$noStates)) {
hcovariates[[j]] <- as.formula(formulaStr)
val <- crudeMean[j]
#Substitution necessary, as hmmPois does lazy evaluation of rate argument
hmodel[[j]] <- eval(substitute(msm::hmmPois(rate=val),list(val=crudeMean[j])))
}
#Any constraints on the parameters of the covariates for the different states
hconstraint <- list()
if (control$covEffectEqual) {
hconstraint <- list(t=rep(1,control$noStates))
for (j in seq_len(control$noHarmonics)) {
hconstraint[[paste("sin",j,"t",sep="")]] <- rep(1,control$noStates)
hconstraint[[paste("cos",j,"t",sep="")]] <- rep(1,control$noStates)
}
}
#Prepare object for msm fitting
msm.args <- list(formula = observed ~ t,
data = counts,
#HMM with "noStates" states having equal initial values
qmatrix = matrix(1/control$noStates,control$noStates,control$noStates),
#y|x \sim Po( \mu[t] ) with some initial values
hmodel = hmodel,
#Models for \log \mu_t^1 and \log \mu_t^2
hcovariates = hcovariates,
#Force the effects of the trend and harmonics to be equal for all states
hconstraint=hconstraint)
#Add additional msm arguments
msm.args <- modifyList(msm.args, control$extraMSMargs)
# fit the HMM
hmm <- do.call(what=msm::msm, args=msm.args)
#In case the model fits should be saved.
if (control$saveHMMs) {
control$hmms[[i]] <- hmm
}
#If most probable state of current time point (i.e. last obs) equals the
#highest state then do alarm
# print(observed)
# print(matrix(viterbi.msm(hmm)$fitted,ncol=1))
alarm[i] <- msm::viterbi.msm(hmm)$fitted[length(t)] == control$noStates
#Upperbound does not have any meaning -- compute posterior probability!
upperbound[i] <- 0
}
#Add name and data name to control object.
control$name <- paste("hmm:", control$trans)
control$data <- paste(deparse(substitute(disProgObj)))
#no need for hmm object -- control$hmm <- hmm
# return alarm and upperbound vectors
result <- list(alarm = alarm, upperbound = upperbound, disProgObj=disProgObj,control=control)
class(result) = "survRes" # for surveillance system result
return(result)
}
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.