Nothing
#' ErlangHMM for MAP
#'
#' A special case of MAP.
#'
ERHMMClass <- R6::R6Class(
"ERHMMClass",
private = list(
param.alpha = NULL,
param.shape = NULL,
param.rate = NULL,
param.P = NULL,
param.xi = NULL,
matclass = "dgeMatrix",
make.matrix = function() {
as(as.matrix(private$param.P), private$matclass)
}
),
public = list(
#' @description
#' Get alpha
#' @return A vector of alpha
alpha = function() private$param.alpha,
#' @description
#' Get shape
#' @return A vector of shapes
shape = function() private$param.shape,
#' @description
#' Get rate
#' @return A vector of rates
rate = function() private$param.rate,
#' @description
#' Get P
#' @return A matrix of P
P = function() private$param.P,
#' @description
#' Get exit rates
#' @return A vector of exit rates
xi = function() private$param.xi,
#' @description
#' Create an ERHMM
#' @param alpha A vector of initial probability
#' @param shape A vector of shape parameters
#' @param rate A vector of rate parameters
#' @param P A matrix of transition probabilities
#' @param xi An exit rate vector
#' @return An instance of ERHMM
initialize = function(alpha, shape, rate, P, xi) {
# check
if (!(length(alpha) == length(shape)
&& length(alpha) == length(rate)
&& length(alpha) == dim(P)[1]
&& length(alpha) == dim(P)[2]
&& length(alpha) == length(xi))) {
stop(sprintf("Error: vector and matrix size are wrong"))
}
private$param.alpha <- alpha
private$param.shape <- shape
private$param.rate <- rate
private$param.P <- as(P, private$matclass)
private$param.xi <- xi
},
#' @description
#' copy
#' @return A new instance
copy = function() {
erhmm(shape=as.vector(self$shape()),
alpha=as.vector(self$alpha()),
rate = as.vector(self$rate()),
P = as.matrix(self$P()))
},
#' @description
#' The number of components
#' @return The number of components
size = function() {
length(self$alpha())
},
#' @description
#' Degrees of freedom
#' @return The degrees of freedom
df = function() {
self$size() * self$size() - self$size() + 2*self$size() - 1
},
#' @description
#' Print
#' @param ... Others
print = function(...) {
cat(gettextf("Size : %d\n", self$size()))
cat("Shape : ", self$shape(), "\n")
cat("Initial : ", self$alpha(), "\n")
cat("Rate : ", self$rate(), "\n")
cat("Transition probability : \n")
print(self$P())
},
#' @description
#' Marginal moments
#' @param k An integer of degree
#' @param ... Others
#' @return A vector of moments
mmoment = function(k, ...) {
m <- as.map(self)
m$mmoment(k, ...)
},
#' @description
#' Joint moments
#' @param lag An integer of lag
#' @param ... Others
#' @return A matrix of moments
jmoment = function(lag, ...) {
m <- as.map(self)
m$jmoment(lag, ...)
},
#' @description
#' k-lag correlation
#' @param ... Others
#' @return A vector for k-lag correlation
acf = function(...) {
m <- as.map(self)
m$acf(...)
},
#' @description
#' Run EM
#' @param data A dataframe
#' @param options A list of options
#' @param ... Others
emfit = function(data, options, ...) {
con <- list(...)
nmsC <- names(options)
options[(namc <- names(con))] <- con
alpha <- self$alpha()
shape <- self$shape()
rate <- self$rate()
xi <- self$xi()
P <- self$P()
H <- private$make.matrix()
switch(class(data),
"map.time" = emfit_erhmm_time(alpha, xi, rate, shape, P, data, options, H),
stop("Error: ERHMM can handle time point data only.")
)
},
#' @description
#' Initialize with data
#' @param data A dataframe
#' @param ... Others
init = function(data, ...) {
m <- erhmm.param(data, self, ...)
self$initialize(m$alpha(), m$shape(), m$rate(), m$P(), m$xi())
}
)
)
#' ErlangHMM for MAP with fixed phases
#'
#' A special case of MAP.
#'
AERHMMClass <- R6::R6Class(
"AERHMMClass",
private = list(
param.size = NULL,
param.erhmm = NULL
),
public = list(
#' @description
#' Get alpha
#' @return A vector of alpha
alpha = function() private$param.erhmm$alpha(),
#' @description
#' Get shape
#' @return A vector of shapes
shape = function() private$param.erhmm$shape(),
#' @description
#' Get rate
#' @return A vector of rates
rate = function() private$param.erhmm$rate(),
#' @description
#' Get P
#' @return A matrix of P
P = function() private$param.erhmm$P(),
#' @description
#' Get exit rates
#' @return A vector of exit rates
xi = function() private$param.erhmm$xi(),
#' @description
#' Create an AERHMM
#' @param size An integer of the number of phases
#' @param erhmm An instance of ERHMM
#' @return An instance of AERHMM
initialize = function(size, erhmm) {
if (size != sum(erhmm$shape())) {
stop(sprintf("Error: size and erhmm are not matched."))
}
private$param.size <- size
private$param.erhmm <- erhmm
},
#' @description
#' copy
#' @return A new instance
copy = function() {
AERHMMClass$new(private$param.size, private$param.erhmm$copy())
},
#' @description
#' The number of components
#' @return The number of components
size = function() {
private$param.erhmm$size()
},
#' @description
#' Degrees of freedom
#' @return The degrees of freedom
df = function() {
self$size() * self$size() - self$size() + 3*self$size() - 2
},
#' @description
#' Print
#' @param ... Others
print = function(...) {
private$param.erhmm$print(...)
},
#' @description
#' Marginal moments
#' @param k An integer of degree
#' @param ... Others
#' @return A vector of moments
mmoment = function(k, ...) {
private$param.erhmm$mmoment(k, ...)
},
#' @description
#' Joint moments
#' @param lag An integer of lag
#' @param ... Others
#' @return A matrix of moments
jmoment = function(lag, ...) {
private$param.erhmm$jmoment(lag, ...)
},
#' @description
#' k-lag correlation
#' @param ... Others
#' @return A vector for k-lag correlation
acf = function(...) {
private$param.erhmm$acf(...)
},
#' @description
#' Run EM
#' @param data A dataframe
#' @param options A list of options
#' @param ... Others
emfit = function(data, options, ...) {
method <- options$shape.method
phsize <- private$param.size
lbound <- options$lbound
ubound <- ifelse(is.na(options$ubound), phsize, options$ubound)
verbose <- options$shape.verbose
result <- switch(
method,
"all" = {
maxllf <- -Inf
maxmodel <- NULL
maxresult <- NULL
allshape <- shape.all(phsize, lbound, ubound)
for (s in allshape) {
m <- erhmm(shape=s)
m$init(data)
result <- m$emfit(data, options, ...)
if (verbose)
cat("shape: ", s, gettextf(" llf=%.2f\n", result$llf))
if (is.finite(result$llf)) {
if (result$llf > maxllf) {
maxllf <- result$llf
maxresult <- result
maxmodel <- m
}
}
}
private$param.erhmm <- maxmodel
maxresult
},
"increment" = {
maxllf <- -Inf
maxmodel <- NULL
maxresult <- NULL
shape <- rep(1, lbound)
repeat {
shapelist1 <- shape.increment(shape, ubound, phsize)
shapelist2 <- shape.decrement(shape, lbound)
shapelist <- c(shapelist1, shapelist2)
shape <- NULL
for (s in shapelist) {
m <- erhmm(shape=s)
m$init(data)
result <- m$emfit(data, options, ...)
if (verbose)
cat("shape: ", s, gettextf(" llf=%.2f\n", result$llf))
if (is.finite(result$llf)) {
if (result$llf > maxllf) {
maxllf <- result$llf
maxresult <- result
shape <- s
maxmodel <- m
}
}
}
if (is.null(shape)) {
## warning(message="break")
break
}
}
private$param.erhmm <- maxmodel
maxresult
}
)
# print(result)
},
#' @description
#' Initialize with data
#' @param data A dataframe
#' @param options A list of options
#' @param ... Others
init = function(data, ...) {
private$param.erhmm$init(data, ...)
}
)
)
#' Create ERHMM
#'
#' Create an instance of ERHMM
#'
#' @param size An integer of the number of phases
#' @param shape A vector of shape parameters
#' @param alpha A vector of initial probability (alpha)
#' @param rate A vector of rate parameters
#' @param P A matrix of transition probabilities
#' @return An instance of ERHMM
#'
#' @note
#' If shape is given, shape is used even though size is set.
#'
#' @export
erhmm <- function(size, shape, alpha = rep(1/length(shape),length(shape)),
rate = rep(1,length(shape)), P = matrix(1/length(shape), length(shape), length(shape))) {
if (missing(shape)) {
if (missing(size)) {
stop("Either size or shape is needed.")
} else {
# create an erhmm
AERHMMClass$new(size, erhmm(shape=c(size)))
}
} else {
if (missing(size)) {
size <- length(shape)
}
ERHMMClass$new(alpha=alpha, shape=shape, rate=rate, P=P, xi=rep(1, size))
}
}
#' Convert from ERHMM to MAP
#'
#' Convert from ERHMM to the general MAP
#'
#' @param x An instance of ERHMM
#' @return An instance of MAP
#'
#' @examples
#' ## create a hyper Erlang with specific parameters
#' (param <- erhmm(shape=c(2,3), alpha=c(0.3,0.7), rate=c(1.0,10.0)))
#'
#' ## convert to a general PH
#' as.map(param)
#'
#' @export
as.map <- function(x) {
size <- x$size()
mixrate <- x$alpha()
shape <- x$shape()
rate <- x$rate()
P <- x$P()
phsize <- sum(shape)
index <- cumsum(shape)
sindex <- c(1, index + 1)[1:size]
eindex <- index
alpha <- numeric(phsize)
alpha[sindex] <- mixrate
v <- numeric(0)
i <- numeric(0)
j <- numeric(0)
for (k in 1:size) {
i <- c(i, sindex[k]:eindex[k])
j <- c(j, sindex[k]:eindex[k])
v <- c(v, rep(-rate[k], shape[k]))
}
for (k in 1:size) {
if (shape[k] != 1) {
i <- c(i, sindex[k]:(eindex[k]-1))
j <- c(j, (sindex[k]+1):eindex[k])
v <- c(v, rep(rate[k], shape[k]-1))
}
}
D0 <- sparseMatrix(dims=c(phsize,phsize), i=i, j=j, x=v)
v <- numeric(0)
i <- numeric(0)
j <- numeric(0)
for (k1 in 1:size) {
for (k2 in 1:size) {
i <- c(i, eindex[k1])
j <- c(j, sindex[k2])
v <- c(v, rate[k1] * P[k1,k2])
}
}
D1 <- sparseMatrix(dims=c(phsize,phsize), i=i, j=j, x=v)
map(alpha=alpha, D0=D0, D1=D1)
}
#' Determine ERHMM parameters
#'
#' Determine ERHMM parameters with k-means.
#'
#' @param data A dataframe
#' @param skel An instance of ERHMM used as a skeleton
#' @param ... Others
#' @return An instance of ERHMM
erhmm.param <- function(data, skel, ...) {
size <- skel$size()
shape <- skel$shape()
tmp <- numeric(size)
rate <- numeric(size)
if (size >= 2) {
x <- data$intervals
result <- kmeans(x, size)
for (k in 1:size) {
m <- base::mean(x[result$cluster == k])
s2 <- var(x[result$cluster == k])
tmp[k] <- round(m^2 / s2)
rate[k] <- 1.0 / m
}
rate <- rate[rank(tmp)] * shape
} else {
m <- mean(data)
rate[1] <- shape[1] / m
}
alpha <- skel$alpha() * runif(size)
alpha <- alpha / sum(alpha)
P <- as.matrix(skel$P()) * matrix(runif(size*size), size, size)
P <- P / apply(P, 1, sum)
erhmm(shape=shape, alpha=alpha, rate=rate, P=P)
}
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.