| Title: | PH/MAP Parameter Estimation |
|---|---|
| Description: | Estimation methods for phase-type distribution (PH) and Markovian arrival process (MAP) from empirical data (point and grouped data) and density function. The tool is based on the following researches: Okamura et al. (2009) <doi:10.1109/TNET.2008.2008750>, Okamura and Dohi (2009) <doi:10.1109/QEST.2009.28>, Okamura et al. (2011) <doi:10.1016/j.peva.2011.04.001>, Okamura et al. (2013) <doi:10.1002/asmb.1919>, Horvath and Okamura (2013) <doi:10.1007/978-3-642-40725-3_10>, Okamura and Dohi (2016) <doi:10.15807/jorsj.59.72>. |
| Authors: | Hiroyuki Okamura [aut, cre]
|
| Maintainer: | Hiroyuki Okamura <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.0 |
| Built: | 2024-11-09 06:31:19 UTC |
| Source: | CRAN |
Estimation methods for phase-type distribution (PH) and Markovian arrival process (MAP) from empirical data (point and grouped data) and density function. The tool is based on the following researches: Okamura et al. (2009) doi:10.1109/TNET.2008.2008750, Okamura and Dohi (2009) doi:10.1109/QEST.2009.28, Okamura et al. (2011) doi:10.1016/j.peva.2011.04.001, Okamura et al. (2013) doi:10.1002/asmb.1919, Horvath and Okamura (2013) doi:10.1007/978-3-642-40725-3_10, Okamura and Dohi (2016) doi:10.15807/jorsj.59.72.
Maintainer: Hiroyuki Okamura [email protected] (ORCID)
Useful links:
ErlangHMM for MAP with fixed phases
ErlangHMM for MAP with fixed phases
A special case of MAP.
alpha()
Get alpha
AERHMMClass$alpha()
A vector of alpha
shape()
Get shape
AERHMMClass$shape()
A vector of shapes
rate()
Get rate
AERHMMClass$rate()
A vector of rates
P()
Get P
AERHMMClass$P()
A matrix of P
xi()
Get exit rates
AERHMMClass$xi()
A vector of exit rates
new()
Create an AERHMM
AERHMMClass$new(size, erhmm)
sizeAn integer of the number of phases
erhmmAn instance of ERHMM
An instance of AERHMM
copy()
copy
AERHMMClass$copy()
A new instance
size()
The number of components
AERHMMClass$size()
The number of components
df()
Degrees of freedom
AERHMMClass$df()
The degrees of freedom
print()
AERHMMClass$print(...)
...Others
mmoment()
Marginal moments
AERHMMClass$mmoment(k, ...)
kAn integer of degree
...Others
A vector of moments
jmoment()
Joint moments
AERHMMClass$jmoment(lag, ...)
lagAn integer of lag
...Others
A matrix of moments
acf()
k-lag correlation
AERHMMClass$acf(...)
...Others
A vector for k-lag correlation
emfit()
Run EM
AERHMMClass$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
init()
Initialize with data
AERHMMClass$init(data, ...)
dataA dataframe
...Others
optionsA list of options
clone()
The objects of this class are cloneable with this method.
AERHMMClass$clone(deep = FALSE)
deepWhether to make a deep clone.
Hyper-Erlang distribution with a fixed phase
Hyper-Erlang distribution with a fixed phase
A mixture of Erlang distributions. A subclass of PH distributions.
mixrate()
Get mixrate
AHerlangClass$mixrate()
A vector of mixrate
shape()
Get shape
AHerlangClass$shape()
A vector of shapes
rate()
Get rate
AHerlangClass$rate()
A vector of rates
new()
Create a hyper-Erlang distribution with fixed phases
AHerlangClass$new(size, herlang)
sizeAn integer of the number of phases
herlangAn instance of HErlang
An instance of AHerlang
copy()
copy
AHerlangClass$copy()
A new instance
size()
The number of components
AHerlangClass$size()
The number of components
df()
Degrees of freedom
AHerlangClass$df()
The degrees of freedom
moment()
Moments of HErlang
AHerlangClass$moment(k, ...)
kA value to indicate the degrees of moments. k-th moment
...Others
A vector of moments from 1st to k-th moments
print()
AHerlangClass$print(...)
...Others
pdf()
AHerlangClass$pdf(x, ...)
xA vector of points
...Others
A vector of densities.
cdf()
CDF
AHerlangClass$cdf(q, ...)
qA vector of points
...Others
A vector of probabilities
ccdf()
Complementary CDF
AHerlangClass$ccdf(q, ...)
qA vector of points
...Others
A vector of probabilities
sample()
Make a sample
AHerlangClass$sample(...)
...Others
A sample of HErlang
emfit()
Run EM
AHerlangClass$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
init()
Initialize with data
AHerlangClass$init(data, ...)
dataA dataframe
...Others
optionsA list of options
clone()
The objects of this class are cloneable with this method.
AHerlangClass$clone(deep = FALSE)
deepWhether to make a deep clone.
Convert from hyper-Erlang distribution to the general PH distribution
as.gph(h)as.gph(h)
h |
An instance of HErlang |
An instance of GPH
#' ## create a hyper Erlang with specific parameters (param <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## convert to a general PH as.gph(param)#' ## create a hyper Erlang with specific parameters (param <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## convert to a general PH as.gph(param)
Convert from ERHMM to the general MAP
as.map(x)as.map(x)
x |
An instance of ERHMM |
An instance of MAP
## 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)## 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)
The data contains packet arrivals seen on an Ethernet at the Bellcore Morristown Research and Engineering facility. Two of the traces are LAN traffic (with a small portion of transit WAN traffic), and two are WAN traffic. The original trace BC-pAug89 began at 11:25 on August 29, 1989, and ran for about 3142.82 seconds (until 1,000,000 packets had been captured). The trace BC-pOct89 began at 11:00 on October 5, 1989, and ran for about 1759.62 seconds. These two traces captured all Ethernet packets. The number of arrivals in the original trace is one million.
BCpAug89 is a vector for the inter-arrival time in seconds for 1000 arrivals.
The original trace data are published in http://ita.ee.lbl.gov/html/contrib/BC.html.
Create an instance of CF1.
cf1(size, alpha, rate)cf1(size, alpha, rate)
size |
An integer of the number of phases |
alpha |
A vector of initial probabilities |
rate |
A vector of rates |
An instance of CF1.
## create a CF1 with 5 phases (param1 <- cf1(5)) ## create a CF1 with 5 phases (param1 <- cf1(size=5)) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0)))## create a CF1 with 5 phases (param1 <- cf1(5)) ## create a CF1 with 5 phases (param1 <- cf1(size=5)) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0)))
Crate CF1 with the first moment of a given data. This function calls cf1.param.linear and cf1.param.power to determine CF1. After execute 5 EM steps, the model with the smallest LLF is selected.
cf1.param(data, size, options, ...)cf1.param(data, size, options, ...)
data |
A dataframe |
size |
An integer for the number of phases |
options |
A list of options for EM steps |
... |
Others. This can provide additional options for EM steps. |
## Generate group data dat <- data.frame.phase.group(c(1,2,0,4), seq(0,10,length.out=5)) ## Create an instance of CF1 p <- cf1.param(data=dat, size=5)## Generate group data dat <- data.frame.phase.group(c(1,2,0,4), seq(0,10,length.out=5)) ## Create an instance of CF1 p <- cf1.param(data=dat, size=5)
Determine CF1 parameters based on the linear rule.
cf1.param.linear(size, mean, s)cf1.param.linear(size, mean, s)
size |
An integer of the number of phases |
mean |
A value of mean of data |
s |
A value of fraction of minimum and maximum rates |
A list of alpha and rate
Determine CF1 parameters based on the power rule.
cf1.param.power(size, mean, s)cf1.param.power(size, mean, s)
size |
An integer of the number of phases |
mean |
A value of mean of data |
s |
A value of fraction of minimum and maximum rates |
A list of alpha and rate
Canonical phase-type distribution
Canonical phase-type distribution
A continuous distribution dominated by a continuous-time Markov chain. A random time is given by an absorbing time. In the CF1 (canonical form 1), the infinitesimal generator is given by a bi-diagonal matrix, and whose order is determined by the ascending order.
mapfit::GPHClass -> CF1Class
rate()
Get rate
CF1Class$rate()
An instance of rate
new()
Create a CF1
CF1Class$new(alpha, rate)
alphaA vector of initial probability
rateA vector of rates
An instance of CF1
copy()
copy
CF1Class$copy()
A new instance
print()
CF1Class$print(...)
...Others
sample()
Generate a sample of CF1
CF1Class$sample(...)
...Others
A sample of CF1
emfit()
Run EM
CF1Class$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
init()
Initialize with data
CF1Class$init(data, options, ...)
dataA dataframe
optionsA list of options
...Others
clone()
The objects of this class are cloneable with this method.
CF1Class$clone(deep = FALSE)
deepWhether to make a deep clone.
Compute the stationary vector with GTH
ctmc.st(Q)ctmc.st(Q)
Q |
DTMC/CTMC kernel |
The stationary vector of DTMC/CTMC
Provide the data.frame for group data.
data.frame.map.group(counts, breaks, intervals, instants)data.frame.map.group(counts, breaks, intervals, instants)
counts |
A vector of the number of samples |
breaks |
A vector of break points |
intervals |
A vector of differences of time |
instants |
A vector meaning whether a sample is observed at the end of break. |
A dataframe
t <- c(1,1,1,1,1) n <- c(1,3,0,0,1) dat <- data.frame.map.group(counts=n, intervals=t) mean(dat) print(dat)t <- c(1,1,1,1,1) n <- c(1,3,0,0,1) dat <- data.frame.map.group(counts=n, intervals=t) mean(dat) print(dat)
Provide a data.frame with samples.
data.frame.map.time(time, intervals)data.frame.map.time(time, intervals)
time |
A vector for cumulative time |
intervals |
A vector for time intervals |
A dataframe
If both time and intervals are used, time is used.
map.time is given by a special case of map.group.
x <- runif(10) dat <- data.frame.map.time(time=x) mean(dat) print(dat)x <- runif(10) dat <- data.frame.map.time(time=x) mean(dat) print(dat)
Provide the data.frame for group data.
data.frame.phase.group(counts, breaks, intervals, instants)data.frame.phase.group(counts, breaks, intervals, instants)
counts |
A vector of the number of samples |
breaks |
A vector of break points |
intervals |
A vector of differences of time |
instants |
A vector meaning whether a sample is observed at the end of break. |
A dataframe
dat <- data.frame.phase.group(counts=c(1,2,1,1,0,0,1,4)) print(dat) mean(dat)dat <- data.frame.phase.group(counts=c(1,2,1,1,0,0,1,4)) print(dat) mean(dat)
Provide a data.frame with weighted samples.
data.frame.phase.time(x, weights)data.frame.phase.time(x, weights)
x |
A vector of point (quantiles) |
weights |
A vector of weights |
A dataframe
The point time is sorted and their differences are stored as the column of time
x <- runif(10) w <- runif(10) dat <- data.frame.phase.time(x=x, weights=w) print(dat) mean(dat)x <- runif(10) w <- runif(10) dat <- data.frame.phase.time(x=x, weights=w) print(dat) mean(dat)
Compute the probability density function (p.d.f.) for a given PH distribution
dphase(x, ph, log = FALSE, ...)dphase(x, ph, log = FALSE, ...)
x |
A numeric vector of quantiles. |
ph |
An instance of PH distribution. |
log |
logical; if TRUE, densities y are returned as log(y) |
... |
Others. |
A vector of densities.
## create a PH with specific parameters (phdist <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## p.d.f. for 0, 0.1, ..., 1 dphase(x=seq(0, 1, 0.1), ph=phdist)## create a PH with specific parameters (phdist <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## p.d.f. for 0, 0.1, ..., 1 dphase(x=seq(0, 1, 0.1), ph=phdist)
A list of options for EM
emoptions()emoptions()
A list of options with default values
Create an instance of ERHMM
erhmm( size, shape, alpha = rep(1/length(shape), length(shape)), rate = rep(1, length(shape)), P = matrix(1/length(shape), length(shape), length(shape)) )erhmm( size, shape, alpha = rep(1/length(shape), length(shape)), rate = rep(1, length(shape)), P = matrix(1/length(shape), length(shape), length(shape)) )
size |
An integer of the number of phases |
shape |
A vector of shape parameters |
alpha |
A vector of initial probability (alpha) |
rate |
A vector of rate parameters |
P |
A matrix of transition probabilities |
An instance of ERHMM
If shape is given, shape is used even though size is set.
Determine ERHMM parameters with k-means.
erhmm.param(data, skel, ...)erhmm.param(data, skel, ...)
data |
A dataframe |
skel |
An instance of ERHMM used as a skeleton |
... |
Others |
An instance of ERHMM
ErlangHMM for MAP
ErlangHMM for MAP
A special case of MAP.
alpha()
Get alpha
ERHMMClass$alpha()
A vector of alpha
shape()
Get shape
ERHMMClass$shape()
A vector of shapes
rate()
Get rate
ERHMMClass$rate()
A vector of rates
P()
Get P
ERHMMClass$P()
A matrix of P
xi()
Get exit rates
ERHMMClass$xi()
A vector of exit rates
new()
Create an ERHMM
ERHMMClass$new(alpha, shape, rate, P, xi)
alphaA vector of initial probability
shapeA vector of shape parameters
rateA vector of rate parameters
PA matrix of transition probabilities
xiAn exit rate vector
An instance of ERHMM
copy()
copy
ERHMMClass$copy()
A new instance
size()
The number of components
ERHMMClass$size()
The number of components
df()
Degrees of freedom
ERHMMClass$df()
The degrees of freedom
print()
ERHMMClass$print(...)
...Others
mmoment()
Marginal moments
ERHMMClass$mmoment(k, ...)
kAn integer of degree
...Others
A vector of moments
jmoment()
Joint moments
ERHMMClass$jmoment(lag, ...)
lagAn integer of lag
...Others
A matrix of moments
acf()
k-lag correlation
ERHMMClass$acf(...)
...Others
A vector for k-lag correlation
emfit()
Run EM
ERHMMClass$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
init()
Initialize with data
ERHMMClass$init(data, ...)
dataA dataframe
...Others
clone()
The objects of this class are cloneable with this method.
ERHMMClass$clone(deep = FALSE)
deepWhether to make a deep clone.
Create an instance of GMMPP
gmmpp(size, alpha, D0, D1)gmmpp(size, alpha, D0, D1)
size |
An integer for the number of phases |
alpha |
A vector of initial probability |
D0 |
An infinitesimal generator without arrivals |
D1 |
An infinitesimal generator with arrivals |
An instance of GMMPP
This function can omit several patterns of arguments. For example, map(5)
omit the arguments alpha, Q and xi. In this case, the default values are
assigned to them.
## create a map (full matrix) with 5 phases (param1 <- gmmpp(5)) ## create a map with specific parameters (param2 <- gmmpp(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), D1=rbind(c(2,0,0),c(0,2,0),c(0,0,0))))## create a map (full matrix) with 5 phases (param1 <- gmmpp(5)) ## create a map with specific parameters (param2 <- gmmpp(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), D1=rbind(c(2,0,0),c(0,2,0),c(0,0,0))))
GMMPP: Approximation for MAP
GMMPP: Approximation for MAP
A point process dominated by a continuous-time Markov chain.
mapfit::MAPClass -> GMMPPClass
new()
Create a MAP
GMMPPClass$new(alpha, D0, D1, xi)
alphaA vector of initial probability
D0An infinitesimal generator
D1An infinitesimal generator
xiAn exit rate vector
An instance of MAP
copy()
copy
GMMPPClass$copy()
A new instance
emfit()
Run EM
GMMPPClass$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
clone()
The objects of this class are cloneable with this method.
GMMPPClass$clone(deep = FALSE)
deepWhether to make a deep clone.
Generate GPH randomly and adjust parameters to fit its first moment to the first moment of data.
gph.param(data, skel, ...)gph.param(data, skel, ...)
data |
A dataframe |
skel |
An instance of skeleton of GPH. |
... |
Others |
An instance of GPH
## Create data wsample <- rweibull(10, shape=2) (dat <- data.frame.phase.time(x=wsample)) ## Generate PH that is fitted to dat (model <- gph.param(data=dat, skel=ph(5)))## Create data wsample <- rweibull(10, shape=2) (dat <- data.frame.phase.time(x=wsample)) ## Generate PH that is fitted to dat (model <- gph.param(data=dat, skel=ph(5)))
General phase-type distribution
General phase-type distribution
A continuous distribution dominated by a continuous-time Markov chain. A random time is given by an absorbing time.
alpha()
Get alpha
GPHClass$alpha()
A vector of alpha
Q()
Get Q
GPHClass$Q()
A matrix of Q
xi()
Get xi
GPHClass$xi()
A vector of xi
new()
Create a GPH
GPHClass$new(alpha, Q, xi)
alphaA vector of initial probability
QAn infinitesimal generator
xiAn exit rate vector
An instance of GPH
copy()
copy
GPHClass$copy()
A new instance
size()
The number of phases
GPHClass$size()
The number of phases
df()
Degrees of freedom
GPHClass$df()
The degrees of freedom
moment()
Moments of GPH
GPHClass$moment(k, ...)
kA value to indicate the degrees of moments. k-th moment
...Others
A vector of moments from 1st to k-th moments
print()
GPHClass$print(...)
...Others
pdf()
GPHClass$pdf(x, poisson.eps = 1e-08, ufactor = 1.01, ...)
xA vector of points
poisson.epsA value of tolerance error for uniformization
ufactorA value of uniformization factor
...Others
A vector of densities.
cdf()
CDF
GPHClass$cdf(x, poisson.eps = 1e-08, ufactor = 1.01, ...)
xA vector of points
poisson.epsA value of tolerance error for uniformization
ufactorA value of uniformization factor
...Others
A vector of probabilities
ccdf()
Complementary CDF
GPHClass$ccdf(x, poisson.eps = 1e-08, ufactor = 1.01, ...)
xA vector of points
poisson.epsA value of tolerance error for uniformization
ufactorA value of uniformization factor
...Others
A vector of probabilities
sample()
Make a sample
GPHClass$sample(...)
...Others
A sample of GPH
emfit()
Run EM
GPHClass$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
init()
Initialize with data
GPHClass$init(data, ...)
dataA dataframe
...Others
optionsA list of options
clone()
The objects of this class are cloneable with this method.
GPHClass$clone(deep = FALSE)
deepWhether to make a deep clone.
This function provides the values of p.d.f. for PH distribution with the uniformization technique.
This function provides the values of c.d.f. for PH distribution with the uniformization technique.
This function provides the values of complementary c.d.f. for PH distribution with the uniformization technique.
Create an instance of Hyper-Erlang distribution
herlang( size, shape, mixrate = rep(1/length(shape), length(shape)), rate = rep(1, length(shape)) )herlang( size, shape, mixrate = rep(1/length(shape), length(shape)), rate = rep(1, length(shape)) )
size |
An integer of the number of phases |
shape |
A vector of shape parameters |
mixrate |
A vector of initial probability (mixrate) |
rate |
A vector of rate parameters |
An instance of HErlang
If shape is given, shape is used even though size is set.
## create a hyper Erlang consisting of two Erlang ## with shape parameters 2 and 3. (param1 <- herlang(shape=c(2,3))) ## create a hyper Erlang with specific parameters (param2 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## convert to a general PH as.gph(param2) ## p.d.f. for 0, 0.1, ..., 1 (dphase(x=seq(0, 1, 0.1), ph=param2)) ## c.d.f. for 0, 0.1, ..., 1 (pphase(q=seq(0, 1, 0.1), ph=param2)) ## generate 10 samples (rphase(n=10, ph=param2))## create a hyper Erlang consisting of two Erlang ## with shape parameters 2 and 3. (param1 <- herlang(shape=c(2,3))) ## create a hyper Erlang with specific parameters (param2 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## convert to a general PH as.gph(param2) ## p.d.f. for 0, 0.1, ..., 1 (dphase(x=seq(0, 1, 0.1), ph=param2)) ## c.d.f. for 0, 0.1, ..., 1 (pphase(q=seq(0, 1, 0.1), ph=param2)) ## generate 10 samples (rphase(n=10, ph=param2))
Determine the hyper-Erlang parameters with k-means.
herlang.param(data, shape, ...)herlang.param(data, shape, ...)
data |
A dataframe |
shape |
A vector of shape parameters |
... |
Others |
An instance of HErlang
## Create data wsample <- rweibull(10, shape=2) (dat <- data.frame.phase.time(x=wsample)) ## Generate PH that is fitted to dat (model <- herlang.param(data=dat, shape=c(1,2,3)))## Create data wsample <- rweibull(10, shape=2) (dat <- data.frame.phase.time(x=wsample)) ## Generate PH that is fitted to dat (model <- herlang.param(data=dat, shape=c(1,2,3)))
Hyper-Erlang distribution
Hyper-Erlang distribution
A mixture of Erlang distributions. A subclass of PH distributions.
mixrate()
Get mixrate
HErlangClass$mixrate()
A vector of mixrate
shape()
Get shape
HErlangClass$shape()
A vector of shapes
rate()
Get rate
HErlangClass$rate()
A vector of rates
new()
Create a hyper-Erlang distribution
HErlangClass$new(mixrate, shape, rate)
mixrateA vector of initial probability
shapeA vector of shape parameters
rateA vector of rate parameters
An instance of HErlang
copy()
copy
HErlangClass$copy()
A new instance
size()
The number of components
HErlangClass$size()
The number of components
df()
Degrees of freedom
HErlangClass$df()
The degrees of freedom
moment()
Moments of HErlang
HErlangClass$moment(k, ...)
kA value to indicate the degrees of moments. k-th moment
...Others
A vector of moments from 1st to k-th moments
print()
HErlangClass$print(...)
...Others
pdf()
HErlangClass$pdf(x, ...)
xA vector of points
...Others
A vector of densities.
cdf()
CDF
HErlangClass$cdf(q, ...)
qA vector of points
...Others
A vector of probabilities
ccdf()
Complementary CDF
HErlangClass$ccdf(q, ...)
qA vector of points
...Others
A vector of probabilities
sample()
Make a sample
HErlangClass$sample(...)
...Others
A sample of HErlang
emfit()
Run EM
HErlangClass$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
init()
Initialize with data
HErlangClass$init(data, ...)
dataA dataframe
...Others
optionsA list of options
clone()
The objects of this class are cloneable with this method.
HErlangClass$clone(deep = FALSE)
deepWhether to make a deep clone.
Create an instance of MAP
map(size, alpha, D0, D1)map(size, alpha, D0, D1)
size |
An integer for the number of phases |
alpha |
A vector of initial probability |
D0 |
An infinitesimal generator without arrivals |
D1 |
An infinitesimal generator with arrivals |
An instance of MAP
This function can omit several patterns of arguments. For example, map(5)
omit the arguments alpha, D0 D1 and xi. In this case, the default values are
assigned to them.
## create a map (full matrix) with 5 phases (param1 <- map(5)) ## create a map with specific parameters (param2 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), D1=rbind(c(2,0,0),c(0,2,0),c(0,0,0))))## create a map (full matrix) with 5 phases (param1 <- map(5)) ## create a map with specific parameters (param2 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), D1=rbind(c(2,0,0),c(0,2,0),c(0,0,0))))
Compute k-lag correlation
map.acf(map, ...)map.acf(map, ...)
map |
An instance of MAP |
... |
Others |
A vector of k-lag correlation
## create an MAP with specific parameters (param1 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-4)), D1=rbind(c(1,1,0),c(1,0,1),c(2,0,1)))) ## create an ER-HMM with specific parameters (param2 <- erhmm(shape=c(2,3), alpha=c(0.3,0.7), rate=c(1.0,10.0), P=rbind(c(0.3, 0.7), c(0.1, 0.9)))) map.acf(map=param1) map.acf(map=param2)## create an MAP with specific parameters (param1 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-4)), D1=rbind(c(1,1,0),c(1,0,1),c(2,0,1)))) ## create an ER-HMM with specific parameters (param2 <- erhmm(shape=c(2,3), alpha=c(0.3,0.7), rate=c(1.0,10.0), P=rbind(c(0.3, 0.7), c(0.1, 0.9)))) map.acf(map=param1) map.acf(map=param2)
Compute joint moments for a given MAP
map.jmoment(lag, map, ...)map.jmoment(lag, map, ...)
lag |
An integer for lag |
map |
An instance of MAP |
... |
Others |
A vector of moments
## create an MAP with specific parameters (param1 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-4)), D1=rbind(c(1,1,0),c(1,0,1),c(2,0,1)))) ## create an ER-HMM with specific parameters (param2 <- erhmm(shape=c(2,3), alpha=c(0.3,0.7), rate=c(1.0,10.0), P=rbind(c(0.3, 0.7), c(0.1, 0.9)))) map.jmoment(lag=1, map=param1) map.jmoment(lag=1, map=param2)## create an MAP with specific parameters (param1 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-4)), D1=rbind(c(1,1,0),c(1,0,1),c(2,0,1)))) ## create an ER-HMM with specific parameters (param2 <- erhmm(shape=c(2,3), alpha=c(0.3,0.7), rate=c(1.0,10.0), P=rbind(c(0.3, 0.7), c(0.1, 0.9)))) map.jmoment(lag=1, map=param1) map.jmoment(lag=1, map=param2)
Compute up to k-th marginal moments for a given MAP
map.mmoment(k, map, ...)map.mmoment(k, map, ...)
k |
An integer for the moments to be computed |
map |
An instance of MAP |
... |
Others |
A vector of moments
## create an MAP with specific parameters (param1 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-4)), D1=rbind(c(1,1,0),c(1,0,1),c(2,0,1)))) ## create an ER-HMM with specific parameters (param2 <- erhmm(shape=c(2,3), alpha=c(0.3,0.7), rate=c(1.0,10.0), P=rbind(c(0.3, 0.7), c(0.1, 0.9)))) map.mmoment(k=3, map=param1) map.mmoment(k=3, map=param2)## create an MAP with specific parameters (param1 <- map(alpha=c(1,0,0), D0=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-4)), D1=rbind(c(1,1,0),c(1,0,1),c(2,0,1)))) ## create an ER-HMM with specific parameters (param2 <- erhmm(shape=c(2,3), alpha=c(0.3,0.7), rate=c(1.0,10.0), P=rbind(c(0.3, 0.7), c(0.1, 0.9)))) map.mmoment(k=3, map=param1) map.mmoment(k=3, map=param2)
Generate MAP randomly and adjust parameters to fit its first moment to the first moment of data.
map.param(data, skel, ...)map.param(data, skel, ...)
data |
A dataframe |
skel |
An instance of skeleton of MAP |
... |
Others |
An instance of MAP
General Markovian arrival process
General Markovian arrival process
A point process dominated by a continuous-time Markov chain.
alpha()
Get alpha
MAPClass$alpha()
A vector of alpha
D0()
Get D0
MAPClass$D0()
A matrix of D0
D1()
Get D1
MAPClass$D1()
A matrix of D1
xi()
Get xi
MAPClass$xi()
A vector of xi
new()
Create a MAP
MAPClass$new(alpha, D0, D1, xi)
alphaA vector of initial probability
D0An infinitesimal generator
D1An infinitesimal generator
xiAn exit rate vector
An instance of MAP
copy()
copy
MAPClass$copy()
A new instance
size()
The number of phases
MAPClass$size()
The number of phases
df()
Degrees of freedom
MAPClass$df()
The degrees of freedom
print()
MAPClass$print(...)
...Others
mmoment()
Marginal moments
MAPClass$mmoment(k, ...)
kAn integer of degree
...Others
A vector of moments
jmoment()
Joint moments
MAPClass$jmoment(lag, ...)
lagAn integer of lag
...Others
A matrix of moments
acf()
k-lag correlation
MAPClass$acf(...)
...Others
A vector for k-lag correlation
emfit()
Run EM
MAPClass$emfit(data, options, ...)
dataA dataframe
optionsA list of options
...Others
init()
Initialize with data
MAPClass$init(data, ...)
dataA dataframe
...Others
clone()
The objects of this class are cloneable with this method.
MAPClass$clone(deep = FALSE)
deepWhether to make a deep clone.
Estimates MAP parameters from grouped data.
mapfit.group(map, counts, breaks, intervals, instants, ...)mapfit.group(map, counts, breaks, intervals, instants, ...)
map |
An object of R6 class. The estimation algorithm is selected depending on this class. |
counts |
A vector of the number of points in intervals. |
breaks |
A vector for a sequence of points of boundaries of intervals.
This is equivalent to |
intervals |
A vector of time lengths for intervals.
This is equivalent to |
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 |
... |
Further options for EM steps. |
Returns a list with components, which is an object of S3 class mapfit.result;
model |
an object for estimated MAP class. |
llf |
a value of the maximum log-likelihood. |
df |
a value of degrees of freedom of the model. |
aic |
a value of Akaike information criterion. |
iter |
the number of iterations. |
convergence |
a logical value for the convergence of estimation algorithm. |
ctime |
computation time (user time). |
data |
an object for data class |
aerror |
a value of absolute error for llf at the last step of algorithm. |
rerror |
a value of relative error for llf at the last step of algorithm. |
options |
a list of options used in the fitting. |
call |
the matched call. |
## load trace data data(BCpAug89) BCpAug89s <- head(BCpAug89, 50) ## make grouped data BCpAug89.group <- hist(cumsum(BCpAug89s), breaks=seq(0, 0.15, 0.005), plot=FALSE) ## MAP fitting for general MAP (result1 <- mapfit.group(map=map(2), counts=BCpAug89.group$counts, breaks=BCpAug89.group$breaks)) ## MAP fitting for MMPP (result2 <- mapfit.group(map=mmpp(2), counts=BCpAug89.group$counts, breaks=BCpAug89.group$breaks)) ## MAP fitting with approximate MMPP (result3 <- mapfit.group(map=gmmpp(2), counts=BCpAug89.group$counts, breaks=BCpAug89.group$breaks)) ## marginal moments for estimated MAP map.mmoment(k=3, map=result1$model) map.mmoment(k=3, map=result2$model) map.mmoment(k=3, map=result3$model) ## joint moments for estimated MAP map.jmoment(lag=1, map=result1$model) map.jmoment(lag=1, map=result2$model) map.jmoment(lag=1, map=result3$model) ## lag-k correlation map.acf(map=result1$model) map.acf(map=result2$model) map.acf(map=result3$model)## load trace data data(BCpAug89) BCpAug89s <- head(BCpAug89, 50) ## make grouped data BCpAug89.group <- hist(cumsum(BCpAug89s), breaks=seq(0, 0.15, 0.005), plot=FALSE) ## MAP fitting for general MAP (result1 <- mapfit.group(map=map(2), counts=BCpAug89.group$counts, breaks=BCpAug89.group$breaks)) ## MAP fitting for MMPP (result2 <- mapfit.group(map=mmpp(2), counts=BCpAug89.group$counts, breaks=BCpAug89.group$breaks)) ## MAP fitting with approximate MMPP (result3 <- mapfit.group(map=gmmpp(2), counts=BCpAug89.group$counts, breaks=BCpAug89.group$breaks)) ## marginal moments for estimated MAP map.mmoment(k=3, map=result1$model) map.mmoment(k=3, map=result2$model) map.mmoment(k=3, map=result3$model) ## joint moments for estimated MAP map.jmoment(lag=1, map=result1$model) map.jmoment(lag=1, map=result2$model) map.jmoment(lag=1, map=result3$model) ## lag-k correlation map.acf(map=result1$model) map.acf(map=result2$model) map.acf(map=result3$model)
Estimates MAP parameters from point data.
mapfit.point(map, x, intervals, ...)mapfit.point(map, x, intervals, ...)
map |
An object for MAP. The estimation algorithm is selected depending on this class. |
x |
A vector for point data. |
intervals |
A vector for intervals. |
... |
Further options for fitting methods. |
Returns a list with components, which is an object of S3 class mapfit.result;
model |
an object for estimated PH class. |
llf |
a value of the maximum log-likelihood. |
df |
a value of degrees of freedom of the model. |
aic |
a value of Akaike information criterion. |
iter |
the number of iterations. |
convergence |
a logical value for the convergence of estimation algorithm. |
ctime |
computation time (user time). |
data |
an object for data class |
aerror |
a value of absolute error for llf at the last step of algorithm. |
rerror |
a value of relative error for llf at the last step of algorithm. |
options |
a list of options used for fitting. |
call |
the matched call. |
## load trace data data(BCpAug89) BCpAug89s <- head(BCpAug89, 50) ## MAP fitting for general MAP (result1 <- mapfit.point(map=map(2), x=cumsum(BCpAug89s))) ## MAP fitting for MMPP (result2 <- mapfit.point(map=mmpp(2), x=cumsum(BCpAug89s))) ## MAP fitting for ER-HMM (result3 <- mapfit.point(map=erhmm(3), x=cumsum(BCpAug89s))) ## marginal moments for estimated MAP map.mmoment(k=3, map=result1$model) map.mmoment(k=3, map=result2$model) map.mmoment(k=3, map=result3$model) ## joint moments for estimated MAP map.jmoment(lag=1, map=result1$model) map.jmoment(lag=1, map=result2$model) map.jmoment(lag=1, map=result3$model) ## lag-k correlation map.acf(map=result1$model) map.acf(map=result2$model) map.acf(map=result3$model)## load trace data data(BCpAug89) BCpAug89s <- head(BCpAug89, 50) ## MAP fitting for general MAP (result1 <- mapfit.point(map=map(2), x=cumsum(BCpAug89s))) ## MAP fitting for MMPP (result2 <- mapfit.point(map=mmpp(2), x=cumsum(BCpAug89s))) ## MAP fitting for ER-HMM (result3 <- mapfit.point(map=erhmm(3), x=cumsum(BCpAug89s))) ## marginal moments for estimated MAP map.mmoment(k=3, map=result1$model) map.mmoment(k=3, map=result2$model) map.mmoment(k=3, map=result3$model) ## joint moments for estimated MAP map.jmoment(lag=1, map=result1$model) map.jmoment(lag=1, map=result2$model) map.jmoment(lag=1, map=result3$model) ## lag-k correlation map.acf(map=result1$model) map.acf(map=result2$model) map.acf(map=result3$model)
Create an instance of MMPP (Markov-Modulated Poisson Process)
mmpp(size)mmpp(size)
size |
An integer for the number of phases |
An instance of MMPP
MMPP is a MAP whose D1 is given by a diagonal matrix.
## create an MMPP with 5 phases (param1 <- mmpp(5))## create an MMPP with 5 phases (param1 <- mmpp(5))
Create an instance of GPH
ph(size, alpha, Q, xi)ph(size, alpha, Q, xi)
size |
An integer for the number of phases |
alpha |
A vector of initial probability |
Q |
An infinitesimal generator |
xi |
An exit rate vector |
An instance of GPH
This function can omit several patterns of arguments. For example, ph(5)
omit the arguments alpha, Q and xi. In this case, the default values are
assigned to them.
## create a PH (full matrix) with 5 phases (param1 <- ph(5)) ## create a PH (full matrix) with 5 phases (param1 <- ph(size=5)) ## create a PH with specific parameters (param2 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0)))## create a PH (full matrix) with 5 phases (param1 <- ph(5)) ## create a PH (full matrix) with 5 phases (param1 <- ph(size=5)) ## create a PH with specific parameters (param2 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0)))
Create an instance of bi-diagonal PH distribution.
ph.bidiag(size)ph.bidiag(size)
size |
An integer for the number of phases |
An instance of bi-diagonal PH distribution
Bi-diagonal PH distribution is the PH distribution whose infinitesimal generator is given by a upper bi-diagonal matrix. This is similar to canonical form 1. But there is no restriction on the order for diagonal elements.
## create a bidiagonal PH with 5 phases (param1 <- ph.bidiag(5))## create a bidiagonal PH with 5 phases (param1 <- ph.bidiag(5))
Create an instance of coxian PH distribution.
ph.coxian(size)ph.coxian(size)
size |
An integer for the number of phases |
An instance of coxian PH distribution
Coxian PH distribution is the PH distribution whose infinitesimal generator is given by a upper bi-diagonal matrix. This is also called canonical form 3.
## create a Coxian PH with 5 phases (param1 <- ph.coxian(5))## create a Coxian PH with 5 phases (param1 <- ph.coxian(5))
Compute the mean of a given PH distribution.
ph.mean(ph, ...)ph.mean(ph, ...)
ph |
An instance of PH distribution |
... |
Others |
A value of mean
## create a PH with specific parameters (param1 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0))) ## create a hyper Erlang with specific parameters (param3 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## mean ph.mean(param1) ph.mean(param2) ph.mean(param3)## create a PH with specific parameters (param1 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0))) ## create a hyper Erlang with specific parameters (param3 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## mean ph.mean(param1) ph.mean(param2) ph.mean(param3)
Generate moments up to k-th moments for a given PH distribution.
ph.moment(k, ph, ...)ph.moment(k, ph, ...)
k |
An integer for the moments to be computed |
ph |
An instance of PH distribution |
... |
Others |
A vector of moments
## create a PH with specific parameters (param1 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0))) ## create a hyper Erlang with specific parameters (param3 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## up to 5 moments ph.moment(5, param1) ph.moment(5, param2) ph.moment(5, param3)## create a PH with specific parameters (param1 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0))) ## create a hyper Erlang with specific parameters (param3 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## up to 5 moments ph.moment(5, param1) ph.moment(5, param2) ph.moment(5, param3)
Create an instance of tri-diagonal PH distribution.
ph.tridiag(size)ph.tridiag(size)
size |
An integer for the number of phases |
An instance of tri-diagonal PH distribution
Tri-diagonal PH distribution is the PH distribution whose infinitesimal generator is given by a tri-diagonal matrix (band matrix).
## create a tridiagonal PH with 5 phases (param1 <- ph.tridiag(5))## create a tridiagonal PH with 5 phases (param1 <- ph.tridiag(5))
Compute the variance of a given PH distribution.
ph.var(ph, ...)ph.var(ph, ...)
ph |
An instance of PH distribution |
... |
Others |
A value of variance
## create a PH with specific parameters (param1 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0))) ## create a hyper Erlang with specific parameters (param3 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## variance ph.var(param1) ph.var(param2) ph.var(param3)## create a PH with specific parameters (param1 <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## create a CF1 with specific parameters (param2 <- cf1(alpha=c(1,0,0), rate=c(1.0,2.0,3.0))) ## create a hyper Erlang with specific parameters (param3 <- herlang(shape=c(2,3), mixrate=c(0.3,0.7), rate=c(1.0,10.0))) ## variance ph.var(param1) ph.var(param2) ph.var(param3)
Estimates PH parameters from three moments.
phfit.3mom( m1, m2, m3, method = c("Osogami06", "Bobbio05"), max.phase = 50, epsilon = sqrt(.Machine$double.eps) )phfit.3mom( m1, m2, m3, method = c("Osogami06", "Bobbio05"), max.phase = 50, epsilon = sqrt(.Machine$double.eps) )
m1 |
A value of the first moment. |
m2 |
A value of the second moment. |
m3 |
A value of the third moment. |
method |
The name of moment matching method. |
max.phase |
An integer for the maximum number of phases in the method "Osogami06". |
epsilon |
A value of precision in the method "Osogami06". |
An object of GPH.
The method "Osogami06" checks the first three moments on whether there exists a PH whose three moments match to them. In such case, the method "Bobbio05" often returns an error.
Osogami, T. and Harchol-Balter, M. (2006) Closed Form Solutions for Mapping General Distributions to Minimal PH Distributions. Performance Evaluation, 63(6), 524–552.
Bobbio, A., Horvath, A. and Telek, M. (2005) Matching Three Moments with Minimal Acyclic Phase Type Distributions. Stochastic Models, 21(2-3), 303–326.
## Three moment matching ## Moments of Weibull(shape=2, scale=1); (0.886227, 1.0, 1.32934) (result1 <- phfit.3mom(0.886227, 1.0, 1.32934)) ## Three moment matching ## Moments of Weibull(shape=2, scale=1); (0.886227, 1.0, 1.32934) (result2 <- phfit.3mom(0.886227, 1.0, 1.32934, method="Bobbio05")) ## mean ph.mean(result1) ph.mean(result2) ## variance ph.var(result1) ph.var(result2) ## up to 5 moments ph.moment(5, result1) ph.moment(5, result2)## Three moment matching ## Moments of Weibull(shape=2, scale=1); (0.886227, 1.0, 1.32934) (result1 <- phfit.3mom(0.886227, 1.0, 1.32934)) ## Three moment matching ## Moments of Weibull(shape=2, scale=1); (0.886227, 1.0, 1.32934) (result2 <- phfit.3mom(0.886227, 1.0, 1.32934, method="Bobbio05")) ## mean ph.mean(result1) ph.mean(result2) ## variance ph.var(result1) ph.var(result2) ## up to 5 moments ph.moment(5, result1) ph.moment(5, result2)
Estimates PH parameters from density function.
phfit.density( ph, f, deformula = deformula.zeroinf, weight.zero = 1e-12, weight.reltol = 1e-08, start.divisions = 8, max.iter = 12, ... )phfit.density( ph, f, deformula = deformula.zeroinf, weight.zero = 1e-12, weight.reltol = 1e-08, start.divisions = 8, max.iter = 12, ... )
ph |
An object of R6 class. The estimation algorithm is selected depending on this class. |
f |
A function object for a density function. |
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). |
weight.zero |
A absolute value which is regarded as zero in numerical integration. |
weight.reltol |
A value for precision of numerical integration. |
start.divisions |
A value for starting value of divisions in deformula. |
max.iter |
A value for the maximum number of iterations to increase divisions in deformula. |
... |
Options for EM steps, which is also used to send the arguments to density function. |
Returns a list with components, which is an object of S3 class phfit.result;
model |
an object for estimated PH class. |
llf |
a value of the maximum log-likelihood (a negative value of the cross entropy). |
df |
a value of degrees of freedom of the model. |
KL |
a value of Kullback-Leibler divergence. |
iter |
the number of iterations. |
convergence |
a logical value for the convergence of estimation algorithm. |
ctime |
computation time (user time). |
data |
an object for data class |
aerror |
a value of absolute error for llf at the last step of algorithm. |
rerror |
a value of relative error for llf at the last step of algorithm. |
options |
a list of options. |
call |
the matched call. |
Any of density function can be applied to the argument f, where
f should be defined f <- function(x, ...).
The first argument of f should be an integral parameter.
The other parameters are set in the argument ... of phfit.density.
The truncated density function can also be used directly.
#################### ##### 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)#################### ##### 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)
Estimates PH parameters from grouped data.
phfit.group(ph, counts, breaks, intervals, instants, ...)phfit.group(ph, counts, breaks, intervals, instants, ...)
ph |
An object of R6 class. The estimation algorithm is selected depending on this class. |
counts |
A vector of the number of points in intervals. |
breaks |
A vector for a sequence of points of boundaries of intervals.
This is equivalent to |
intervals |
A vector of time lengths for intervals.
This is equivalent to |
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 |
... |
Further options for EM steps. |
Returns a list with components, which is an object of S3 class phfit.result;
model |
an object for estimated PH class. |
llf |
a value of the maximum log-likelihood. |
df |
a value of degrees of freedom of the model. |
aic |
a value of Akaike information criterion. |
iter |
the number of iterations. |
convergence |
a logical value for the convergence of estimation algorithm. |
ctime |
computation time (user time). |
data |
an object for data class |
aerror |
a value of absolute error for llf at the last step of algorithm. |
rerror |
a value of relative error for llf at the last step of algorithm. |
options |
a list of options used in the fitting. |
call |
the matched call. |
In this method, we can handle truncated data using NA and Inf;
phfit.group(ph=cf1(5), counts=c(countsdata, NA), breaks=c(breakdata, +Inf))
NA means missing of count data at the corresponding interval, and Inf is allowed to put
the last of breaks or intervals which represents a special interval [the last break point,infinity).
## 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)## 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)
Estimates PH parameters from point data.
phfit.point(ph, x, weights, ...)phfit.point(ph, x, weights, ...)
ph |
An object of R6 class for PH. The estimation algorithm is selected depending on this class. |
x |
A vector for point data. |
weights |
A vector of weights for points. |
... |
Further options for fitting methods. |
Returns a list with components, which is an object of S3 class phfit.result;
model |
an object for estimated PH class. |
llf |
a value of the maximum log-likelihood. |
df |
a value of degrees of freedom of the model. |
aic |
a value of Akaike information criterion. |
iter |
the number of iterations. |
convergence |
a logical value for the convergence of estimation algorithm. |
ctime |
computation time (user time). |
data |
an object for data class |
aerror |
a value of absolute error for llf at the last step of algorithm. |
rerror |
a value of relative error for llf at the last step of algorithm. |
options |
a list of options used for fitting. |
call |
the matched call. |
## 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)## 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)
Compute the cumulative distribution function (c.d.f.) for a given PH distribution
pphase(q, ph, lower.tail = TRUE, log.p = FALSE, ...)pphase(q, ph, lower.tail = TRUE, log.p = FALSE, ...)
q |
A numeric vector of quantiles. |
ph |
An instance of PH distribution. |
lower.tail |
logical; if TRUE, probabilities are P(X <= x), otherwise P(X > x). |
log.p |
logical; if TRUE, probabilities p are returned as log(p). |
... |
Others |
A vector of densities
## create a PH with specific parameters (phdist <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## c.d.f. for 0, 0.1, ..., 1 pphase(q=seq(0, 1, 0.1), ph=phdist)## create a PH with specific parameters (phdist <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## c.d.f. for 0, 0.1, ..., 1 pphase(q=seq(0, 1, 0.1), ph=phdist)
Generate a sample from a given PH distribution.
rphase(n, ph, ...)rphase(n, ph, ...)
n |
An integer of the number of samples. |
ph |
An instance of PH distribution. |
... |
Others |
A vector of samples.
## create a PH with specific parameters (phdist <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## generate 10 samples rphase(n=10, ph=phdist)## create a PH with specific parameters (phdist <- ph(alpha=c(1,0,0), Q=rbind(c(-4,2,0),c(2,-5,1),c(1,0,-1)), xi=c(2,2,0))) ## generate 10 samples rphase(n=10, ph=phdist)