fit_dsmm {dsmmR} | R Documentation |
Estimation of a drifting semi-Markov chain
Description
Estimation of a drifting semi-Markov chain, given one sequence of states. This estimation can be parametric or non-parametric and is available for the three types of drifting semi-Markov models.
Usage
fit_dsmm(
sequence,
degree,
f_is_drifting,
p_is_drifting,
states = NULL,
initial_dist = "unif",
estimation = "nonparametric",
f_dist = NULL
)
Arguments
sequence |
Character vector that represents a sequence of states
from the state space |
degree |
Positive integer that represents the polynomial degree
|
f_is_drifting |
Logical. Specifies if |
p_is_drifting |
Logical. Specifies if |
states |
Character vector that represents the state space
|
initial_dist |
Optional. Character that represents the method to estimate the initial distribution.
|
estimation |
Optional. Character. Represents whether the estimation will be nonparametric or parametric.
|
f_dist |
Optional. It can be defined in two ways:
|
Details
This function estimates a drifting semi-Markov model in the parametric and non-parametric case. The parametric estimation can be achieved by the following steps:
We obtain the non-parametric estimation of the sojourn time distributions.
We estimate the parameters for the distributions defined in the attribute
f_dist
through the probabilities that were obtained in step 1.
Three different models are possible for to be estimated for each case.
A normalization technique is used in order to correct estimation errors
from small sequences.
We will use the exponentials to distinguish between
the drifting semi-Markov kernel
and the
semi-Markov kernels
used in
Model 1, Model 2, Model 3.
More about the theory of drifting semi-Markov models in dsmmR.
Non-parametric Estimation
Model 1
When the transition matrix of the embedded Markov chain
and
the conditional sojourn time distribution
are both drifting,
the drifting semi-Markov kernel can be estimated as:
, where
is the maximum
sojourn time that was observed in the sequence and
are
polynomials with degree
(see dsmmR).
The semi-Markov kernels
,
are estimated through Least Squares Estimation (LSE) and are obtained
after solving the following system,
,
and
:
where the matrices are written as:
and we use the following indicator functions:
-
, if at
the previous state is
,
otherwise.
-
, if at
the previous state is
with sojourn time
and next state
,
otherwise.
In order to obtain the estimations of
and
, we use the following formulas:
Model 2
In this case, is drifting and
is not drifting. Therefore,
the estimated drifting semi-Markov kernel will be given by:
, where
is the maximum
sojourn time that was observed in the sequence and
are
polynomials with degree
(see dsmmR).
In order to obtain the estimators
and
,
we use the estimated semi-Markov kernels
from Model 1.
Since
is drifting, we define the estimation
the same way as we did in Model 1.
In total, we have the following estimations,
:
Thus, the estimated semi-Markov kernels for Model 2,
, can be written with
regards to the estimated semi-Markov kernels of Model 1,
, as in the following:
Model 3
In this case, is drifting and
is not drifting. Therefore,
the estimated drifting semi-Markov kernel will be given by:
, where
is the maximum
sojourn time that was observed in the sequence and
are
polynomials with degree
(see dsmmR).
In order to obtain the estimators
and
,
we use the estimated semi-Markov kernels
estimated semi-Markov kernels
from Model 1. Since
is drifting,
we define the estimation
the same way as we did in
Model 1. In total, we have the following estimations,
:
Thus, the estimated semi-Markov kernels for Model 3,
, can be written with
regards to the estimated semi-Markov kernels of Model 1,
, as in the following:
Parametric Estimation
In this package, the parametric estimation of the sojourn time distributions
defined in the attribute f_dist
is achieved as follows:
-
We obtain the non-parametric LSE of the sojourn time distributions
.
-
We estimate the parameters for the distributions defined in
f_dist
through the probabilities of, estimated in previously in
.
The available distributions for the modelling of the conditional sojourn
times of the drifting semi-Markov model, defined from
the argument f_dist
, have their parameters estimated through the
following formulas:
Geometric
:
, where
. We estimate the probability of success
as such:
Poisson
:
, where
. We estimate
as such:
Negative binomial
:
, where
.
is the Gamma function,
is the probability of success and
is the parameter describing the number of successful trials, or the dispersion parameter (the shape parameter of the gamma mixing distribution). We estimate them as such:
Discrete Weibull of type 1
:
, where
,
is the first parameter with
and
the second parameter. We estimate them as such:
Note that we require
for estimating
.
Uniform
:
where
, for
a positive integer. We use a numerical method to obtain an estimator for
in this case.
Value
Returns an object of S3 class (dsmm_fit_nonparametric, dsmm)
or (dsmm_fit_parametric, dsmm)
. It has the following attributes:
-
dist
: List. Contains 2 or 3 arrays.If
estimation = "nonparametric"
we have 2 arrays:-
p_drift
orp_notdrift
, corresponding to whether the definedtransition matrix is drifting or not.
-
f_drift
orf_notdrift
, corresponding to whether the definedsojourn time distribution is drifting or not.
-
If
estimation = "parametric"
we have 3 arrays:-
p_drift
orp_notdrift
, corresponding to whether the definedtransition matrix is drifting or not.
-
f_drift_parametric
orf_notdrift_parametric
, corresponding to whether the definedsojourn time distribution is drifting or not.
-
f_drift_parameters
orf_notdrift_parameters
, which are the definedsojourn time distribution parameters, depending on whether
is drifting or not.
-
-
emc
: Character vector that contains the embedded Markov chainof the original sequence. It is this attribute of the object that describes the size of the model
. Last state is also included, for a total length of
, but it is not used for any calculation.
-
soj_times
: Numerical vector that contains the sojourn times spent for each state inemc
before the jump to the next state. Last state is also included, for a total length of, but it is not used for any calculation.
-
initial_dist
: Numerical vector that contains an estimation for the initial distribution of the realized states insequence
. It always has values betweenand
.
-
states
: Character vector. Passing down from the arguments. It contains the realized states given in the argumentsequence
. -
s
: Positive integer that contains the length of the character vector given in the attributestates
, which is equal to.
-
degree
: Positive integer. Passing down from the arguments. It contains the polynomial degreeconsidered for the drifting of the model.
-
k_max
: Numerical value that contains the maximum sojourn time, which is the maximum value insoj_times
, excluding the last state. -
model_size
: Positive integer that contains the size of the drifting semi-Markov model, which is equal to the length of the embedded Markov chain
, minus the last state. It has a value of
length(emc) - 1
, foremc
as defined above. -
f_is_drifting
: Logical. Passing down from the arguments. Specifies ifis drifting or not.
-
p_is_drifting
: Logical. Passing down from the arguments. Specifies ifis drifting or not.
-
Model
: Character. Possible values:-
"Model_1"
: Bothand
are drifting.
-
"Model_2"
:is drifting and
is not drifting.
-
"Model_3"
:is drifting and
is not drifting.
-
-
estimation
: Character. Specifies whether parametric or nonparametric estimation was used. -
A_i
: Numerical Matrix. Represents the polynomialswith degree
that were used for solving the system
. Used for the methods defined for the object. Not printed when viewing the object.
-
J_i
: Numerical Array. Represents the estimated semi-Markov kernels of the first modelthat were obtained after solving the system
. Not printed when viewing the object.
References
V. S. Barbu, N. Limnios. (2008). semi-Markov Chains and Hidden semi-Markov Models Toward Applications - Their Use in Reliability and DNA Analysis. New York: Lecture Notes in Statistics, vol. 191, Springer.
Vergne, N. (2008). Drifting Markov models with Polynomial Drift and Applications to DNA Sequences. Statistical Applications in Genetics Molecular Biology 7 (1).
Barbu V. S., Vergne, N. (2019). Reliability and survival analysis for Drifting Markov models: modelling and estimation. Methodology and Computing in Applied Probability, 21(4), 1407-1429.
T. Nakagawa and S. Osaki. (1975). The discrete Weibull distribution. IEEE Transactions on Reliability, R-24, 300-301.
See Also
For the theoretical background of drifting semi-Markov models: dsmmR.
For sequence simulation: simulate.dsmm and create_sequence.
For drifting semi-Markov model specification: parametric_dsmm, nonparametric_dsmm
For the retrieval of the drifting semi-Markov kernel: get_kernel.
Examples
# Create a random sequence
sequence <- create_sequence("DNA", len = 2000, seed = 1)
## Alternatively, we could obtain a sequence as follows:
## > data("lambda", package = "dsmmR")
## > sequence <- c(lambda)
states <- sort(unique(sequence))
degree <- 3
# ===========================================================================
# Nonparametric Estimation.
# Fitting a random sequence under distributions of unknown shape.
# ===========================================================================
# ---------------------------------------------------------------------------
# Both p and f are drifting - Model 1.
# ---------------------------------------------------------------------------
obj_model_1 <- fit_dsmm(sequence = sequence,
degree = degree,
f_is_drifting = TRUE,
p_is_drifting = TRUE,
states = states,
initial_dist = "freq",
estimation = "nonparametric", # default value
f_dist = NULL # default value
)
cat(paste0("We fitted a sequence with ", obj_model_1$Model, ",\n",
"model size: n = ", obj_model_1$model_size, ",\n",
"length of state space: s = ", obj_model_1$s, ",\n",
"maximum sojourn time: k_max = ", obj_model_1$k_max, " and\n",
"polynomial (drifting) Degree: d = ", obj_model_1$degree, ".\n"))
# Get the drifting p and f arrays.
p_drift <- obj_model_1$dist$p_drift
f_drift <- obj_model_1$dist$f_drift
cat(paste0("Dimension of p_drift: (s, s, d + 1) = (",
paste(dim(p_drift), collapse = ", "), ").\n",
"Dimension of f_drift: (s, s, k_max, d + 1) = (",
paste(dim(f_drift), collapse = ", "), ").\n"))
# We can even check the embedded Markov chain and the sojourn times
# directly from the returned object, if we wish to do so.
# This is achieved through the `base::rle()` function, used on `sequence`.
model_emc <- obj_model_1$emc
model_sojourn_times <- obj_model_1$soj_times
# ---------------------------------------------------------------------------
# Fitting the sequence when p is drifting and f is not drifting - Model 2.
# ---------------------------------------------------------------------------
obj_model_2 <- fit_dsmm(sequence = sequence,
degree = degree,
f_is_drifting = FALSE,
p_is_drifting = TRUE)
cat(paste0("We fitted a sequence with ", obj_model_2$Model, ".\n"))
# Get the drifting p and non-drifting f arrays.
p_drift_2 <- obj_model_2$dist$p_drift
f_notdrift <- obj_model_2$dist$f_notdrift
all.equal.numeric(p_drift, p_drift_2) # p is the same as in Model 1.
cat(paste0("Dimension of f_notdrift: (s, s, k_max) = (",
paste(dim(f_notdrift), collapse = ", "), ").\n"))
# ---------------------------------------------------------------------------
# Fitting the sequence when f is drifting and p is not drifting - Model 3.
# ---------------------------------------------------------------------------
obj_model_3 <- fit_dsmm(sequence = sequence,
degree = degree,
f_is_drifting = TRUE,
p_is_drifting = FALSE)
cat(paste0("We fitted a sequence with ", obj_model_3$Model, ".\n"))
# Get the drifting f and non-drifting p arrays.
p_notdrift <- obj_model_3$dist$p_notdrift
f_drift_3 <- obj_model_3$dist$f_drift
all.equal.numeric(f_drift, f_drift_3) # f is the same as in Model 1.
cat(paste0("Dimension of f_notdrift: (s, s) = (",
paste(dim(p_notdrift), collapse = ", "), ").\n"))
# ===========================================================================
# Parametric Estimation
# Fitting a random sequence under distributions of known shape.
# ===========================================================================
### Comments
### 1. For the parametric estimation it is recommended to use a common set
### of distributions while only the parameters (of the sojourn times)
### are drifting. This results in (generally) higher accuracy.
### 2. This process is similar to that used in `dsmm_parametric()`.
s <- length(states)
# Getting the distributions for the states.
# Rows correspond to previous state `u`.
# Columns correspond to next state `v`.
f_dist_1 <- matrix(c(NA, "unif", "dweibull", "nbinom",
"pois", NA, "pois", "dweibull",
"geom", "pois", NA, "geom",
"dweibull", 'geom', "pois", NA),
nrow = s, ncol = s, byrow = TRUE)
f_dist <- array(f_dist_1, dim = c(s, s, degree + 1))
dim(f_dist)
# ---------------------------------------------------------------------------
# Both p and f are drifting - Model 1.
# ---------------------------------------------------------------------------
obj_fit_parametric <- fit_dsmm(sequence = sequence,
degree = degree,
f_is_drifting = TRUE,
p_is_drifting = TRUE,
states = states,
initial_dist = 'unif',
estimation = 'parametric',
f_dist = f_dist)
cat("The class of `obj_fit_parametric` is : (",
paste0(class(obj_fit_parametric), collapse = ', '), ").\n")
# Estimated parameters.
f_params <- obj_fit_parametric$dist$f_drift_parameters
# The drifting sojourn time distribution parameters.
f_0 <- f_params[,,,1]
f_1.3 <- f_params[,,,2]
f_2.3 <- f_params[,,,3]
f_1 <- f_params[,,,4]
params <- paste0('q = ', round(f_params["c", "t", 1, ], 3),
', beta = ', round(f_params["c", "t", 2, ], 3))
f_names <- c("f_0", paste0("f_", 1:(degree-1), "/", degree), "f_1")
all_names <- paste(f_names, ":", params)
cat("The drifting of the parameters for passing from \n",
"`u` = 'c' to `v` = 't' under a discrete Weibull distribution is:",
"\n", all_names[1], "\n", all_names[2],
"\n", all_names[3], "\n", all_names[4])
# ---------------------------------------------------------------------------
# f is not drifting, only p is drifting - Model 2.
# ---------------------------------------------------------------------------
obj_fit_parametric_2 <- fit_dsmm(sequence = sequence,
degree = degree,
f_is_drifting = FALSE,
p_is_drifting = TRUE,
initial_dist = 'unif',
estimation = 'parametric',
f_dist = f_dist_1)
cat("The class of `obj_fit_parametric_2` is : (",
paste0(class(obj_fit_parametric_2), collapse = ', '), ").\n")
# Estimated parameters.
f_params_2 <- obj_fit_parametric_2$dist$f_notdrift_parameters
params_2 <- paste0('q = ', round(f_params_2["c", "t", 1], 3),
', beta = ', round(f_params_2["c", "t", 2], 3))
cat("Not-drifting parameters for passing from ",
"`u` = 'c' to `v` = 't' \n under a discrete Weibull distribution are:\n",
paste("f :", params_2))
# ===========================================================================
# `simulate()` and `get_kernel()` can be used for the two objects,
# `dsmm_fit_nonparametric` and `dsmm_fit_parametric`.
# ===========================================================================
sim_seq_nonparametric <- simulate(obj_model_1, nsim = 10)
str(sim_seq_nonparametric)
kernel_drift_parametric <- get_kernel(obj_fit_parametric, klim = 10)
str(kernel_drift_parametric)