fit_dsmm: Estimation of a drifting semi-Markov chain
In dsmmR: Estimation and Simulation of Drifting Semi-Markov Models
fit_dsmm 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 E.
degree
Positive integer that represents the polynomial degree
d for the drifting semi-Markov model.
f_is_drifting
Logical. Specifies if f is drifting or not.
p_is_drifting
Logical. Specifies if p is drifting or not.
states
Character vector that represents the state space
E, with length equal to s = |E|. Default value is set
equal to the sorted, unique states present in the given sequence.
initial_dist
Optional. Character that represents the method to
estimate the initial distribution.
-
"unif" : The initial distribution of each state is
equal to 1/s (default value).
-
"freq" : The initial distribution of each state is
equal to the frequency that it has in the sequence.
estimation
Optional. Character. Represents whether the
estimation will be nonparametric or parametric.
-
"nonparametric" : The estimation will be
non-parametric (default value).
-
"parametric" : The estimation will be
parametric.
f_dist
Optional. It can be defined in two ways:
If estimation = "nonparametric", it is
equal to NULL (default value).
If estimation = "parametric", it is a
character array that specifies the distributions of the sojourn
times, for every state transition.
The list of possible values is:
["unif", "geom", "pois", "dweibull", "nbinom", NA].
It can be defined in two ways:
-
If f is not drifting, it has dimensions of
s \times s.
-
If f is drifting, it has dimensions of
s \times s \times (d+1)
(see more in Details, Parametric Estimation.)
It is defined similarly to the attribute f_dist
in dsmm_parametric.
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 (1), (2), (3) to distinguish between
the drifting semi-Markov kernel \widehat{q}_{\frac{t}{n}} and the
semi-Markov kernels \widehat{q}_\frac{i}{d} 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 p of the embedded Markov chain
(J_{t})_{t\in \{0,\dots,n\}} and
the conditional sojourn time distribution f are both drifting,
the drifting semi-Markov kernel can be estimated as:
\widehat{q}_{\frac{t}{n}}^{\ (1)}(u,v,l) =
\sum_{i = 0}^{d}A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),
\forall t \in \{0,\dots,n\}, \forall u,v\in E,
\forall l \in \{1,\dots, k_{max} \} , where k_{max} is the maximum
sojourn time that was observed in the sequence and
A_i, i = 0, \dots, d are d + 1 polynomials with degree
d (see dsmmR).
The semi-Markov kernels
\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l), i = 0, \dots, d,
are estimated through Least Squares Estimation (LSE) and are obtained
after solving the following system, \forall t \in \{0, \dots, n\},
\forall u, v \in E and \forall l \in \{1, \dots, k_{max}\}:
MJ = P,
where the matrices are written as:
-
M = (M_{ij})_{i,j \in \{0, \dots, d\} } =
\left(\sum_{t=1}^{n}1_{u}(t)A_{i}(t)A_{j}(t)\right)_{
i,j \in \{0, \dots, d\}}
-
J = (J_i)_{i \in \{0, \dots, d\} } =
\left(\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)\right)_{
i \in \{0, \dots, d\}}
-
P=(P_i)_{i\in \{0, \dots, d\} }=
\left(\sum_{t=1}^{n}1_{uvl}(t)A_{i}(t)\right)_{
i \in \{0, \dots, d\}}
and we use the following indicator functions:
-
1_{u}(t) = 1_{ \{J_{t-1} = u \} } = 1,
if at t the previous state is u,
0 otherwise.
-
1_{uvl}(t) = 1_{ \{J_{t-1} = u, J_{t} = v, X_{t} = l \} } = 1,
if at t the previous state is u
with sojourn time l and next state v, 0 otherwise.
In order to obtain the estimations of \widehat{p}_{\frac{i}{d}}(u,v)
and \widehat{f}_{\frac{i}{d}}(u,v,l), we use the following formulas:
\widehat{p}_{\frac{i}{d}}(u,v) =
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),
\widehat{f}_{\frac{i}{d}}(u,v,l) =
\frac{\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Model 2
In this case, p is drifting and f is not drifting. Therefore,
the estimated drifting semi-Markov kernel will be given by:
\widehat{q}_{\frac{t}{n}}^{\ (2)}(u,v,l) =
\sum_{i=0}^{d} A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l),
\forall t \in \{0,\dots,n\}, \forall u,v\in E,
\forall l\in \{1,\dots, k_{max} \}, where k_{max} is the maximum
sojourn time that was observed in the sequence and
A_i, i = 0, \dots, d are d + 1 polynomials with degree
d (see dsmmR).
In order to obtain the estimators \widehat{p} and \widehat{f},
we use the estimated semi-Markov kernels
\widehat{q}_{\frac{i}{d}}^{\ (1)} from Model 1.
Since p is drifting, we define the estimation \widehat{p}
the same way as we did in Model 1.
In total, we have the following estimations,
\forall u,v \in E, \forall l \in \{1,\dots, k_{max} \}:
\widehat{p}_{\frac{i}{d}}(u,v) =
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),
\widehat{f}(u,v,l) =
\frac{\sum_{i=0}^{d}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{
\sum_{i=0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Thus, the estimated semi-Markov kernels for Model 2,
\widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l) =
\widehat{p}_{\frac{i}{d}}(u,v)\widehat{f}(u,v,l), can be written with
regards to the estimated semi-Markov kernels of Model 1,
\widehat{q}_{\frac{i}{d}}^{\ (1)}, as in the following:
\widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l) = \frac{
(\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))
(\sum_{i = 0}^{d}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))}{
\sum_{i = 0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Model 3
In this case, f is drifting and p is not drifting. Therefore,
the estimated drifting semi-Markov kernel will be given by:
\widehat{q}_{\frac{t}{n}}^{\ (3)}(u,v,l) =
\sum_{i=0}^{d} A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l),
\forall t \in \{0,\dots,n\}, \forall u,v\in E,
\forall l\in \{1,\dots, k_{max} \}, where k_{max} is the maximum
sojourn time that was observed in the sequence and
A_i, i = 0, \dots, d are d + 1 polynomials with degree
d (see dsmmR).
In order to obtain the estimators \widehat{p} and \widehat{f},
we use the estimated semi-Markov kernels
estimated semi-Markov kernels \widehat{q}_{\frac{i}{d}}^{\ (1)}
from Model 1. Since f is drifting,
we define the estimation \widehat{f} the same way as we did in
Model 1. In total, we have the following estimations,
\forall u,v \in E, \forall l \in \{1,\dots, k_{max} \}:
\widehat{p}\ (u,v) =
\frac{\sum_{i=0}^{d}\sum_{l=1}^{k_{max}}
\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{d+1},
\widehat{f}_{\frac{i}{d}}(u,v,l) =
\frac{\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Thus, the estimated semi-Markov kernels for Model 3,
\widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l) =
\widehat{p}\ (u,v)\widehat{f}_{\frac{i}{d}}(u,v,l), can be written with
regards to the estimated semi-Markov kernels of Model 1,
\widehat{q}_{\frac{i}{d}}^{\ (1)}, as in the following:
\widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l) = \frac{
\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)
\sum_{i=0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}
{(d+1)\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
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
f.
-
We estimate the parameters for the distributions defined in
f_dist through the probabilities of f, estimated
in previously in 1.
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 (p):
f(x) = p (1-p)^{x-1}, where x = 1, 2, \dots,k_{max}.
We estimate the probability of success \widehat{p}
as such:
\widehat{p} = \frac{1}{E(X)}
Poisson (\lambda):
f(x) = \frac{\lambda^{x-1} exp(-\lambda)}{(x-1)!}, where
x = 1, 2,\dots,k_{max}. We estimate \widehat{\lambda} > 0
as such:
\widehat{\lambda} = E(X)
Negative binomial (\alpha, p):
f(x)=\frac{\Gamma(x + \alpha - 1)}{\Gamma(\alpha)(x-1)!}
p^{\alpha}(1-p)^{x-1}, where
x = 1, 2,\ldots,k_{max}.
\Gamma is the Gamma function,
p is the probability of success and
\alpha \in (0, +\infty) 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:
\widehat{p} = \frac{E(X)}{Var(X)},
\widehat{\alpha} = E(X)\frac{\widehat{p}}{1 - \widehat{p}} =
\frac{E(X)^2}{Var(X) - E(X)}.
Discrete Weibull of type 1 (q, \beta):
f(x)=q^{(x-1)^{\beta}}-q^{x^{\beta}}, where
x= 1, 2, \dots,k_{max}, q is the first parameter with
0 < q < 1 and \beta \in (0, +\infty)
the second parameter. We estimate them as such:
\widehat{q} = 1 - f(1),
\widehat{\beta} = \frac{\sum_{i = 2}^{k_{max}}
\log_{i}(\log_{\widehat{q}}(\sum_{j = 1}^{i}f(j)))}{k_{max} - 1}.
Note that we require k_{max} \geq 2 for estimating
\widehat{\beta}.
Uniform (n):
f(x) = 1/n where x = 1, 2, \dots, n, for n a
positive integer. We use a numerical method to obtain an estimator
for \widehat{n} 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 or p_notdrift, corresponding to whether the
defined p transition matrix is drifting or not.
-
f_drift or f_notdrift, corresponding to whether the
defined f sojourn time distribution is drifting or not.
If estimation = "parametric" we have 3 arrays:
-
p_drift or p_notdrift, corresponding to whether the
defined p transition matrix is drifting or not.
-
f_drift_parametric or f_notdrift_parametric,
corresponding to whether the
defined f sojourn time distribution is drifting or not.
-
f_drift_parameters or f_notdrift_parameters,
which are the defined f sojourn time distribution parameters,
depending on whether f is drifting or not.
-
emc : Character vector that contains the
embedded Markov chain (J_{t})_{t\in \{0,\dots,n\}}
of the original sequence.
It is this attribute of the object that describes the size of the model
n. Last state is also included, for a total length of n+1,
but it is not used for any calculation.
-
soj_times : Numerical vector that contains the sojourn times
spent for each state in emc before the jump to the next state.
Last state is also included, for a total length of n+1,
but it is not used for any calculation.
-
initial_dist : Numerical vector that contains an estimation
for the initial distribution of the realized states in sequence.
It always has values between 0 and 1.
-
states : Character vector. Passing down from the arguments.
It contains the realized states given in the argument sequence.
-
s : Positive integer that contains the length of the character
vector given in the attribute states, which is equal to s = |E|.
-
degree : Positive integer. Passing down from the arguments.
It contains the polynomial degree
d considered for the drifting of the model.
-
k_max : Numerical value that contains the maximum sojourn
time, which is the maximum value in soj_times, excluding the last
state.
-
model_size : Positive integer that contains the size of the
drifting semi-Markov model n, which is equal to the length of the
embedded Markov chain (J_{t})_{t\in \{0,\dots,n\}},
minus the last state.
It has a value of length(emc) - 1, for emc as defined above.
-
f_is_drifting : Logical. Passing down from the arguments.
Specifies if f is drifting or not.
-
p_is_drifting : Logical. Passing down from the arguments.
Specifies if p is drifting or not.
-
Model : Character. Possible values:
-
"Model_1" : Both p and f are drifting.
-
"Model_2" : p is drifting and
f is not drifting.
-
"Model_3" : f is drifting and
p is not drifting.
-
estimation : Character. Specifies whether parametric or
nonparametric estimation was used.
-
A_i : Numerical Matrix. Represents the polynomials
A_i(t) with degree d that were used for solving
the system MJ = P. 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 model
(\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))_{i\in\{0,\dots,d\}}
that were obtained after solving the system MJ = P.
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)
dsmmR documentation built on Sept. 14, 2024, 9:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| fit_dsmm | 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_distthrough 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 (1), (2), (3) to distinguish between
the drifting semi-Markov kernel \widehat{q}_{\frac{t}{n}} and the
semi-Markov kernels \widehat{q}_\frac{i}{d} 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 p of the embedded Markov chain
(J_{t})_{t\in \{0,\dots,n\}} and
the conditional sojourn time distribution f are both drifting,
the drifting semi-Markov kernel can be estimated as:
\widehat{q}_{\frac{t}{n}}^{\ (1)}(u,v,l) =
\sum_{i = 0}^{d}A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),
\forall t \in \{0,\dots,n\}, \forall u,v\in E,
\forall l \in \{1,\dots, k_{max} \} , where k_{max} is the maximum
sojourn time that was observed in the sequence and
A_i, i = 0, \dots, d are d + 1 polynomials with degree
d (see dsmmR).
The semi-Markov kernels
\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l), i = 0, \dots, d,
are estimated through Least Squares Estimation (LSE) and are obtained
after solving the following system, \forall t \in \{0, \dots, n\},
\forall u, v \in E and \forall l \in \{1, \dots, k_{max}\}:
MJ = P,
where the matrices are written as:
-
M = (M_{ij})_{i,j \in \{0, \dots, d\} } = \left(\sum_{t=1}^{n}1_{u}(t)A_{i}(t)A_{j}(t)\right)_{ i,j \in \{0, \dots, d\}} -
J = (J_i)_{i \in \{0, \dots, d\} } = \left(\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)\right)_{ i \in \{0, \dots, d\}} -
P=(P_i)_{i\in \{0, \dots, d\} }= \left(\sum_{t=1}^{n}1_{uvl}(t)A_{i}(t)\right)_{ i \in \{0, \dots, d\}}
and we use the following indicator functions:
-
1_{u}(t) = 1_{ \{J_{t-1} = u \} } = 1, if attthe previous state isu,0otherwise. -
1_{uvl}(t) = 1_{ \{J_{t-1} = u, J_{t} = v, X_{t} = l \} } = 1, if attthe previous state isuwith sojourn timeland next statev,0otherwise.
In order to obtain the estimations of \widehat{p}_{\frac{i}{d}}(u,v)
and \widehat{f}_{\frac{i}{d}}(u,v,l), we use the following formulas:
\widehat{p}_{\frac{i}{d}}(u,v) =
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),
\widehat{f}_{\frac{i}{d}}(u,v,l) =
\frac{\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Model 2
In this case, p is drifting and f is not drifting. Therefore,
the estimated drifting semi-Markov kernel will be given by:
\widehat{q}_{\frac{t}{n}}^{\ (2)}(u,v,l) =
\sum_{i=0}^{d} A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l),
\forall t \in \{0,\dots,n\}, \forall u,v\in E,
\forall l\in \{1,\dots, k_{max} \}, where k_{max} is the maximum
sojourn time that was observed in the sequence and
A_i, i = 0, \dots, d are d + 1 polynomials with degree
d (see dsmmR).
In order to obtain the estimators \widehat{p} and \widehat{f},
we use the estimated semi-Markov kernels
\widehat{q}_{\frac{i}{d}}^{\ (1)} from Model 1.
Since p is drifting, we define the estimation \widehat{p}
the same way as we did in Model 1.
In total, we have the following estimations,
\forall u,v \in E, \forall l \in \{1,\dots, k_{max} \}:
\widehat{p}_{\frac{i}{d}}(u,v) =
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l),
\widehat{f}(u,v,l) =
\frac{\sum_{i=0}^{d}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{
\sum_{i=0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Thus, the estimated semi-Markov kernels for Model 2,
\widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l) =
\widehat{p}_{\frac{i}{d}}(u,v)\widehat{f}(u,v,l), can be written with
regards to the estimated semi-Markov kernels of Model 1,
\widehat{q}_{\frac{i}{d}}^{\ (1)}, as in the following:
\widehat{q}_{\frac{i}{d}}^{\ (2)}(u,v,l) = \frac{
(\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))
(\sum_{i = 0}^{d}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))}{
\sum_{i = 0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Model 3
In this case, f is drifting and p is not drifting. Therefore,
the estimated drifting semi-Markov kernel will be given by:
\widehat{q}_{\frac{t}{n}}^{\ (3)}(u,v,l) =
\sum_{i=0}^{d} A_{i}(t)\ \widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l),
\forall t \in \{0,\dots,n\}, \forall u,v\in E,
\forall l\in \{1,\dots, k_{max} \}, where k_{max} is the maximum
sojourn time that was observed in the sequence and
A_i, i = 0, \dots, d are d + 1 polynomials with degree
d (see dsmmR).
In order to obtain the estimators \widehat{p} and \widehat{f},
we use the estimated semi-Markov kernels
estimated semi-Markov kernels \widehat{q}_{\frac{i}{d}}^{\ (1)}
from Model 1. Since f is drifting,
we define the estimation \widehat{f} the same way as we did in
Model 1. In total, we have the following estimations,
\forall u,v \in E, \forall l \in \{1,\dots, k_{max} \}:
\widehat{p}\ (u,v) =
\frac{\sum_{i=0}^{d}\sum_{l=1}^{k_{max}}
\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{d+1},
\widehat{f}_{\frac{i}{d}}(u,v,l) =
\frac{\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}{
\sum_{l = 1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
Thus, the estimated semi-Markov kernels for Model 3,
\widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l) =
\widehat{p}\ (u,v)\widehat{f}_{\frac{i}{d}}(u,v,l), can be written with
regards to the estimated semi-Markov kernels of Model 1,
\widehat{q}_{\frac{i}{d}}^{\ (1)}, as in the following:
\widehat{q}_{\frac{i}{d}}^{\ (3)}(u,v,l) = \frac{
\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)
\sum_{i=0}^{d}\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}
{(d+1)\sum_{l=1}^{k_{max}}\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l)}.
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
f. -
We estimate the parameters for the distributions defined in
f_distthrough the probabilities off, estimated in previously in1.
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
(p):f(x) = p (1-p)^{x-1}, wherex = 1, 2, \dots,k_{max}. We estimate the probability of success\widehat{p}as such:\widehat{p} = \frac{1}{E(X)}Poisson
(\lambda):f(x) = \frac{\lambda^{x-1} exp(-\lambda)}{(x-1)!}, wherex = 1, 2,\dots,k_{max}. We estimate\widehat{\lambda} > 0as such:\widehat{\lambda} = E(X)Negative binomial
(\alpha, p):f(x)=\frac{\Gamma(x + \alpha - 1)}{\Gamma(\alpha)(x-1)!} p^{\alpha}(1-p)^{x-1}, wherex = 1, 2,\ldots,k_{max}.\Gammais the Gamma function,pis the probability of success and\alpha \in (0, +\infty)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:\widehat{p} = \frac{E(X)}{Var(X)},\widehat{\alpha} = E(X)\frac{\widehat{p}}{1 - \widehat{p}} = \frac{E(X)^2}{Var(X) - E(X)}.Discrete Weibull of type 1
(q, \beta):f(x)=q^{(x-1)^{\beta}}-q^{x^{\beta}}, wherex= 1, 2, \dots,k_{max},qis the first parameter with0 < q < 1and\beta \in (0, +\infty)the second parameter. We estimate them as such:\widehat{q} = 1 - f(1),\widehat{\beta} = \frac{\sum_{i = 2}^{k_{max}} \log_{i}(\log_{\widehat{q}}(\sum_{j = 1}^{i}f(j)))}{k_{max} - 1}.Note that we require
k_{max} \geq 2for estimating\widehat{\beta}.Uniform
(n):f(x) = 1/nwherex = 1, 2, \dots, n, forna positive integer. We use a numerical method to obtain an estimator for\widehat{n}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_driftorp_notdrift, corresponding to whether the definedptransition matrix is drifting or not. -
f_driftorf_notdrift, corresponding to whether the definedfsojourn time distribution is drifting or not.
-
If
estimation = "parametric"we have 3 arrays:-
p_driftorp_notdrift, corresponding to whether the definedptransition matrix is drifting or not. -
f_drift_parametricorf_notdrift_parametric, corresponding to whether the definedfsojourn time distribution is drifting or not. -
f_drift_parametersorf_notdrift_parameters, which are the definedfsojourn time distribution parameters, depending on whetherfis drifting or not.
-
-
emc: Character vector that contains the embedded Markov chain(J_{t})_{t\in \{0,\dots,n\}}of the original sequence. It is this attribute of the object that describes the size of the modeln. Last state is also included, for a total length ofn+1, but it is not used for any calculation. -
soj_times: Numerical vector that contains the sojourn times spent for each state inemcbefore the jump to the next state. Last state is also included, for a total length ofn+1, 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 between0and1. -
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 tos = |E|. -
degree: Positive integer. Passing down from the arguments. It contains the polynomial degreedconsidered 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 modeln, which is equal to the length of the embedded Markov chain(J_{t})_{t\in \{0,\dots,n\}}, minus the last state. It has a value oflength(emc) - 1, foremcas defined above. -
f_is_drifting: Logical. Passing down from the arguments. Specifies iffis drifting or not. -
p_is_drifting: Logical. Passing down from the arguments. Specifies ifpis drifting or not. -
Model: Character. Possible values:-
"Model_1": Bothpandfare drifting. -
"Model_2":pis drifting andfis not drifting. -
"Model_3":fis drifting andpis not drifting.
-
-
estimation: Character. Specifies whether parametric or nonparametric estimation was used. -
A_i: Numerical Matrix. Represents the polynomialsA_i(t)with degreedthat were used for solving the systemMJ = P. 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 model(\widehat{q}_{\frac{i}{d}}^{\ (1)}(u,v,l))_{i\in\{0,\dots,d\}}that were obtained after solving the systemMJ = P. 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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.