Discretised version of coliform counts in sea-water samples

Description

Discretised version of counts of faecal coliform bacteria in sea water samples collected at seven locations near Sydney NSW, Australia. There were four “controls”: Longreef, Bondi East, Port Hacking “50”, and Port Hacking “100” and three “outfalls”: Bondi Offshore, Malabar Offshore and North Head Offshore. At each location measurements were made at four depths: 0, 20, 40, and 60 meters. A large fraction of the counts are missing values.

Usage

SydColDisc

Format

A data frame with 5432 observations on the following 6 variables.

y

A factor consisting of a discretisation of counts of faecal coliform count bacteria in sea water samples. The original measures were obtained by a repeated dilution process.

The data were discretised using the cut() function with breaks given by c(0,1,5,25,200,Inf) and labels equal to c("lo","mlo","m","mhi","hi").

locn

a factor with levels Longreef, Bondi East, Port Hacking 50, Port Hacking 100, Bondi Offshore, Malabar Offshore and North Head Offshore.

depth

a factor with levels 0 (0 metres), 20 (20 metres), 40 (40 metres), 60 (60 metres).

ma.com

A factor with levels no and yes, indicating whether the Malabar sewage outfall had been commissioned.

nh.com

A factor with levels no and yes, indicating whether the North Head sewage outfall had been commissioned.

bo.com

A factor with levels no and yes, indicating whether the Bondi Offshore sewage outfall had been commissioned.

Details

The observations corresponding to each location-depth combination constitute a (discrete valued) time series. The sampling interval is ostensibly 1 week; distinct time series are ostensibly synchronous. The measurements were made over a 194 week period. Due to exigencies of weather, the unreliabitity of boats and other factors the collection times were actually highly irregular and have been rounded to the neares week. Often no sample was obtained at a given site within a week of the putative collection time, in which the observed count is given as a missing value. In fact over 75% of the counts are missing. See Turner et al. (1998) for more detail.

Modelling

The hidden Markov models applied in the paper Turner et al. (1998) and in the paper Turner (2008) used a numeric version of the response in this data set. The numeric response was essentially a square root transformation of the original data, and the resulting values were modelled in terms of a Poisson distribution. See the references for details.

Source

The original data were kindly supplied by Geoff Coade, of the New South Wales Environment Protection Authority (Australia)

References

T. Rolf Turner, Murray A. Cameron, and Peter J. Thomson. Hidden Markov chains in generalized linear models. Canadian J. Statist. 26 (1998) 107 – 125.

Rolf Turner. Direct maximization of the likelihood of a hidden Markov model. Comp. Statist. Data Anal. 52 (2008) 4147–4160.

Anova for hmm.discnp models

Description

Performs a likelihood ratio test to compare two discrete non-parametric hidden Markov models.

Usage

## S3 method for class 'hmm.discnp'
anova(object, ...)

Arguments

Value

A list with entries

This list has an attribute “details” which is a vector consisting of the first and second log likelihoods and the associated numbers of parameters, in order of these numbers of parameters. (See Warning.)

Warning

Hidden Markov models can be numerically delicate and the fitting algorithm can converge to a local maximum of the likelihood surface which is not the global maximum. Thus it is entirely possible for the log likelihood of the model with the greater number of parameters to be smaller than that of the model with the lesser number of parameters.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

hmm()

Examples

xxx  <- with(SydColDisc,split(y,f=list(locn,depth)))
fit1 <- hmm(xxx,K=1,itmax=10)
fit2 <- hmm(xxx,K=2,itmax=10)
anova(fit1,fit2)

Simulated monocyte counts and psychosis symptoms.

Description

Discretised values of monocyte counts, and ratings of level of psychosis simulated from a model fitted to a data set consisting of observations made on a number of patients from the Northern District Health Board system. The real data must be kept confidential due to ethics constraints.

Usage

    data("ccprSim")

Format

The object ccprSim is a list of length 1258. Each entry of this list is to be considered to correspond to an individual subject. The entries consist of matrices having two columns named cellCount and psychosisRating. The number of rows of these matrices varies from entry to entry of the list (i.e. from subject to subject).

