Nothing
#' PH fitting with point data
#'
#' Estimates PH parameters from point data.
#'
#' @param ph An object of R6 class for PH. The estimation algorithm is selected depending on this class.
#' @param x A vector for point data.
#' @param weights A vector of weights for points.
#' @param ... Further options for fitting methods.
#' @return
#' Returns a list with components, which is an object of S3 class \code{phfit.result};
#' \item{model}{an object for estimated PH class.}
#' \item{llf}{a value of the maximum log-likelihood.}
#' \item{df}{a value of degrees of freedom of the model.}
#' \item{aic}{a value of Akaike information criterion.}
#' \item{iter}{the number of iterations.}
#' \item{convergence}{a logical value for the convergence of estimation algorithm.}
#' \item{ctime}{computation time (user time).}
#' \item{data}{an object for data class}
#' \item{aerror}{a value of absolute error for llf at the last step of algorithm.}
#' \item{rerror}{a value of relative error for llf at the last step of algorithm.}
#' \item{options}{a list of options used for fitting.}
#' \item{call}{the matched call.}
#'
#' @examples
#' ## make sample
#' wsample <- rweibull(n=100, shape=2, scale=1)
#'
#' ## PH fitting for general PH
#' (result1 <- phfit.point(ph=ph(2), x=wsample))
#'
#' ## PH fitting for CF1
#' (result2 <- phfit.point(ph=cf1(2), x=wsample))
#'
#' ## PH fitting for hyper Erlang
#' (result3 <- phfit.point(ph=herlang(3), x=wsample))
#'
#' ## mean
#' ph.mean(result1$model)
#' ph.mean(result2$model)
#' ph.mean(result3$model)
#'
#' ## variance
#' ph.var(result1$model)
#' ph.var(result2$model)
#' ph.var(result3$model)
#'
#' ## up to 5 moments
#' ph.moment(5, result1$model)
#' ph.moment(5, result2$model)
#' ph.moment(5, result3$model)
#'
#' @export
phfit.point <- function(ph, x, weights, ...) {
call <- match.call()
ph <- ph$copy()
options <- emoptions()
con <- list(...)
nmsC <- names(options)
options[(namc <- names(con))] <- con
if (length(noNms <- namc[!namc %in% nmsC]))
warning("unknown names in control: ", paste(noNms, collapse = ", "))
data <- data.frame.phase.time(x=x, weights=weights)
if (options$initialize == TRUE) {
ph$init(data, options)
}
tres <- system.time(result <- ph$emfit(data, options, ...))
result <- c(result, list(model=ph, aic=-2*(result$llf - ph$df()), df=ph$df(),
data=data, ctime=tres[1], options=options, call=call))
class(result) <- "phfit.result"
result
}
#' PH fitting with grouped data
#'
#' Estimates PH parameters from grouped data.
#'
#' @param ph An object of R6 class. The estimation algorithm is selected depending on this class.
#' @param counts A vector of the number of points in intervals.
#' @param breaks A vector for a sequence of points of boundaries of intervals.
#' This is equivalent to \code{c(0,cumsum(intervals))}.
#' If this is missing, it is assigned to \code{0:length(counts)}.
#' @param intervals A vector of time lengths for intervals.
#' This is equivalent to \code{diff(breaks)}).
#' If this is missing, it is assigned to \code{rep(1,length(counts))}.
#' @param instants A vector of integers to indicate whether sample is drawn at
#' the last of interval. If instant is 1, a sample is drawn at the last of interval.
#' If instant is 0, no sample is drawn at the last of interval.
#' By using instant, point data can be expressed by grouped data.
#' If instant is missing, it is given by \code{rep(0L,length(counts))}, i.e.,
#' there are no samples at the last of interval.
#' @param ... Further options for EM steps.
#' @return
#' Returns a list with components, which is an object of S3 class \code{phfit.result};
#' \item{model}{an object for estimated PH class.}
#' \item{llf}{a value of the maximum log-likelihood.}
#' \item{df}{a value of degrees of freedom of the model.}
#' \item{aic}{a value of Akaike information criterion.}
#' \item{iter}{the number of iterations.}
#' \item{convergence}{a logical value for the convergence of estimation algorithm.}
#' \item{ctime}{computation time (user time).}
#' \item{data}{an object for data class}
#' \item{aerror}{a value of absolute error for llf at the last step of algorithm.}
#' \item{rerror}{a value of relative error for llf at the last step of algorithm.}
#' \item{options}{a list of options used in the fitting.}
#' \item{call}{the matched call.}
#'
#' @note
#' In this method, we can handle truncated data using \code{NA} and \code{Inf};
#' \code{phfit.group(ph=cf1(5), counts=c(countsdata, NA), breaks=c(breakdata, +Inf))}
#' \code{NA} means missing of count data at the corresponding interval, and \code{Inf} is allowed to put
#' the last of breaks or intervals which represents a special interval [the last break point,infinity).
#'
#' @examples
#' ## make sample
#' wsample <- rweibull(n=100, shape=2, scale=1)
#' wgroup <- hist(x=wsample, breaks="fd", plot=FALSE)
#'
#' ## PH fitting for general PH
#' (result1 <- phfit.group(ph=ph(2), counts=wgroup$counts, breaks=wgroup$breaks))
#'
#' ## PH fitting for CF1
#' (result2 <- phfit.group(ph=cf1(2), counts=wgroup$counts, breaks=wgroup$breaks))
#'
#' ## PH fitting for hyper Erlang
#' (result3 <- phfit.group(ph=herlang(3), counts=wgroup$counts, breaks=wgroup$breaks))
#'
#' ## mean
#' ph.mean(result1$model)
#' ph.mean(result2$model)
#' ph.mean(result3$model)
#'
#' ## variance
#' ph.var(result1$model)
#' ph.var(result2$model)
#' ph.var(result3$model)
#'
#' ## up to 5 moments
#' ph.moment(5, result1$model)
#' ph.moment(5, result2$model)
#' ph.moment(5, result3$model)
#'
#' @export
phfit.group <- function(ph, counts, breaks, intervals, instants, ...) {
call <- match.call()
ph <- ph$copy()
options <- emoptions()
con <- list(...)
nmsC <- names(options)
options[(namc <- names(con))] <- con
if (length(noNms <- namc[!namc %in% nmsC]))
warning("unknown names in control: ", paste(noNms, collapse = ", "))
data <- data.frame.phase.group(counts=counts, breaks=breaks,
intervals=intervals, instants=instants)
if (options$initialize == TRUE) {
ph$init(data, options)
}
tres <- system.time(result <- ph$emfit(data, options, ...))
result <- c(result, list(model=ph, aic=-2*(result$llf - ph$df()), df=ph$df(),
data=data, ctime=tres[1], options=options, call=call))
class(result) <- "phfit.result"
result
}
#' PH fitting with density function
#'
#' Estimates PH parameters from density function.
#'
#' @param ph An object of R6 class. The estimation algorithm is selected depending on this class.
#' @param f A function object for a density function.
#' @param deformula An object for formulas of numerical integration.
#' It is not necessary to change it when the density function is defined on
#' the positive domain [0,infinity).
#' @param weight.zero A absolute value which is regarded as zero in numerical integration.
#' @param weight.reltol A value for precision of numerical integration.
#' @param start.divisions A value for starting value of divisions in deformula.
#' @param max.iter A value for the maximum number of iterations to increase divisions in deformula.
#' @param ... Options for EM steps, which is also used to send the arguments to density function.
#'
#' @return
#' Returns a list with components, which is an object of S3 class \code{phfit.result};
#' \item{model}{an object for estimated PH class.}
#' \item{llf}{a value of the maximum log-likelihood (a negative value of the cross entropy).}
#' \item{df}{a value of degrees of freedom of the model.}
#' \item{KL}{a value of Kullback-Leibler divergence.}
#' \item{iter}{the number of iterations.}
#' \item{convergence}{a logical value for the convergence of estimation algorithm.}
#' \item{ctime}{computation time (user time).}
#' \item{data}{an object for data class}
#' \item{aerror}{a value of absolute error for llf at the last step of algorithm.}
#' \item{rerror}{a value of relative error for llf at the last step of algorithm.}
#' \item{options}{a list of options.}
#' \item{call}{the matched call.}
#'
#' @note
#' Any of density function can be applied to the argument \code{f}, where
#' \code{f} should be defined \code{f <- function(x, ...)}.
#' The first argument of \code{f} should be an integral parameter.
#' The other parameters are set in the argument \code{...} of \code{phfit.density}.
#' The truncated density function can also be used directly.
#'
#' @examples
#' ####################
#' ##### truncated density
#' ####################
#'
#' ## PH fitting for general PH
#' (result1 <- phfit.density(ph=ph(2), f=dnorm, mean=3, sd=1))
#'
#' ## PH fitting for CF1
#' (result2 <- phfit.density(ph=cf1(2), f=dnorm, mean=3, sd=1))
#'
#' ## PH fitting for hyper Erlang
#' (result3 <- phfit.density(ph=herlang(3), f=dnorm, mean=3, sd=1))
#'
#' ## mean
#' ph.mean(result1$model)
#' ph.mean(result2$model)
#' ph.mean(result3$model)
#'
#' ## variance
#' ph.var(result1$model)
#' ph.var(result2$model)
#' ph.var(result3$model)
#'
#' ## up to 5 moments
#' ph.moment(5, result1$model)
#' ph.moment(5, result2$model)
#' ph.moment(5, result3$model)
#'
#' @export
phfit.density <- function(
ph, f, deformula = deformula.zeroinf, weight.zero = 1.0e-12,
weight.reltol = 1.0e-8, start.divisions = 8, max.iter = 12,
...) {
call <- match.call()
ph <- ph$copy()
x <- deformula(f, ..., zero.eps = weight.zero,
rel.tol = weight.reltol,
start.divisions = start.divisions, max.iter = max.iter)
ll <- sum(x$w * log(f(x$x, ...)))
data <- data.frame.phase.time(x=x$x, weights=x$w)
options <- emoptions()
con <- list(...)
nmsC <- names(options)
options[(namc <- names(con))] <- con
# if (length(noNms <- namc[!namc %in% nmsC]))
# warning("unknown names in control: ", paste(noNms, collapse = ", "))
if (options$initialize == TRUE) {
ph$init(data, options)
}
tres <- system.time(result <- ph$emfit(data, options, ...))
result <- c(result, list(model=ph, KL=x$h * (ll-result$llf), df=ph$df(),
data=data, ctime=tres[1], options=options, call=call))
class(result) <- "phfit.result.density"
result
}
#' @aliases phfit.point phfit.group
#' @export
print.phfit.result <- function (x, ...) {
cat("\n")
cat(sprintf("Maximum LLF: %f\n", x$llf))
cat(sprintf("DF: %d\n", x$df))
cat(sprintf("AIC: %f\n", x$aic))
cat(sprintf("Iteration: %d / %d\n", x$iter, x$options$maxiter))
cat(sprintf("Computation time (user): %f\n", x$ctime))
cat(sprintf("Convergence: %s\n", x$convergence))
cat(sprintf("Error (abs): %e (tolerance %e)\n", x$aerror, x$options$abstol))
cat(sprintf("Error (rel): %e (tolerance %e)\n", x$rerror, x$options$reltol))
cat("\n")
x$model$print(...)
cat("\n\n")
invisible(x)
}
#' @aliases phfit.density
#' @export
print.phfit.result.density <- function (x, ...) {
cat("\n")
cat(sprintf("Maximum LLF: %f\n", x$llf))
cat(sprintf("DF: %d\n", x$df))
cat(sprintf("KL: %f\n", x$KL))
cat(sprintf("Iteration: %d / %d\n", x$iter, x$options$maxiter))
cat(sprintf("Computation time (user): %f\n", x$ctime))
cat(sprintf("Convergence: %s\n", x$convergence))
cat(sprintf("Error (abs): %e (tolerance %e)\n", x$aerror, x$options$abstol))
cat(sprintf("Error (rel): %e (tolerance %e)\n", x$rerror, x$options$reltol))
cat("\n")
x$model$print(...)
cat("\n\n")
invisible(x)
}
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.