Most of the entries of these matrices are NA. The entries are temporally ordered and correspond to the number of weeks from the start of observation. Observations in the real data set were made only when the patient in question visted a physician and so weeks in which no visit was made resulted in an “observation” of NA. The object ccprSim was simulated in such a way as to imitate this characteristic. The fraction of missing observations in each variate (i.e. cellCount and psychosisRating is roughly commensurate with the corresponding fractions in the real data.

The values in the first column of each matrix (the cellCount column) consist of integers from 1 to 5 and are to be interpreted as indicators of cell counts in units of 10^9 cells per litre, discretised according to the following scale:

• 0.0 \leq c \leq 0.3 \leftrightarrow 1

• 0.3 < c \leq 0.5 \leftrightarrow 2

• 0.5 < c \leq 0.7 \leftrightarrow 3

• 0.7 < c \leq 1.0 \leftrightarrow 4

• 1.0 < c \leq 2.0 \leftrightarrow 5

where c represents “count”.

The values in the second column of each matrix (the psychosisRating column consist of integers from 0 to 4 and are to be interpreted as indicators of a physician's assessment of the level of pschosis of the patient. A value of 0 corresponds to “no symptoms”; a value of 4 corresponds to “severe”.

The question of essential interest in respect of the real data was “Is there any association between the cell count values and the psychosis ratings?” More specifically it was “Can the level of psychosis be predicted from the cell counts?”

Source

The real data, on the basis of which these data were simulated, were supplied by Dr. Jonathan Williams, Northern District Health Board.

Examples

## Not run:  # Takes too long.
    fit <- hmm(ccprSim,K=2,indep=FALSE,itmax=5,verbose=TRUE)

## End(Not run)

Convert Rho between forms.

Description

Converts a matrix specification of the emission probabilities (in which the probabilities are simply the entries of the matrix) to a data frame specification (in which the probabilities are a logistic-style function of the parameters) or vice versa.

Usage

cnvrtRho(Rho)

Arguments

Details

The data frame specification of Rho allows for predictor variables x. If Rho is of the data frame form, and is designed to allow for predictor variables, then it will have more than three columns and cannot be converted to the matrix form. In such cases cnvrtRho will throw an error.

Value

A data frame if the argument Rho is a matrix, or a matrix if the argument Rho is a data frame.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

hmm()

Examples

Yval <- LETTERS[1:10]
Tpm  <- matrix(c(0.75,0.25,0.25,0.75),ncol=2,byrow=TRUE)
Rho  <- cbind(c(rep(1,5),rep(0,5)),c(rep(0,5),rep(1,5)))/5
rownames(Rho) <- Yval
newRho <- cnvrtRho(Rho)
oldRho <- cnvrtRho(newRho)

Fitted values of a discrete non-parametric hidden Markov model.

Description

Calculates the fitted values of a discrete non-parametric hidden Markov model. If the data are numeric these are the conditional expectations of the observations, given the entire observation sequence (and the estimated parameters of the model). If the data are categorical (whence “expectations” make no sense) the “fitted values” are taken to be the probabilities of each of the possible values of the observations, at each time point.

Usage

## S3 method for class 'hmm.discnp'
fitted(object, warn=TRUE, drop=TRUE, ...)

Arguments

Details

The observation sequence(s) must be present in object (which will be the case if object was returned by hmm() and if the argument keep.y was set to TRUE). If it is not present an error is thrown.

However, if such an error is thrown, do not despair! You do not have to start from scratch when fitting your model with keep.y==TRUE. If fit is your fitted model that you obtained without setting keep.y==TRUE, then you can just re-fit the model using fit as the starting values:

      fit2 <- hmm(<whatever>,par0=fit,keep.y=TRUE)
  

This will of course converge instantaneously. You could also do:

      fit2 <- update(fit,data=<whatever>,keep.y=TRUE)
  

Value

If the observations consist of a single sequence and if drop is TRUE then the returned value consists of a single object (matrix, list of two matrices, or 3-dimensional array, depending on circumstances; see below). Otherwise the returned value is a list of such objects, one for each observation sequence.

If the observations are numeric then the object corresponding to each observation sequence is a matrix. If the model is univariate (see hmm()) then matrix has a single column constituting the sequence of fitted values corresponding to the observations in the given sequence. The number of rows is the number of observations and the entry in row t is the fitted value (conditional expection) corresponding to the observation made at time t. If the model is bivariate (either independent or dependent) then the matrix has two columns corresponding respectively to the two variables in the bivariate model.

If the observations are categorical then the nature of the object returned changes substantially depending on whether the data are univariate, bivariate independent or bivariate dependent. (See hmm().

In the unvariate case the object corresponding to each sequence is a matrix, the number of rows of which is the number of observations and the number of columns of which is the number of unique possible values of the observations. The entry of this matrix in row t and column j is the conditional probability that an emission, at time t, is equal to u_i where u_1, \dots, u_m are the unique possible values.

In the bivariate independent case the object is a list of two matrices, each of which is of the same nature as that produced in the univariate case, corresponding respectively to the first and second of the two variables.

In the bivariate dependent case the object is a 3-dimensional array of dimension m_1 \times m_2 \times n where m_1 is the number of unique possible values of the first variable, m_2 is the number of unique possible values of the second variable, and n is the number of observations. The (i,j,t)-th entry of this array is the conditional probability that an emission, at time t, is equal to (u_i,v_j) where the u_i are the unique possible values of the first variable and the v_j are the unique possible values of the second variable.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

sp() link{predict.hmm.discnp}()

Examples

P <- matrix(c(0.7,0.3,0.1,0.9),2,2,byrow=TRUE)
R <- matrix(c(0.5,0,0.1,0.1,0.3,
              0.1,0.1,0,0.3,0.5),5,2)
set.seed(42)
lll <- sample(250:350,20,TRUE)
y   <- rhmm(ylengths=lll,nsim=1,drop=TRUE,tpm=P,Rho=R)
fit <- hmm(y,K=2,verb=TRUE,keep.y=TRUE,itmax=10)
fv  <- fitted(fit)

Fit a hidden Markov model to discrete data.

Description

Effects a maximum likelihood fit of a hidden Markov model to discrete data where the observations come from one of a number of finite discrete distributions, depending on the (hidden) state of the Markov chain. These distributions (the “emission probabilities”) are specified non-parametrically. The observations may be univariate, independent bivariate, or dependent bivariate. By default this function uses the EM algorithm. In the univariate setting it may alternatively use a “brute force” method.

Usage

hmm(y, yval=NULL, par0=NULL, K=NULL, rand.start=NULL,
    method=c("EM","bf","LM","SD"), hglmethod=c("fortran","oraw","raw"),
    optimiser=c("nlm","optim"), optimMethod=NULL, stationary=cis,
    mixture=FALSE, cis=TRUE, indep=NULL, tolerance=1e-4, digits=NULL,
    verbose=FALSE, itmax=200, crit=c("PCLL","L2","Linf","ABSGRD"),
    X=NULL,keep.y=FALSE, keep.X=keep.y,
    addIntercept=TRUE, lmc=10, hessian=FALSE,...)

Arguments

Details

• Univariate case: In the univariate case the emission probabilities are specified by means of a data frame Rho. The first column of Rho, named "y", is a factor consisting of the possible values of the emissions, repeated K times (where K is the number of states). The second column, named states, is a factor consisting of integer values 1, 2, ..., K. Each of these values is repeated m times where m is the length of yval. Further columns of Rho are numeric and consist of coefficients of the linear predictor of the probabilities of the various values of y. If X is NULL then Rho has only one further column named Intercept.

  If X is not NULL then the Intercept column is present only if addIntercept is TRUE. There as many (other, in addition to the possible Intercept column) numeric columns as there are columns in X or in the matrices in the list X. The names of these columns are taken to be the column names of X or the first entry of X if such column names are present. Otherwise the names default to V1, V2 ....

  The probabilities of the emissions taking on their various possible values are given by

  \Pr(Y = y_i | \boldsymbol{x}, \textrm{state}=S) = \ell_i/\sum_{j=1}^m \ell_j

  where \ell_j is the j\textrm{th} entry of \boldsymbol{\beta}^{\top}\boldsymbol{x} and where in turn \boldsymbol{x} is the vector of predictors and \boldsymbol{\beta} is the coefficient vector in the linear predicator that corresponds to y_i and the hidden state S. For identifiability the vectors \boldsymbol{\beta} corresponding to the first value of Y (the first level of Rho$y) are set equal to the zero vector for all values of the state S.

  Note that the Rho component of the starting values par0 may be specified as a matrix of probabilities, with rows corresponding to possible values of the observations and columns corresponding to states. That is the Rho component of par0 may be provided in the form \textrm{Rho} = [\rho_{ij}] where \rho_{ij} = \Pr(Y = y_i | S = j). This is permissible as long as X is NULL and may be found to be more convenient and intuitive. If the starting value for Rho is provided in matrix form it is (silently) converted internally into the data frame form, by the (undocumented) function cnvrtRho().

  When argument X is not NULL, it is difficult to specify a “reasonable” value for the Rho component of par0. One might try to specify par0$Rho in the data frame form. The question of how to specify the columns of par0$Rho corresponding to the auxiliary predictors (columns of X or of the entries of X) is a thorny one.

  It is permissible in these circumstances to specify par0$Rho as a matrix of probabilities, just as one would do if X were NULL. In this setting the (undocumented) function checkStartVal() converts the matrix of probabilities to data frame form and then appends columns, all of whose entries are 0, corresponding to the auxiliary predictors. When par0 is unspecified, the (undocumented) function init.all() performs similar construction to accommodate a non-NULL value of X. Whether the resulting starting value for Rho makes any real sense, is questionable. However little else can be done.

• Independent bivariate case: the emission probabilities are specified by a list of two matrices. In this setting \Pr(Y_1,Y_2) = (y_{i1},y_{i2}) | S = j) = \rho^{(1)}_{i_1,j} \rho^{(2)}_{i_2,j} where R^{(k)} = [\rho^{(k)}_{ij}] (k = 1,2) are the two emission probability matrices.

• Dependent bivariate case: the emission probabilities are specified by a three dimensional array. In this setting \Pr((Y_1,Y_2) = (y_{i1},y_{i2}) | S = j) = \rho_{i_1,i_2,j} where R = [\rho_{ijk}] is the emission probability array.

The hard work of calculating the recursive probabilities used to fit the model is done by a Fortran subroutine recurse (actually coded in Ratfor) which is dynamically loaded. In the univariate case, when X is provided, the estimation of the “linear predictor” vectors \boldsymbol{\beta} is handled by the function multinom() from the nnet package. Note that this is a “Recommended” package and is thereby automatically available (i.e. does not have to be installed).

Value

A list with components:

Thanks

A massive nest of bugs was eliminated in the transition from version 3.0-8 to version 3.0-9. These bugs arose in the context of using auxiliary predictor variables (argument X). The handling of such auxiliary predictors was completely messed up. I am grateful to Leah Walker for pointing out the problem to me.

Warnings

The ordering of the (hidden) states can be arbitrary. What the estimation procedure decides to call “state 1” may not be what you think of as being state number 1. The ordering of the states will be affected by the starting values used.

Some experiences with using the "ABSGRD" stopping criterion indicate that it may be problematic in the context of discrete non-parametric distributions. For example a value of 1854.955 was returned after 200 LM steps in one (non-convergent, of course!) attempt at fitting a model. The stopping criterion "PCLL" in this example took the “reasonable” value of 0.03193748 when iterations ceased.

Notes — Various

This function used to have an argument newstyle, a logical scalar (defaulting to TRUE) indicating whether (in the univariate setting) the emission probabilities should be represented in “logistic” form. (See Details, Univariate case:, above.) Now the emission probabilities are always represented in the “logistic” form. The component Rho of the starting parameter values par0 may still be supplied as a matrix of probabilities (with columns summing to 1), but this component is converted (internally, silently) to the logistic form.

The object returned by this function also has (in the univariate setting), in addition to the component Rho, a component Rho.matrix giving the emission probabilities in the more readily interpretable matrix-of-probabilities form. (See Value above.)

The package used to require the argument y to be a matrix in the case of multiple observed sequences. If the series were of unequal length the user was expected to pad them out with NAs to equalize the lengths.

The old matrix format for multiple observation sequences was permitted for a while (and the matrix was internally changed into a list) but this is no longer allowed. If y is indeed given as a matrix then this corresponds to a single observation sequence and it must have one (univariate setting) or two (bivariate setting) columns which constitute the observations of the respective variates.

If K=1 then tpm, ispd, converged, and nstep are all set equal to NA in the list returned by this function.

The estimate of ispd in the non-stationary setting is inevitably very poor, unless the number of sequences of observations (the length of the list y) is very large. We have in effect “less than one” relevant observation for each such sequence.

The returned values of tpm and Rho (or the entries of Rho when Rho is a list) have dimension names. These are formed from the argument yval if this is supplied, otherwise from the sorted unique values of the observations in y. Likewise the returned value of ispd is a named vector, the names being the same as the row (and column) names of tpm.

If method is equal to "EM" there may be a decrease (!!!) in the log likelihood at some EM step. This is “theoretically impossible” but can occur in practice due to an intricacy in the way that the EM algorithm treats ispd when stationary is TRUE. It turns out to be effectively impossible to maximise the expected log likelihood unless the term in that quantity corresponding to ispd is ignored (whence it is ignored). Ignoring this term is “asymptotically negligible” but can have the unfortunate effect of occasionally leading to a decrease in the log likelihood.

If such a decrease is detected, then the algorithm terminates and issues a message to the effect that the decrease occurred. The message suggests that another method be used and that perhaps the results from the penultimate EM step (which are returned by this function) be used as starting values.

It seems to me that it should be the case that such a decrease in the log likelihood can occur only if stationary is TRUE. However I have encountered instances in which a decrease occurred when stationary was FALSE. I have yet to figure out/track down what is going on here.

Note on method

If the method is "EM" it is actually possible for the log likelihood to decrease at some EM step. This is “impossible in an ideal world” but can happen to the fact the EM algorithm, as implemented in this package at least, cannot maximise the expected log likelihood if the component corresponding to the initial state probability distribution is taken into consideration. This component should ideally be maximised subject to the constraint that t(P)%*%ispd = ispd, but this constraint seems to effectively impossible to impose. Lagrangian multipliers don't cut it. Hence the summand in question is ignored at the M-step. This usually works alright since the summand is asymptotically negligible, but things can sometimes go wrong. If such a decrease occurs, an error is thrown.

In previous versions of this package, instead of throwing an error the hmm() function would automatically switch to either the "bf" or the "LM" method, depending whether a matrix X of auxiliary predictors is supplied, starting from the penultimate parameter estimates produced by the EM algorithm. However this appears not to be a good idea; those “penultimate estimates” appear not to be good starting values for the other methods. Hence an error is now thrown and the user is explicitly instructed to invoke a different method, “starting from scratch”.

Fitted Coefficients of the Predictors

It is of course of interest to understand the meaning of the coefficients that are fitted to the predictors in the model. If X is supplied then the number of predictors is (as a rule) one (for the intercept) plus the number of columns in each entry of X. We say “as a rule” because, e.g., the entries of X could each have an “intercept” column, or the addIntercept argument could be FALSE. If X is not supplied there is only one predictor, named Intercept.

The interpretation of these predictor coefficients is a bit subtle. To get an idea of what it's all about, consider the output from example 4. (See Examples). The fitted coefficients in question are to be found in columns 3 and onward of the component Rho of the object returned by hmm(). In the context of example 4, this object is fit.wap. (The suffix wap stands for “with auxiliary predictors”.)

  fit.wap$Rho
       y state Intercept     ma.com      nh.com     bo.com
  1   lo     1 1.3810463  0.4527982 -3.27161353 -1.9563915
  2  mlo     1 0.1255631 -1.1402546 -1.37713744  0.5946980
  3    m     1 0.7356526  0.1523734 -2.70841817 -0.1794645
  4  mhi     1 0.8479798 -0.2438988 -1.12544989 -0.9650320
  5   hi     1 0.0000000  0.0000000  0.00000000  0.0000000
  6   lo     2 3.9439410 -0.8355306 -0.77702276  1.4963631
  7  mlo     2 2.6189880 -1.9373885 -0.09190623  0.8316870
  8    m     2 2.1457317 -1.7276183  0.19524655 -0.3249485
  9  mhi     2 1.8834139 -1.3760011 -0.59806309  1.2828365
  10  hi     2 0.0000000  0.0000000  0.00000000  0.0000000

If you multiply the matrix consisting of the predictor coefficients (columns 3 to 6 of Rho in this instance) times a vector of predictors you get, for each state, the “exponential form” of the probabilities (“pre-probabilities”) for each of the possible y-values, given the vector of predictors.

E.g. set x <- c(1,1,0,0). This vector picks up the intercept and indicates that the Malabar outfall has been commissioned, the North Head outfall has not been commissioned, and the Bondi Offshore outfall has not been commissioned.

Now set:

    pp1 <- (as.matrix(fit.wap$Rho)[,3:6]%*%x)[1:5]
    pp2 <- (as.matrix(fit.wap$Rho)[,3:6]%*%x)[6:10]

Note that pp1 consists of “exponential probabilities” corresponding to state 1, and pp2 consists of “exponential probabilities” corresponding to state 2. To convert the foregoing pre-probabilities to the actual probabilities of the y-values, we apply the — undocumented — function expForm2p():

    p1 <- expForm2p(pp1)
    p2 <- expForm2p(pp2)

The value of p1 is

[1] 0.52674539 0.03051387 0.20456767 0.15400019 0.08417288

and that of p2 is

[1] 0.78428283 0.06926632 0.05322204 0.05819340 0.03503541

Note that p1 and p2 each sum to 1, as they should/must do. This says, e.g., that when the system is in state 2, and Malabar has been commissioned but North Head and Bondi Offshore have not, the (estimated) probability that y is "mhi" (medium-high) is 0.05819340.

It may be of some interest to test the hypothesis that the predictors have any actual predictive power at all:

    fit.nap <- hmm(xxx,yval=Yval,K=2,verb=TRUE)
    # "nap" <--> no aux. preds

There is a bit of a problem here, in that the likelihood decreases at EM step 65. (See the warning message.)

We can check on this problem by refitting using method="LM".

    fit.nap.lm <- hmm(xxx,yval=Yval,par0=fit.nap,method="LM",verb=TRUE)

Doing so produces only a small improvement in the log likelihood (from -1821.425 to -1820.314), so we really could have ignored the problem. We can now do anova(fit.wap,fit.nap) which gives

    $stat
    [1] 153.5491

    $df
    [1] 24

    $pvalue
    [1] 7.237102e-21

Thus the p-value is effectively zero, saying that in this instance the auxiliary predictors appear to have a “significant” impact on the fit.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

References

Rabiner, L. R., "A tutorial on hidden Markov models and selected applications in speech recognition," Proc. IEEE vol. 77, pp. 257 – 286, 1989.

Zucchini, W. and Guttorp, P., "A hidden Markov model for space-time precipitation," Water Resources Research vol. 27, pp. 1917-1923, 1991.

MacDonald, I. L., and Zucchini, W., "Hidden Markov and Other Models for Discrete-valued Time Series", Chapman & Hall, London, 1997.

Liu, Limin, "Hidden Markov Models for Precipitation in a Region of Atlantic Canada", Master's Report, University of New Brunswick, 1997.

See Also

rhmm(), mps(), viterbi()

Examples

# TO DO: Create one or more bivariate examples.
#
# The value of itmax in the following examples is so much
# too small as to be risible.  This is just to speed up the
# R CMD check process.
# 1.
Yval <- LETTERS[1:10]
Tpm  <- matrix(c(0.75,0.25,0.25,0.75),ncol=2,byrow=TRUE)
Rho  <- cbind(c(rep(1,5),rep(0,5)),c(rep(0,5),rep(1,5)))/5
rownames(Rho) <- Yval
set.seed(42)
xxx  <- rhmm(ylengths=rep(1000,5),nsim=1,tpm=Tpm,Rho=Rho,yval=Yval,drop=TRUE)
fit  <- hmm(xxx,par0=list(tpm=Tpm,Rho=Rho),itmax=10)
print(fit$Rho) # A data frame
print(cnvrtRho(fit$Rho)) # A matrix of probabilities
                         # whose columns sum to 1.

# 2.
# See the help for logLikHmm() for how to generate y.num.
## Not run: 
   fit.num     <- hmm(y.num,K=2,verb=TRUE,itmax=10)
   fit.num.mix <- hmm(y.num,K=2,verb=TRUE,mixture=TRUE,itmax=10)
   print(fit.num[c("tpm","Rho")])

## End(Not run)
# Note that states 1 and 2 get swapped.

# 3.
xxx <- with(SydColDisc,split(y,f=list(locn,depth)))
Yval <- c("lo","mlo","m","mhi","hi")
# Two states: above and below the thermocline.
fitSydCol <- hmm(xxx,yval=Yval,K=2,verb=TRUE,itmax=10)

# 4.
X <- split(SydColDisc[,c("ma.com","nh.com","bo.com")],
           f=with(SydColDisc,list(locn,depth)))
X <- lapply(X,function(x){
                 as.matrix(as.data.frame(lapply(x,as.numeric)))-1})
fit.wap <- hmm(xxx,yval=Yval,K=2,X=X,verb=TRUE,itmax=10)
# wap <--> with auxiliary predictors.

# 5.
## Not run:  # Takes too long.
fitlm <- hmm(xxx,yval=Yval,K=2,method="LM",verb=TRUE)
fitem <- hmm(xxx,yval=Yval,K=2,verb=TRUE)
# Algorithm terminates due to a decrease in the log likelihood
# at EM step 64.
newfitlm <- hmm(xxx,yval=Yval,par0=fitem,method="LM",verb=TRUE)
# The log likelihood improves from -1900.988 to -1820.314

## End(Not run)

# 6.
fitLesCount <- hmm(lesionCount,K=2,itmax=10) # Two states: relapse and remission.

Internal hmm.discnp functions.

Description

Internal hmm.discnp functions.

Usage

## S3 method for class 'multipleHmmDataSets'
x[i]
check.yval(yval, Rho, type, warn=TRUE)
checkConstColumns(y,prednames)
checkStartVal(par0,K,indep,yval,rand.start,mixture,prednames)
checkyXoK(y, X)
derivf(theta, K)
derivp(theta, K)
derivpi(ispd, tpm, npar, dp)
ell(phi, G)
expForm2p(x)
ffun(Dat, Rho, type)
forgethgl(fy, y, ymiss, tpm, ispd, d1pi, d2pi, npar, d1p,
          d2p, m, d1f, d2f)
orgethgl(fy, y, ymiss, tpm, xispd, d1pi, d2pi, npar, d1p,
          d2p, m, d1f, d2f)
rgethgl(fy, y, ymiss, tpm, xispd, d1pi, d2pi, npar, d1p,
          d2p, m, d1f, d2f)
get.gl(theta, K, y)
get.hgl(theta, K, y, hglmethod)
get.l(theta, K, y)
getIspd(pars, K)
getRho(pars, K, rhovals, stationary, prednames)
getTpm(pars, K, stationary)
hmmBD(y, par0, K, stationary,
      mixture, cis, tolerance, digits, verbose,
      itmax, crit, bicm)
hmmBI(y, par0, K, stationary,
      mixture, cis, tolerance, digits, verbose,
      itmax, crit, bicm)
hmmNumOpt(Dat, par0, stationary, verbose, itmax, bicm, rhovals, npar, optimiser,
          optimMethod, hessian=FALSE, ...)
hmmLM(y, par0, itmax=200, crit, lmc=10, tolerance,
      bicm, rhovals, hglmethod, digits=NULL, verbose=FALSE) 
hmmSD(y, par0, itmax=200, crit, tolerance,
      bicm, rhovals, hglmethod, digits=NULL, verbose=FALSE) 
hmmUV(y, par0, K, method,
      hglmethod, optimiser, optimMethod, stationary,
      mixture, cis, tolerance, digits, verbose,
      itmax, crit, bicm, X, addIntercept, lmc, hessian, ...)
init.all(K, rand.start, mixture, indep, yval, prednames)
lmstep(theta, K, y, lmc, hglmethod)
lse(z)
msRho(Rho0, G)
makeDat(y, X, addIntercept)
ordinal(k)
ordinalsuffix(k)
paramExtract(Rho)
p2expForm(x)
phi2rho(phi, ijk)
predictEngineHmmDiscnp(stateProbs,Rho,numb,drop)
recurse(fy, tpm, ispd, lns)
reparam(object, expForm=TRUE, stationary=NULL)
revise.ispd(tpm=NULL, gamma=NULL, lns=NULL, cis=TRUE)
revise.rho(Dat, gamma, type)
revise.tpm(xi, mixture)
rho2phi(Rho)
simference(object, expForm=TRUE, seed=NULL,nsim=100,verbose=TRUE)
steepest(K, y, theta)
tidyList(y, rp=c("response","predictor"),addIntercept=NULL, yval=NULL)

Details

These functions are auxiliary and are not intended to be called by the user.

Canadian hydrological data sets.

Description

Five data sets obtained from the “HYDAT” database, Environment and Climate Change Canada's database of historical hydrometric data. The data were obtained using the tidyhydat package. The data have been trimmed so that there are no gaps in the observation dates and are presented in “raw” form and in discretised form as deciles of the residuals (difference between raw values and the daily mean over years).

Usage

    data("linLandFlows")
    data("ftLiardFlows")
    data("portMannFlows")
    data("portMannSedLoads")
    data("portMannSedCon")

Format

Data frames with observations on the following 3 variables.

Date

Dates on which observations were made.

Value

Numeric vector of observation values.

mean

The mean over years of Value.

resid

The difference Value - mean.

deciles

A factor with levels d1, ..., d10, which are the deciles of the variable resid

Details

The variable mean was calculated as follows:

        yday <- as.POSIXlt(X$Date)$yday
        mn   <- tapply(X$Value,yday,mean,na.rm=TRUE)
        mean <- mn[as.character(yday)]
   

where X is the data set being processed.

The data set linLandFlows originally consisted of 2008 observations; there were 1980 observations after “trimming”. The data set ftLiardFlows originally consisted of 22364 observations; there were 11932 observations after “trimming”. The data set portMannFlows originally consisted of 6455 observations; there were 3653 observations after “trimming”. The data set portMannSedLoads consists of 2771 observations; no observations were trimmed. The data set portMannSedCon consists of 4597 observations; no observations were trimmed.

The units of the “Flows” variables are cubic metres per second (m^3/s); the units of “portMannSedLoads” are tonnes; the units of “portMannSedCon” are milligrams per litre (mg/l).

The “linLandFlows” data were obtained at the Lindberg Landing hydrometric station on the Liard River in the Northwest Territories of Canada. The “ftLiardFlows” data were obtained at the Fort Liard hydrometric station on the Liard River in the Northwest Territories of Canada. The “portMann” data were obtained at the hydrometric station located at the Port Mann pumping station on the Fraser River in the Province of British Columbia in Canada.

Source

Environment and Climate Change Canada's database “HYDAT”, a database of historical hydrometric data. The data were obtained vis the tidyhydat package, which is available from “CRAN”, https://cran.r-project.org

Examples

fit <- hmm(linLandFlows$deciles,K=4,itmax=10)

Multiple sclerosis lesion counts for three patients.

Description

Lesion counts for three multiple sclerosis patients. The counts were obtained by magnetic resonance imaging, and were observed at monthly intervals.

Usage

  lesionCount

Format

A list with three components each component being the sequence of counts for a given patient and consisting of a vector with non-negative integer entries.

Modelling

The hidden Markov models applied to these data by Albert et al. and by MacKay and Petkau were much more complex and elaborate than those fitted in the examples in this package. See the references for details.

Source

The data were originally studied by Albert et al., (1994). They are were also analyzed by Altman and Petkau (2005). The data were kindly provided by Prof. Altman.

References

Albert, P. S., McFarland, H. F., Smith, M. E., and Frank, J. A. Time series for modelling counts from a relapsing-remitting disease: application to modelling disease activity in multiple sclerosis. Statistics in Medicine 13 (1994) 453–466.

Altman, Rachel MacKay, and Petkau, A. John. Application of hidden Markov models to multiple sclerosis lesion count data. Statistics in Medicine 24 (2005) 2335–2344.

Log likelihood of a hidden Markov model

Description

Calculate the log likelihood of a hidden Markov model with discrete non-parametric observation distributions.

Usage

logLikHmm(y, model=NULL, tpm=NULL, ispd=NULL, Rho=NULL, X=NULL,
          addIntercept=NULL, warn=TRUE)

Arguments

Details

If y is not provided the function simply returns the log.like component of model (which could be NULL if model was not produced by hmm().

The observation values (the entries of the vector or matrix y, or of the vectors or matrices which constitute the entries of y if y is a list) must be consistent with the appropriate dimension names of Rho or of its entries when Rho is a list. More specifically, if Rho has dimension names (or its entries have dimension names) then the observation values must all be found as entries of the appropriate dimension name vector. If a vector of dimension names is NULL then the corresponding dimension must be equal to the number of unique observations of the appropriate variate. integers between 1 and nrow(Rho).

Value

The loglikehood of y given the parameter values specified in par.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

References

See hmm() for references.

See Also

hmm(), pr(), sp()

Examples

# TO DO:  One or more bivariate examples.
P <- matrix(c(0.7,0.3,0.1,0.9),2,2,byrow=TRUE)
R <- matrix(c(0.5,0,0.1,0.1,0.3,
              0.1,0.1,0,0.3,0.5),5,2)
set.seed(42)
lll   <- sample(250:350,20,TRUE)
set.seed(909)
y.num <- rhmm(ylengths=lll,nsim=1,tpm=P,Rho=R,drop=TRUE)
set.seed(909)
y.let <- rhmm(ylengths=lll,nsim=1,tpm=P,Rho=R,yval=letters[1:5],drop=TRUE)
row.names(R) <- 1:5
ll1   <- logLikHmm(y.num,tpm=P,Rho=R)
row.names(R) <- letters[1:5]
ll2   <- logLikHmm(y.let,tpm=P,Rho=R)
ll3   <- logLikHmm(y.let,tpm=P,Rho=R,ispd=c(0.5,0.5))
fit   <- hmm(y.num,K=2,itmax=10)
ll4   <- logLikHmm(y.num,fit) # Use the fitted rather than the "true" parameters.

Insert missing values.

Description

Insert missing values into data simulated by rhmm.

Usage

misstify(y, nafrac, fep = NULL)

Arguments

Value

An object with a structure similar to that of y, containing the same data as y but with some of these data having been replaced by missing values (NA). In particular, if y is of class "multipleHmmDataSets" then so is the returned value.

Note that rhmm() calls upon misstify() to effect the replacement of a certain fraction of the simulated observations by missing values. If rhmm() is applied to a fitted model, then by default, this “certain fraction” is determined, using nafracCalc(), from the data set to which the model was fitted.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

rhmm() nafracCalc()

Examples

P <- matrix(c(0.7,0.3,0.1,0.9),2,2,byrow=TRUE)
R <- matrix(c(0.5,0,0.1,0.1,0.3,
              0.1,0.1,0,0.3,0.5),5,2)
set.seed(42)
lll   <- sample(250:350,20,TRUE)
y1    <- rhmm(ylengths=lll,nsim=1,tpm=P,Rho=R)
y1m   <- misstify(y1,nafrac=0.5,fep=list(TRUE))
y2    <- rhmm(ylengths=lll,nsim=5,tpm=P,Rho=R)
set.seed(127)
y2m   <- misstify(y2,nafrac=0.5,fep=list(TRUE))
nafracCalc(y2m) # A list all of whose entries are close to 0.5.
set.seed(127)
y2ma  <- lapply(y2,misstify,nafrac=0.5,fep=list(TRUE))
## Not run: 
    nafracCalc(y2ma) # Throws an error.

## End(Not run)
sapply(y2ma,nafracCalc) # Effectively the same as nafracCalc(y2m).

Most probable states.

Description

Calculates the most probable hidden state underlying each observation.

Usage

mps(y, model = NULL, tpm, Rho, ispd=NULL, warn=TRUE)

Arguments

Details

For each t the maximum value of \gamma_t(i), i.e. of the (estimated) probability that the state at time t is equal to i, is calculated, and the value of the state with the corresponding index is returned.

Value

If y is a single observation sequence, then the value is a vector of corresponding most probable states.

If y is a list of replicate sequences, then the value is a list, the j-th entry of which constitutes the vector of most probable states underlying the j-th replicate sequence.

If y is of class "multipleHmmDataSets" then the value returned is a list of lists of the sort described above.

Warning

The sequence of most probable states as calculated by this function will not in general be the most probable sequence of states. It may not even be a possible sequence of states. This function looks at the state probabilities separately for each time t, and not at the states in their sequential context.

To obtain the most probable sequence of states use viterbi().

Author(s)

Rolf Turner r.turner@auckland.ac.nz

References

Rabiner, L. R., "A tutorial on hidden Markov models and selected applications in speech recognition," Proc. IEEE vol. 77, pp. 257 – 286, 1989.

See Also

hmm(), rhmm(), viterbi()

Examples

## Not run: 
P <- matrix(c(0.7,0.3,0.1,0.9),2,2,byrow=TRUE)
rownames(P) <- 1:2
R <- matrix(c(0.5,0,0.1,0.1,0.3,
              0.1,0.1,0,0.3,0.5),5,2)
set.seed(42)
lll   <- sample(250:350,20,TRUE)
set.seed(909)
y.num <- rhmm(ylengths=lll,nsim=1,tpm=P,Rho=R,drop=TRUE)
fit.num <- hmm(y.num,K=2,verb=TRUE)
s.1 <- mps(y.num,fit.num)
s.2 <- mps(y.num,tpm=P,ispd=c(0.25,0.75),Rho=R)
# The order of the states has got swapped; 
# note that ifelse(s.1[[1]]=="1","2","1") is much
# more similar to s.2[[1]] than is s.1[[1]].

## End(Not run)

Calculate fractions of missing values.

Description

Calculate the fraction (univariate case) or fractions (bivariate case) of missing values in the data or in each component of the data.

Usage

nafracCalc(y,drop=TRUE)

Arguments

Value

If y is not of class "multipleHmmDataSets", then the returned value is a scalar (between 0 and 1) if the data are univariate or a pair (2-vector) of such scalars if the data are bivariate. The values are equal to the ratios of the total count of missing values in the appropriate column to the total number of observations.

If y is of class "multipleHmmDataSets", and if y has length greater than 1 or drop is FALSE, then the returned value is a list of such scalars or 2-vectors, each corresponding to one of the data sets constituting y. If y has length equal to 1 and drop is TRUE, then the returned value is the same as if codey were not of class "multipleHmmDataSets".

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

rhmm() misstify()

Examples

xxx <- with(SydColDisc,split(y,f=list(locn,depth)))
nafracCalc(xxx) # 0.7185199

Probability of state sequences.

Description

Calculates the conditional probability of one or more state sequences, given the corresponding observations sequences (and the model parameters.

Usage

pr(s, y, model=NULL, tpm, Rho, ispd=NULL, warn=TRUE)

Arguments

Value

The probability of s given y, or a vector of such probabilities if s and y are lists, or a list of such vectors if y is of class "multipleHmmDataSets".

Warning

The conditional probabilities will be tiny if the sequences involved are of any substantial length. Underflow may be a problem. The implementation of the calculations is not sophisticated.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

hmm(), mps(), viterbi(), sp(), fitted.hmm.discnp()

Examples

## Not run: 
P <- matrix(c(0.7,0.3,0.1,0.9),2,2,byrow=TRUE)
R <- matrix(c(0.5,0,0.1,0.1,0.3,
              0.1,0.1,0,0.3,0.5),5,2)
set.seed(42)
lll   <- sample(250:350,20,TRUE)
set.seed(909)
y.num <- rhmm(ylengths=lll,nsim=1,tpm=P,Rho=R,drop=TRUE)
fit.num <- hmm(y.num,K=2,keep.y=TRUE,verb=TRUE)
# Using fitted parmeters.
s.vit.1   <- viterbi(y.num,fit.num)
pr.vit.1  <- pr(s.vit.1,model=fit.num)
# Using true parameters from which y.num was generated.
s.vit.2   <- viterbi(y.num,tpm=P,Rho=R)
pr.vit.2  <- pr(s.vit.2,y.num,tpm=P,Rho=R)
set.seed(202)
y.mult <- rhmm(fit.num,nsim=4)
s.vit.3 <- viterbi(y.mult,tpm=fit.num$tpm,Rho=fit.num$Rho)
pr.vit.3  <- pr(s.vit.3,y.mult,tpm=fit.num$tpm,Rho=fit.num$Rho)

## End(Not run)

Predicted values of a discrete non-parametric hidden Markov model.

Description

Calculates predicted values given a specification of a discrete non-parametric hidden Markov model. The specification may be provided in the form of a hmm.discnp object as returned by hmm() or in the form of “components” of such a model: the data y, the transition probability matrix tpm, the emission probabilities Rho, etc. If the data are numeric then these predicted values are the conditional expectations of the observations, given the entire observation sequence (and the — possibly estimated — parameters of the model). If the data are categorical (whence “expectations” make no sense) the “predicted values” are taken to be the probabilities of each of the possible values of the observations, at each time point.

Usage

## S3 method for class 'hmm.discnp'
predict(object, y = NULL, tpm=NULL, Rho=NULL,
                             ispd=NULL, X=NULL,addIntercept=NULL,
                             warn=TRUE, drop=TRUE, ...)

Arguments

Details

This function is essentially the same as fitted.hmm.discnp(). The main difference is that it allows the calculation of fitted/predicted values for a data object y possibly different from that to which the model was fitted. Note that if both the argument y and object[["y"]] are present, the “argument” value takes precedence. This function also allows the model to be specfied in terms of individual components rather than as a fitted model of class "hmm.discnp". These components, (tpm, Rho, ispd, X, addIntercept) if supplied, take precedence over the corresponding components of object. The opposite applies with sp(). The function fitted.hmm.discnp() makes use only of the components of object.

Value

See the help for fitted.hmm.discnp().

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

sp() link{fitted.hmm.discnp}()

Examples

P <- matrix(c(0.7,0.3,0.1,0.9),2,2,byrow=TRUE)
R <- matrix(c(0.5,0,0.1,0.1,0.3,
              0.1,0.1,0,0.3,0.5),5,2)
set.seed(42)
ll1  <- sample(250:350,20,TRUE)
y1   <- rhmm(ylengths=ll1,nsim=1,tpm=P,Rho=R,drop=TRUE)
fit  <- hmm(y1,K=2,verb=TRUE,keep.y=TRUE,itmax=10)
fv   <- fitted(fit)
set.seed(176)
ll2  <- sample(250:350,20,TRUE)
y2   <- rhmm(ylengths=ll2,nsim=1,tpm=P,Rho=R,drop=TRUE)
pv   <- predict(fit,y=y2)
yval <- letters[1:5]
set.seed(171)
y3   <- rhmm(ylengths=ll2,yval=yval,nsim=1,tpm=P,Rho=R,drop=TRUE)
fit3 <- hmm(y3,K=2,verb=TRUE,keep.y=TRUE,itmax=10)
pv3  <- predict(fit3) # Same as fitted(fit3).

Simulate discrete data from a non-parametric hidden Markov model.

Description

Simulates one or more replicates of discrete data from a model such as is fitted by the function hmm().

Usage

rhmm(model,...,nsim,verbose=FALSE)
## Default S3 method:
rhmm(model, ..., nsim=1, verbose=FALSE, ylengths,
                       nafrac=NULL, fep=NULL, tpm, Rho, ispd=NULL, yval=NULL,
                       drop=TRUE, forceNumeric=TRUE)
## S3 method for class 'hmm.discnp'
rhmm(model, ..., nsim=1, verbose=FALSE, inMiss=TRUE,
                          fep=NULL, drop=TRUE, forceNumeric=TRUE)

Arguments

Value

If nsim>1 or drop is FALSE then the value returned is a list of length nsim. Each entry of this list is in turn a list of the same length as ylengths, each component of which is an independent vector or matrix of simulated observations. The length or number of rows of component i of this list is equal to ylengths[i]. The values of the observations are entries of yval or of its entries when yval is a list.

If nsim=1 and drop is TRUE then the (“outer”) list described above is replaced by its first and only entry

If the length of ylengths is 1 and drop is TRUE then each “inner” list described above is replaced by its first and only entry.

Note

You may find it useful to avail yourself of the function nafracCalc() to determine the fraction of missing values in a given existing (presumably “real”) data set.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

hmm() nafracCalc() misstify()

Examples

# To do: one or more bivariate examples.
## Not run: 
    y <- list(linLandFlows$deciles,ftLiardFlows$deciles)
    fit <- hmm(y,K=3)
    simX <- rhmm(fit)

## End(Not run)

Simulation based covariance matrix.

Description

Produces an estimate of the covariance matrix of the parameter estimates in a model fitted by hmm.discnp. Uses a method based on simulation (or “parametric bootstrapping”).

Usage

scovmat(object, expForm=TRUE, seed = NULL, nsim=100, verbose = TRUE)

Arguments

Details

This function is currently applicable only to models fitted to univariate data. If there are predictors in the model, then only the exponential form of the parameters may be used, i.e. expForm must be TRUE.

Value

A (positive definite) matrix which is an estimate of the covariance of the parameter estimates from the fitted model specified by object. It has row and column labels which indicate the parameters to which its entries pertain, in a reasonably perspicuous manner.

This matrix has an attribute seed (the random number generation seed that was used) so that the calculations can be reproduced.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

squantCI() link{rhmm}() link{hmm)}()

Examples

## Not run: 
y   <- list(lindLandFlows$deciles,ftLiardFlows$deciles)
fit <- hmm(y,K=3)
ccc <- scovmat(fit,nsim=100)

## End(Not run)

Calculate the conditional state probabilities.

Description

Returns the probabilities that the underlying hidden state is equal to each of the possible state values, at each time point, given the observation sequence.

Usage

sp(y, model = NULL, tpm=NULL, Rho=NULL, ispd=NULL, X=NULL,
   addIntercept=NULL, warn=TRUE, drop=TRUE)

Arguments

Details

Note that in contrast to predict.hmm.discnp(), components in model take precendence over individually supplied components (tpm, Rho, ispd, X and addIntercept).

Value

If y is a single matrix of observations or a list of length 1, and if drop is TRUE then the returned value is a matrix whose rows correspond to the states of the hidden Markov chain, and whose columns correspond to the observation times. Otherwise the returned value is a list of such matrices, one for each matrix of observations.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

hmm(), mps(), viterbi(), pr(), fitted.hmm.discnp()

Examples

P <- matrix(c(0.7,0.3,0.1,0.9),2,2,byrow=TRUE)
R <- matrix(c(0.5,0,0.1,0.1,0.3,
              0.1,0.1,0,0.3,0.5),5,2)
set.seed(42)
y    <- rhmm(ylengths=rep(300,20),nsim=1,tpm=P,Rho=R,drop=TRUE)
fit  <- hmm(y,K=2,verb=TRUE,keep.y=TRUE,itmax=10)
cpe1 <- sp(model=fit) # Using the estimated parameters.
cpe2 <- sp(y,tpm=P,Rho=R,warn=FALSE) # Using the ``true'' parameters.
# The foregoing would issue a warning that Rho had no row names
# were it not for the fact that "warn" has been set to FALSE.

Simulation-quantile based confidence intervals.

Description

Calculates estimates of confidence intervals for the parameters of a model fitted by hmm.discnp. Uses a method based quantiles of estimates produced by simulation (or “parametric bootstrapping”).

Usage

squantCI(object, expForm = TRUE, seed = NULL, alpha = 0.05,
         nsim=100, verbose = TRUE)

Arguments

Details

This function is currently applicable only to models fitted to univariate data. If there are predictors in the model, then only the exponential form of the parameters may be used, i.e. expForm must be TRUE.

Value

A npar-by-2 matrix (where npar is the number of “independent” parameters in the model) whose rows form the estimated confidence intervals. (The first entry of each row is the lower bound of a confidence interval for the corresponding parameter, and the second entry is the upper bound. The row labels indicate the parameters to which each row pertains, in a reasonably perspicuous manner. The column labels indicate the relevant quantiles in percentages.

This matrix has an attribute seed (the random number generation seed that was used) so that the calculations can be reproduced.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

scovmat() link{rhmm}() link{hmm)}()

Examples

## Not run: 
y   <- list(lindLandFlows$deciles,ftLiardFlows$deciles)
fit <- hmm(y,K=3)
CIs <- squantCI(fit,nsim=100)

## End(Not run)

Update a fitted hmm.discnp model.

Description

An update() method for objects of class hmm.discnp.

Usage

## S3 method for class 'hmm.discnp'
update(object,..., data, Kplus1=FALSE,
        tpm2=NULL, verbose=FALSE, method=NULL, optimiser=NULL,
        stationary=NULL, mixture=NULL, cis=NULL, tolerance=NULL,
        itmax=NULL, crit=NULL, X=NULL, addIntercept=NULL)

Arguments

Details

Except for argument X, any arguments that are left NULL have their values supplied from the args component of object.

Value

An object of class hmm.discnp with an additional component init.log.like which is the initial log likelihood calculated at the starting values of the parameters (which may be modified from the parameters returned in the object being updated, if Kplus1 is TRUE). The calculation is done by the function logLikHmm(). Barring the strange and unforeseen, init.log.like should be (reassuringly) equal to object$log.like. See hmm() for details of the other components of the returned value.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

See Also

hmm() rhmm.hmm.discnp()

Examples

set.seed(294)
fit  <- hmm(WoodPeweeSong,K=2,rand.start=list(tpm=TRUE,Rho=TRUE),itmax=10)
xxx  <- rhmm(fit,nsim=1)
sfit <- update(fit,data=xxx,itmax=10)
yyy  <- with(SydColDisc,split(y,f=list(locn,depth)))
f1   <- hmm(yyy,K=1)
f2   <- update(f1,data=yyy,Kplus1=TRUE) # Big improvement, but ...
## Not run: 
g2   <- hmm(yyy,K=2) # Substantially better than f2. 

## End(Not run)

Most probable state sequence.

Description

Calculates “the” most probable state sequence underlying each of one or more replicate observation sequences.

Usage

viterbi(y, model = NULL, tpm, Rho, ispd=NULL,log=FALSE, warn=TRUE)

Arguments

Details

Applies the Viterbi algorithm to calculate “the” most probable robable state sequence underlying each observation sequences.

Value

If y consists of a single observation sequence, the value is the underlying most probable observation sequence, or a matrix whose columns consist of such sequences if there is more than one (equally) most probable sequence.

If y consists of a list of observation sequences, the value is a list each entry of which is of the form described above.

If y is of class "multipleHmmDataSets" then the value returned is a list of lists of the sort described above.

Warning

There may be more than one equally most probable state sequence underlying a given observation sequence. This phenomenon can occur but appears to be unlikely to do so in practice.

Thanks

The correction made to the code so as to avoid underflow problems was made due to an inquiry and suggestion from Owen Marshall.

Author(s)

Rolf Turner r.turner@auckland.ac.nz

References

Rabiner, L. R., "A tutorial on hidden Markov models and selected applications in speech recognition," Proc. IEEE vol. 77, pp. 257 – 286, 1989.

See Also

hmm(), rhmm(), mps(), pr()

Examples

# See the help for logLikHmm() for how to generate y.num and y.let.
## Not run: 
fit.num     <- hmm(y.num,K=2,verb=TRUE,keep.y=TRUE)
v.1         <- viterbi(model=fit.num)
rownames(R) <- 1:5 # Avoids a (harmless) warning.
v.2         <- viterbi(y.num,tpm=P,Rho=R)
# P and R as in the help for logLikHmm() and for sp().

# Note that the order of the states has gotten swapped; 3-v.1[[1]]
# is identical to v.2[[1]]; for other k = 2, ..., 20, 3-v.1[[k]]
# is much more similar to v.2[[k]] than is v.1[[k]].

fit.let <- hmm(y.let,K=2,verb=TRUE,keep.y=TRUE))
v.3 <- viterbi(model=fit.let)
rownames(R) <- letters[1:5]
v.4 <- viterbi(y.let,tpm=P,Rho=R)

## End(Not run)

Data from “An Introduction to Discrete-Valued Time Series”

Description

Data sets from the book “An Introduction to Discrete-Valued Time Series” by Christian H. Weiß.

Usage

    data(Bovine)
    data(Cryptosporidiosis)
    data(Downloads)
    data(EricssonB_Jul2)
    data(FattyLiver)
    data(FattyLiver2)
    data(goldparticle380)
    data(Hanta)
    data(InfantEEGsleepstates)
    data(IPs)
    data(LegionnairesDisease)
    data(OffshoreRigcountsAlaska)
    data(PriceStability)
    data(Strikes)
    data(WoodPeweeSong)

Format

• Bovine A character vector of length 8419.

• Cryptosporidiosis A numeric (integer) vector of length 365.

• Downloads A numeric (integer) vector of length 267.

• EricssonB_Jul2 A numeric (integer) vector of length 460.

• FattyLiver2 A numeric (integer) vector of length 449.

• FattyLiver A numeric (integer) vector of length 928.

• goldparticle380 A numeric (integer) vector of length 380.

• Hanta A numeric (integer) vector of length 52.

• InfantEEGsleepstates A character vector of length 107.

• IPs A numeric (integer) vector of length 241.

• LegionnairesDisease A numeric (integer) vector of length 365.

• OffshoreRigcountsAlaska A numeric (integer) vector of length 417.

• PriceStability A numeric (integer) vector of length 152.

• Strikes A numeric (integer) vector of length 108.

• WoodPeweeSong A numeric (integer) vector of length 1327.

Details

For detailed information about each of these data sets, see the book cited in the References.

Note that the data sets Cryptosporidiosis and LegionnairesDisease are actually called
Cryptosporidiosis_02-08 and LegionnairesDisease_02-08 in the given reference. The
“suffixes” were removed since the minus sign causes problems in a variable name in R.

Source

These data sets were kindly provided by Prof. Christian H. Weiß. The package author is also pleased to acknowledge the kind permission granted by Prof. Kurt Brännäs (Professor Emeritus of Economics at Umeå University) to include the Ericsson time series data set (EricssonB_Jul2).

References

Christian H. Weiß (2018). An Introduction to Discrete-Valued Time Series. Chichester: John Wiley & Sons.

Examples

## Not run: 
fit1 <- hmm(WoodPeweeSong,K=2,verbose=TRUE)
# EM converges in 6 steps --- suspicious.
set.seed(321)
fit2 <- hmm(WoodPeweeSong,K=2,verbose=TRUE,rand.start=list(tpm=TRUE,Rho=TRUE))
# 52 steps --- note the huge difference between fit1$log.like and fit2$log.like!
set.seed(321)
fit3 <- hmm(WoodPeweeSong,K=2,verbose=TRUE,method="bf",
            rand.start=list(tpm=TRUE,Rho=TRUE))
# log likelihood essentially the same as for fit2

## End(Not run)