| Title: | Inference for Hidden Markov and Semi-Markov Models |
|---|---|
| Description: | Parameter estimation and prediction for hidden Markov and semi-Markov models for data with multiple observation sequences. Suitable for equidistant time series data, with multivariate and/or missing data. Allows user defined emission distributions. |
| Authors: | Jared O'Connell <[email protected]>, Søren Højsgaard <[email protected]> |
| Maintainer: | Jared O'Connell <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.4.21 |
| Built: | 2024-12-15 07:20:06 UTC |
| Source: | CRAN |
Add a colour coded horizontal bar representing the state sequence to a plot of (presumably time-series) data.
addStates(states, x=NULL,ybot = axTicks(2)[1], ytop=ybot + (axTicks(2)[2] - axTicks(2)[1])/5,dy = ytop - ybot, greyscale = FALSE,leg = NA, J = length(unique(states)), time.scale = 1, shiftx = 0)addStates(states, x=NULL,ybot = axTicks(2)[1], ytop=ybot + (axTicks(2)[2] - axTicks(2)[1])/5,dy = ytop - ybot, greyscale = FALSE,leg = NA, J = length(unique(states)), time.scale = 1, shiftx = 0)
states |
A vector of integers representing the states traversed |
x |
The time values where the states are observed ((1:length(states)-shiftx)/time.scale if NULL) |
ybot |
Vertical bottom limit of the bar. |
ytop |
Vertical top limit of the bar. |
dy |
Height of the bar. |
greyscale |
If TRUE produces a bar in greyscale. |
leg |
Array of state names, if present, produces a legend. |
J |
Number of states |
time.scale |
Resolution of the timescale |
shiftx |
Shift the bar forward or backwards horizontal by shiftx distance. |
Soren Hojsgaard [email protected]
addStates
plot(rnorm(100),type='l') addStates(rep(c(1,2),each=50)) plot(seq(0.01,1,.01),rnorm(100),type='l') addStates(rep(c(1,2),each=50),seq(0.01,1,.01))plot(rnorm(100),type='l') addStates(rep(c(1,2),each=50)) plot(seq(0.01,1,.01),rnorm(100),type='l') addStates(rep(c(1,2),each=50),seq(0.01,1,.01))
Calculates the density of observations x for state j given the parameters in model. This is used for
a multivariate Gaussian emission distribution of a HMM or HSMM and is a suitable prototype for user's to make their own custom distributions.
dmvnorm.hsmm(x, j, model)dmvnorm.hsmm(x, j, model)
x |
Observed value |
j |
State |
model |
A |
This is used by hmm and hsmm to calculate densities for use in the E-step of the EM algorithm.
It can also be used as a template for users wishing to building their own emission distributions
A vector of probability densities.
Jared O'Connell [email protected]
J<-2 initial <- rep(1/J,J) P <- matrix(c(.3,.5,.7,.5),nrow=J) b <- list(mu=list(c(-3,0),c(1,2)),sigma=list(diag(2),matrix(c(4,2,2,3), ncol=2))) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dmvnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rmvnorm.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.mvnorm)J<-2 initial <- rep(1/J,J) P <- matrix(c(.3,.5,.7,.5),nrow=J) b <- list(mu=list(c(-3,0),c(1,2)),sigma=list(diag(2),matrix(c(4,2,2,3), ncol=2))) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dmvnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rmvnorm.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.mvnorm)
Calculates the density of observations x for state j given the parameters in model. This is used for
the Gaussian emission distribution of a HMM or HSMM and is a suitable prototype for user's to make their own custom distributions.
dnorm.hsmm(x, j, model)dnorm.hsmm(x, j, model)
x |
Observed value |
j |
State |
model |
A |
This is used by hmm and hsmm to calculate densities for use in the E-step of the EM algorithm.
It can also be used as a template for users wishing to building their own emission distributions
A vector of probability densities.
Jared O'Connell [email protected]
Calculates the density of observations x for state j given the parameters in model. This is used for
a Poisson emission distribution of a HMM or HSMM and is a suitable prototype for user's to make their own custom distributions.
dpois.hsmm(x, j, model)dpois.hsmm(x, j, model)
x |
Observed value |
j |
State |
model |
A |
This is used by hmm and hsmm to calculate densities for use in the E-step of the EM algorithm.
It can also be used as a template for users wishing to building their own emission distributions
A vector of probability densities.
Jared O'Connell [email protected]
J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(lambda=c(1,3,6)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dpois.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rpois.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.pois)J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(lambda=c(1,3,6)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dpois.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rpois.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.pois)
Estimates parameters for the Gamma distribution using the Method of Maximum Likelihood, works with weighted data.
gammafit(x, wt = NULL)gammafit(x, wt = NULL)
x |
A vector of observations |
wt |
Optional set of weights |
shape |
The shape parameter |
scale |
The scale parameter (equal to 1/rate) |
Jared O'Connell [email protected]
Choi, S. and Wette, R. (1969), Maximum likelihood estimation of the parameters of the gamma distribution and their bias, Technometrics, 11, 683-96-690.
gammafit(rgamma(1000,shape=10,scale=13))gammafit(rgamma(1000,shape=10,scale=13))
Estimates parameters of a HMM using the EM algorithm.
hmmfit(x,start.val,mstep=mstep.norm,lock.transition=FALSE,tol=1e-08,maxit=1000)hmmfit(x,start.val,mstep=mstep.norm,lock.transition=FALSE,tol=1e-08,maxit=1000)
x |
A hsmm.data object (see Details) |
start.val |
Starting parameters for the model (see Details) |
mstep |
Re-estimates the parameters of density function on each iteration |
lock.transition |
If TRUE will not re-estimate the transition matrix |
maxit |
Maximum number of iterations |
tol |
Convergence tolerance |
start |
A vector of the starting probabilities for each state |
a |
The transition matrix of the embedded Markov chain |
emission |
A list of the parameters of the emission distribution |
Jared O'Connell [email protected]
Jared O'Connell, Soren Hojsgaard (2011). Hidden Semi Markov Models for Multiple Observation Sequences: The mhsmm Package for R., Journal of Statistical Software, 39(4), 1-22., URL http://www.jstatsoft.org/v39/i04/.
Rabiner, L. (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE, 77, 257-286.
J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rnorm.hsmm) plot(train,xlim=c(0,100)) init0 <- rep(1/J,J) P0 <- matrix(1/J,nrow=J,ncol=J) b0 <- list(mu=c(-3,1,3),sigma=c(1,1,1)) startval <- hmmspec(init=init0, trans=P0,parms.emission=b0,dens.emission=dnorm.hsmm) h1 = hmmfit(train,startval,mstep=mstep.norm) plot(h1$loglik,type='b',ylab='Log-likelihood',xlab='Iteration') summary(h1) #proportion of incorrect states mean(train$s!=predict(h1,train)$s) #simulate a new test set test <- simulate(model, nsim=c(100,200,300), seed=1234,rand.emis=rnorm.hsmm) mean(test$s!=predict(h1,test)$s)J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rnorm.hsmm) plot(train,xlim=c(0,100)) init0 <- rep(1/J,J) P0 <- matrix(1/J,nrow=J,ncol=J) b0 <- list(mu=c(-3,1,3),sigma=c(1,1,1)) startval <- hmmspec(init=init0, trans=P0,parms.emission=b0,dens.emission=dnorm.hsmm) h1 = hmmfit(train,startval,mstep=mstep.norm) plot(h1$loglik,type='b',ylab='Log-likelihood',xlab='Iteration') summary(h1) #proportion of incorrect states mean(train$s!=predict(h1,train)$s) #simulate a new test set test <- simulate(model, nsim=c(100,200,300), seed=1234,rand.emis=rnorm.hsmm) mean(test$s!=predict(h1,test)$s)
Creates a model specficiation for a hidden Markov model
hmmspec(init, trans, parms.emission, dens.emission, rand.emission=NULL,mstep=NULL)hmmspec(init, trans, parms.emission, dens.emission, rand.emission=NULL,mstep=NULL)
init |
Distribution of states at t=1 ie. P(S=s) at t=1 |
trans |
The transition matrix of the Markov chain |
parms.emission |
A list containing the parameters of the emission distribution |
dens.emission |
Density function of the emission distribution. |
rand.emission |
The function used to generate observations from the emission distribution |
mstep |
Re-estimates the parameters of density function on each iteration |
A hmmspec object
Jared O'Connell [email protected]
Jared O'Connell, Soren Hojsgaard (2011). Hidden Semi Markov Models for Multiple Observation Sequences: The mhsmm Package for R., Journal of Statistical Software, 39(4), 1-22., URL http://www.jstatsoft.org/v39/i04/.
Rabiner, L. (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE, 77, 257-286.
simulate.hmmspec, simulate.hmmspec, hmmfit, predict.hmm
Estimates parameters of a HSMM using the EM algorithm.
hsmmfit(x,model,mstep=NULL,M=NA,maxit=100, lock.transition=FALSE,lock.d=FALSE,graphical=FALSE)hsmmfit(x,model,mstep=NULL,M=NA,maxit=100, lock.transition=FALSE,lock.d=FALSE,graphical=FALSE)
x |
A hsmm.data object (see Details) |
model |
Starting parameters for the model (see |
mstep |
Re-estimates the parameters of density function on each iteration |
maxit |
Maximum number of iterations |
M |
Maximum number of time spent in a state (truncates the waiting distribution) |
lock.transition |
If TRUE will not re-estimate the transition matrix |
lock.d |
If TRUE will not re-estimate the sojourn time density |
graphical |
If TRUE will plot the sojourn densities on each iteration |
start |
A vector of the starting probabilities for each state |
a |
The transition matrix of the embedded Markov chain |
emission |
A list of the parameters of the emission distribution |
waiting |
A list of the parameters of the waiting distribution |
Jared O'Connell [email protected]
Jared O'Connell, Soren Hojsgaard (2011). Hidden Semi Markov Models for Multiple Observation Sequences: The mhsmm Package for R., Journal of Statistical Software, 39(4), 1-22., URL http://www.jstatsoft.org/v39/i04/.
Guedon, Y. (2003), Estimating hidden semi-Markov chains from discrete sequences, Journal of Computational and Graphical Statistics, Volume 12, Number 3, page 604-639 - 2003
hsmmspec,simulate.hsmmspec,predict.hsmm
J <- 3 init <- c(0,0,1) P <- matrix(c(0,.1,.4,.5,0,.6,.5,.9,0),nrow=J) B <- list(mu=c(10,15,20),sigma=c(2,1,1.5)) d <- list(lambda=c(10,30,60),shift=c(10,100,30),type='poisson') model <- hsmmspec(init,P,parms.emission=B,sojourn=d,dens.emission=dnorm.hsmm) train <- simulate(model,r=rnorm.hsmm,nsim=100,seed=123456) plot(train,xlim=c(0,400)) start.poisson <- hsmmspec(init=rep(1/J,J), transition=matrix(c(0,.5,.5,.5,0,.5,.5,.5,0),nrow=J), parms.emission=list(mu=c(4,12,23), sigma=c(1,1,1)), sojourn=list(lambda=c(9,25,40),shift=c(5,95,45),type='poisson'), dens.emission=dnorm.hsmm) M=500 # not run (takes some time) # h.poisson <- hsmmfit(train,start.poisson,mstep=mstep.norm,M=M) # plot(h.poisson$loglik,type='b',ylab='Log-likelihood',xlab='Iteration') ##has it converged? # summary(h.poisson) # predicted <- predict(h.poisson,train) # table(train$s,predicted$s) ##classification matrix # mean(predicted$s!=train$s) ##misclassification rate d <- cbind(dunif(1:M,0,50),dunif(1:M,100,175),dunif(1:M,50,130)) start.np <- hsmmspec(init=rep(1/J,J), transition=matrix(c(0,.5,.5,.5,0,.5,.5,.5,0),nrow=J), parms.emission=list(mu=c(4,12,23), sigma=c(1,1,1)), sojourn=list(d=d,type='nonparametric'), dens.emission=dnorm.hsmm) # not run (takes some time) # h.np <- hsmmfit(train,start.np,mstep=mstep.norm,M=M,graphical=TRUE) # mean(predicted$s!=train$s) ##misclassification rate #J <- 2 #init <- c(1, 0) #P <- matrix(c(0, 1, 1, 0), nrow = J) #B <- list(mu = list(c(2, 3), c(3, 4)), sigma = list(matrix(c(4, 2, 2, 3), ncol = 2), diag(2))) #d <- list(shape = c(10, 25), scale = c(2, 2), type = "gamma") #model <- hsmmspec(init, P, parms.emis = B, sojourn = d, dens.emis = dmvnorm.hsmm) #train <- simulate(model, c(1000,100), seed = 123, rand.emis = rmvnorm.hsmm) #yhat <- predict(model, train) #mean(predict(model,train)$s==train$s)J <- 3 init <- c(0,0,1) P <- matrix(c(0,.1,.4,.5,0,.6,.5,.9,0),nrow=J) B <- list(mu=c(10,15,20),sigma=c(2,1,1.5)) d <- list(lambda=c(10,30,60),shift=c(10,100,30),type='poisson') model <- hsmmspec(init,P,parms.emission=B,sojourn=d,dens.emission=dnorm.hsmm) train <- simulate(model,r=rnorm.hsmm,nsim=100,seed=123456) plot(train,xlim=c(0,400)) start.poisson <- hsmmspec(init=rep(1/J,J), transition=matrix(c(0,.5,.5,.5,0,.5,.5,.5,0),nrow=J), parms.emission=list(mu=c(4,12,23), sigma=c(1,1,1)), sojourn=list(lambda=c(9,25,40),shift=c(5,95,45),type='poisson'), dens.emission=dnorm.hsmm) M=500 # not run (takes some time) # h.poisson <- hsmmfit(train,start.poisson,mstep=mstep.norm,M=M) # plot(h.poisson$loglik,type='b',ylab='Log-likelihood',xlab='Iteration') ##has it converged? # summary(h.poisson) # predicted <- predict(h.poisson,train) # table(train$s,predicted$s) ##classification matrix # mean(predicted$s!=train$s) ##misclassification rate d <- cbind(dunif(1:M,0,50),dunif(1:M,100,175),dunif(1:M,50,130)) start.np <- hsmmspec(init=rep(1/J,J), transition=matrix(c(0,.5,.5,.5,0,.5,.5,.5,0),nrow=J), parms.emission=list(mu=c(4,12,23), sigma=c(1,1,1)), sojourn=list(d=d,type='nonparametric'), dens.emission=dnorm.hsmm) # not run (takes some time) # h.np <- hsmmfit(train,start.np,mstep=mstep.norm,M=M,graphical=TRUE) # mean(predicted$s!=train$s) ##misclassification rate #J <- 2 #init <- c(1, 0) #P <- matrix(c(0, 1, 1, 0), nrow = J) #B <- list(mu = list(c(2, 3), c(3, 4)), sigma = list(matrix(c(4, 2, 2, 3), ncol = 2), diag(2))) #d <- list(shape = c(10, 25), scale = c(2, 2), type = "gamma") #model <- hsmmspec(init, P, parms.emis = B, sojourn = d, dens.emis = dmvnorm.hsmm) #train <- simulate(model, c(1000,100), seed = 123, rand.emis = rmvnorm.hsmm) #yhat <- predict(model, train) #mean(predict(model,train)$s==train$s)
Creates a model specification of a hidden semi-Markov model.
hsmmspec(init,transition,parms.emission,sojourn,dens.emission, rand.emission=NULL,mstep=NULL)hsmmspec(init,transition,parms.emission,sojourn,dens.emission, rand.emission=NULL,mstep=NULL)
init |
Distribution of states at t=1 ie. P(S=s) at t=1 |
transition |
The transition matrix of the embedded Markov chain (diagonal must be 0) |
parms.emission |
A list containing the parameters of the emission distribution |
sojourn |
A list containining the parameters and type of sojourn distribtuion (see Details) |
dens.emission |
Density function of the emission distribution |
rand.emission |
The function used to generate observations from the emission distribution |
mstep |
Re-estimates the parameters of density function on each iteration |
The sojourn argument provides a list containing the parameters for the available sojourn distributions. Available sojourn distributions are shifted Poisson, Gamma and non-parametric.
In the case of the Gamma distribution, sojourn is a list with vectors shape and scale (the Gamma parameters in dgamma), both of length J. Where J is the number of states. See reprocows for an example using Gamma sojourn distributions.
In the case of shifted Poisson, sojourn is list with vectors shift and lambda, both of length J. See hsmmfit for an example using shifted Poisson sojourn distributions.
In the case of non-parametric, sojourn is a list containing a M x J matrix. Where entry (i,j) is the probability of a sojourn of length i in state j. See hsmmfit for an example using shifted non-parametric sojourn distributions.
An object of class hsmmspec
Jared O'Connell [email protected]
Jared O'Connell, Soren Hojsgaard (2011). Hidden Semi Markov Models for Multiple Observation Sequences: The mhsmm Package for R., Journal of Statistical Software, 39(4), 1-22., URL http://www.jstatsoft.org/v39/i04/.
Guedon, Y. (2003), Estimating hidden semi-Markov chains from discrete sequences, Journal of Computational and Graphical Statistics, Volume 12, Number 3, page 604-639 - 2003
hsmmfit ,simulate.hsmmspec, predict.hsmm
Re-estimates the parameters of a multivariate normal emission
distribution as part of the EM algorithm for HMMs and HSMMs.
This is called by the hmm and hsmm functions. It is a
suitable prototype function for users wishing to design their own
emission distributions.
mstep.mvnorm(x, wt)mstep.mvnorm(x, wt)
x |
A vector of observed values |
wt |
A T x J matrix of weights. Column entries are the weights for respective states. |
Users may write functions that take the same arguments and return the same values for their own custom emission distributions.
Returns the emission slot of a hmmspec or hsmmspec object
mu |
A list of length J contain the mean vectors |
sigma |
A list of length J containing the covariance matrices |
Jared O'Connell [email protected]
J<-2 initial <- rep(1/J,J) P <- matrix(c(.3,.5,.7,.5),nrow=J) b <- list(mu=list(c(-3,0),c(1,2)),sigma=list(diag(2),matrix(c(4,2,2,3), ncol=2))) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dmvnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rmvnorm.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.mvnorm)J<-2 initial <- rep(1/J,J) P <- matrix(c(.3,.5,.7,.5),nrow=J) b <- list(mu=list(c(-3,0),c(1,2)),sigma=list(diag(2),matrix(c(4,2,2,3), ncol=2))) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dmvnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rmvnorm.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.mvnorm)
Re-estimates the parameters of a normal emission distribution as part of the EM algorithm for HMMs and HSMMs.
This is called by the hmm and hsmm functions. It is a suitable prototype function for users
wishing to design their own emission distributions.
mstep.norm(x, wt)mstep.norm(x, wt)
x |
A vector of observed values |
wt |
A T x J matrix of weights. Column entries are the weights for respective states. |
Users may write functions that take the same arguments and return the same values for their own custom emission distributions.
Returns the emission slot of a hmmspec or hsmmspec object
mu |
Vector of length J contain the means |
sigma |
Vector of length J containing the variances |
Jared O'Connell [email protected]
Re-estimates the parameters of a Poisson emission distribution as part
of the EM algorithm for HMMs and HSMMs. This is called by the
hmm and hsmm functions. It is a suitable prototype
function for users wishing to design their own emission
distributions.
mstep.pois(x, wt)mstep.pois(x, wt)
x |
A vector of observed values |
wt |
A T x J matrix of weights. Column entries are the weights for respective states. |
Users may write functions that take the same arguments and return the same values for their own custom emission distributions.
Returns the emission slot of a hmmspec or hsmmspec object
lambda |
Vector of length J containing the Poisson parameters for each state j |
Jared O'Connell [email protected]
J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(lambda=c(1,3,6)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dpois.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rpois.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.pois)J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(lambda=c(1,3,6)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dpois.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rpois.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.pois)
Displays the densities for the sojourn distributions of each state.
## S3 method for class 'hsmm' plot(x, ...)## S3 method for class 'hsmm' plot(x, ...)
x |
A |
... |
Arguments passed to |
Jared O'Connell [email protected]
Produces a plot of the observed sequences, and displays a coloured bar signifying the hidden states (if available)
## S3 method for class 'hsmm.data' plot(x, ...)## S3 method for class 'hsmm.data' plot(x, ...)
x |
A |
... |
Arguments passed to |
Jared O'Connell [email protected]
J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b, dens.emission=dnorm.hsmm) train <- simulate(model, nsim=300, seed=1234, rand.emis=rnorm.hsmm) plot(train,xlim=c(0,100))J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b, dens.emission=dnorm.hsmm) train <- simulate(model, nsim=300, seed=1234, rand.emis=rnorm.hsmm) plot(train,xlim=c(0,100))
Predicts the underlying state sequence for an observed sequence newdata given a hmm model
## S3 method for class 'hmm' predict(object, newdata,method = "viterbi", ...)## S3 method for class 'hmm' predict(object, newdata,method = "viterbi", ...)
object |
An object of class |
newdata |
A vector or list of observations |
method |
Prediction method (see details) |
... |
further arguments passed to or from other methods. |
If method="viterbi", this technique applies the Viterbi algorithm for HMMs, producing the most likely sequence of states given the observed data. If method="smoothed", then the individually most likely (or smoothed) state sequence is produced, along with a matrix with the respective probabilities for each state.
Returns a hsmm.data object, suitable for plotting.
newdata |
A vector or list of observations |
s |
A vector containing the reconstructed state sequence |
N |
The lengths of each sequence |
p |
A matrix where the rows represent time steps and the columns are the probability for the respective state (only produced when |
Jared O'Connell [email protected]
Rabiner, L. (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE, 77, 257-286.
hmmfit,hmmspec
##See examples in 'hmmfit'##See examples in 'hmmfit'
Predicts the underlying state sequence for an observed sequence newdata given a hmmspec model
## S3 method for class 'hmmspec' predict(object, newdata,method = "viterbi", ...)## S3 method for class 'hmmspec' predict(object, newdata,method = "viterbi", ...)
object |
An object of class |
newdata |
A vector or list of observations |
method |
Prediction method (see details) |
... |
further arguments passed to or from other methods. |
If method="viterbi", this technique applies the Viterbi
algorithm for HMMs, producing the most likely sequence of states
given the observed data. If method="smoothed", then the
individually most likely (or smoothed) state sequence is produced,
along with a matrix with the respective probabilities for each
state. This function differs from predict.hmm in that it takes the
output from hmmspec ie. this is useful when users already know their
parameters and wish to make predictions.
Returns a hsmm.data object, suitable for plotting.
newdata |
A vector or list of observations |
s |
A vector containing the reconstructed state sequence |
N |
The lengths of each sequence |
p |
A matrix where the rows represent time steps and the columns are the probability for the respective state (only produced when |
Jared O'Connell [email protected]
Rabiner, L. (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE, 77, 257-286.
hmmspec
J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dnorm.hsmm) train <- simulate(model, nsim=300, seed=1234, rand.emis=rnorm.hsmm) mean(predict(model,train)$s!=train$s) #error rate when true model is knownJ<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dnorm.hsmm) train <- simulate(model, nsim=300, seed=1234, rand.emis=rnorm.hsmm) mean(predict(model,train)$s!=train$s) #error rate when true model is known
Predicts the underlying state sequence for an observed sequence newdata given a hsmm model
## S3 method for class 'hsmm' predict(object, newdata, method = "viterbi", ...)## S3 method for class 'hsmm' predict(object, newdata, method = "viterbi", ...)
object |
An object of type |
newdata |
A vector or dataframe of observations |
method |
Prediction method (see details) |
... |
further arguments passed to or from other methods. |
If method="viterbi", this technique applies the Viterbi algorithm for HSMMs, producing the most likely sequence of states given the observed data. If method="smoothed", then the individually most likely (or smoothed) state sequence is produced, along with a matrix with the respective probabilities for each state.
Returns a hsmm.data object, suitable for plotting.
newdata |
A vector or list of observations |
s |
A vector containing the reconstructed state sequence |
N |
The lengths of each sequence |
p |
A matrix where the rows represent time steps and the columns are the probability for the respective state (only produced when |
Jared O'Connell [email protected]
Guedon, Y. (2003), Estimating hidden semi-Markov chains from discrete sequences, Journal of Computational and Graphical Statistics, Volume 12, Number 3, page 604-639 - 2003
##See 'hsmmfit' for examples##See 'hsmmfit' for examples
Predicts the underlying state sequence for an observed sequence newdata given a hsmm model
## S3 method for class 'hsmmspec' predict(object, newdata, method = "viterbi",M=NA, ...)## S3 method for class 'hsmmspec' predict(object, newdata, method = "viterbi",M=NA, ...)
object |
An object of type |
newdata |
A vector or dataframe of observations |
method |
Prediction method (see details) |
M |
Maximum number of time spent in a state (truncates the waiting distribution) |
... |
further arguments passed to or from other methods. |
If method="viterbi", this technique applies the Viterbi
algorithm for HSMMs, producing the most likely sequence of states
given the observed data. If method="smoothed", then the
individually most likely (or smoothed) state sequence is produced,
along with a matrix with the respective probabilities for each
state. This method is different to predict.hsmm in that it takes
the output from hsmmspec as input ie. it is useful for people
who already know their model parameters.
Returns a hsmm.data object, suitable for plotting.
newdata |
A vector or list of observations |
s |
A vector containing the reconstructed state sequence |
N |
The lengths of each sequence |
p |
A matrix where the rows represent time steps and the columns are the probability for the respective state (only produced when |
Jared O'Connell [email protected]
Guedon, Y. (2003), Estimating hidden semi-Markov chains from discrete sequences, Journal of Computational and Graphical Statistics, Volume 12, Number 3, page 604-639 - 2003
J <- 3 init <- c(0,0,1) P <- matrix(c(0,.1,.4,.5,0,.6,.5,.9,0),nrow=J) B <- list(mu=c(10,15,20),sigma=c(2,1,1.5)) d <- list(lambda=c(10,30,60),shift=c(10,100,30),type='poisson') model <- hsmmspec(init,P,parms.emission=B,sojourn=d,dens.emission=dnorm.hsmm) train <- simulate(model,r=rnorm.hsmm,nsim=100,seed=123456) mean(predict(model,train,M=500)$s!=train$s) #error rate when true model is knownJ <- 3 init <- c(0,0,1) P <- matrix(c(0,.1,.4,.5,0,.6,.5,.9,0),nrow=J) B <- list(mu=c(10,15,20),sigma=c(2,1,1.5)) d <- list(lambda=c(10,30,60),shift=c(10,100,30),type='poisson') model <- hsmmspec(init,P,parms.emission=B,sojourn=d,dens.emission=dnorm.hsmm) train <- simulate(model,r=rnorm.hsmm,nsim=100,seed=123456) mean(predict(model,train,M=500)$s!=train$s) #error rate when true model is known
Prints the slots of a hmm object
## S3 method for class 'hmm' print(x, ...)## S3 method for class 'hmm' print(x, ...)
x |
An object of type |
... |
further arguments passed to or from other methods. |
Jared O'Connell [email protected]
Prints the parameters contained in the object
## S3 method for class 'hmmspec' print(x, ...)## S3 method for class 'hmmspec' print(x, ...)
x |
An object of type |
... |
further arguments passed to or from other methods. |
Jared O'Connell [email protected]
Prints the parameters contained in the object
## S3 method for class 'hsmmspec' print(x, ...)## S3 method for class 'hsmmspec' print(x, ...)
x |
An object of type |
... |
further arguments passed to or from other methods. |
Jared O'Connell [email protected]
This is an auxilliary data set to the cows data set containing times of artificial insemination for respective cows.
Only the day of insemination was recorded so time of day is always midday.
reproaireproai
reproai is a dataframe with 12 rows and id being the cow's id and days.from.calving recording the number of days from calving when insemination occurred.
Danish Cattle Research Centre
Peters, A. and Ball, P. (1995), "Reproduction in Cattle," 2nd ed.
This data set contains hourly observations on progesterone and an activity index at hourly intervals since calving on seven dairy cows.
reprocowsreprocows
reprocows is a data frame containing 13040 rows.
id is the cow ID, progesterone is a measurement of the
hormone in ng/L taken from a milk sample, activity is a relative measure
of activity calculated from a pedometer.
There are a large number of missing values as progesterone is measured only at milking time (and at a farm manager's discretion). Missing values in activity occur due to hardware problems can occur with pedometers.
Danish Cattle Research Centre
Peters, A. and Ball, P. (1995), "Reproduction in Cattle," 2nd ed.
data(reprocows) data(reproai) data(reproppa) tm = 1600 J <- 3 init <- c(1,0,0) trans <- matrix(c(0,0,0,1,0,1,0,1,0),nrow=J) emis <- list(mu=c(0,2.5,0),sigma=c(1,1,1)) N <- as.numeric(table(reprocows$id)) train <- list(x=reprocows$activity,N=N) class(train) <- "hsmm.data" tmp <- gammafit(reproppa * 24) M <- max(N) d <- cbind(dgamma(1:M,shape=tmp$shape,scale=tmp$scale), # ppa sojourn directly estimated from ppa data set dunif(1:M,4,30), # oestrus between 4 and 30 hours dunif(1:M,15*24,40*24)) #not-oestrus between 15 and 40 days startval <- hsmmspec(init,trans,parms.emission=emis,list(d=d,type='gamma'), dens.emission=dnorm.hsmm) #not run (takes some time) #h.activity <- hsmmfit(train,startval,mstep=mstep.norm,maxit=10,M=M,lock.transition=TRUE)data(reprocows) data(reproai) data(reproppa) tm = 1600 J <- 3 init <- c(1,0,0) trans <- matrix(c(0,0,0,1,0,1,0,1,0),nrow=J) emis <- list(mu=c(0,2.5,0),sigma=c(1,1,1)) N <- as.numeric(table(reprocows$id)) train <- list(x=reprocows$activity,N=N) class(train) <- "hsmm.data" tmp <- gammafit(reproppa * 24) M <- max(N) d <- cbind(dgamma(1:M,shape=tmp$shape,scale=tmp$scale), # ppa sojourn directly estimated from ppa data set dunif(1:M,4,30), # oestrus between 4 and 30 hours dunif(1:M,15*24,40*24)) #not-oestrus between 15 and 40 days startval <- hsmmspec(init,trans,parms.emission=emis,list(d=d,type='gamma'), dens.emission=dnorm.hsmm) #not run (takes some time) #h.activity <- hsmmfit(train,startval,mstep=mstep.norm,maxit=10,M=M,lock.transition=TRUE)
This data set contains the observed length of post-partum anoestrus (in days) for 73 dairy cattle.
reproppareproppa
reproppa a vector containing 73 integers.
Danish Cattle Research Centre
Peters, A. and Ball, P. (1995), "Reproduction in Cattle," 2nd ed.
This generates values from a multivariate normal distributed emission
state j given parameters in model.
rmvnorm.hsmm(j, model)rmvnorm.hsmm(j, model)
j |
An integer representing the state |
model |
A |
This is essentially a wrapper for rnorm. Users may build functions with the same
arguments and return values so they can use their own custom emission distributions.
A single value from the emission distribution.
Jared O'Connell [email protected]
J<-2 initial <- rep(1/J,J) P <- matrix(c(.3,.5,.7,.5),nrow=J) b <- list(mu=list(c(-3,0),c(1,2)),sigma=list(diag(2),matrix(c(4,2,2,3), ncol=2))) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dmvnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rmvnorm.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.mvnorm)J<-2 initial <- rep(1/J,J) P <- matrix(c(.3,.5,.7,.5),nrow=J) b <- list(mu=list(c(-3,0),c(1,2)),sigma=list(diag(2),matrix(c(4,2,2,3), ncol=2))) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dmvnorm.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rmvnorm.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.mvnorm)
This generates values from a normally distributed emission state j given parameters in model.
rnorm.hsmm(j, model)rnorm.hsmm(j, model)
j |
An integer representing the state |
model |
A |
This is essentially a wrapper for rnorm. Users may build functions with the same
arguments and return values so they can use their own custom emission distributions.
A single value from the emission distribution.
Jared O'Connell [email protected]
This generates values from a Poisson distributed emission state j given parameters in model.
rpois.hsmm(j, model)rpois.hsmm(j, model)
j |
An integer representing the state |
model |
A |
This is essentially a wrapper for rpois. Users may build functions with the same
arguments and return values so they can use their own custom emission distributions.
A single value from the emission distribution.
Jared O'Connell [email protected]
J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(lambda=c(1,3,6)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dpois.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rpois.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.pois)J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(lambda=c(1,3,6)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dpois.hsmm) model train <- simulate(model, nsim=300, seed=1234, rand.emis=rpois.hsmm) plot(train,xlim=c(0,100)) h1 = hmmfit(train,model,mstep=mstep.pois)
Simulates a Markov chain
sim.mc(init, transition, N)sim.mc(init, transition, N)
init |
The distribution of states at the first time step |
transition |
The transition probability matrix of the Markov chain |
N |
The number of observations to simulate |
A vector of integers representing the state sequence.
Jared O'Connell [email protected]
p <- matrix(c(.1,.3,.6,rep(1/3,3),0,.5,.5),ncol=3,byrow=TRUE) init <- rep(1/3,3) sim.mc(init,p,10)p <- matrix(c(.1,.3,.6,rep(1/3,3),0,.5,.5),ncol=3,byrow=TRUE) init <- rep(1/3,3) sim.mc(init,p,10)
Simulates data from a hidden Markov model
## S3 method for class 'hmmspec' simulate(object, nsim, seed = NULL, rand.emission=NULL,...)## S3 method for class 'hmmspec' simulate(object, nsim, seed = NULL, rand.emission=NULL,...)
object |
A |
nsim |
An integer or vector of integers (for multiple sequences) specifying the length of the sequence(s) |
seed |
|
rand.emission |
The function used to generate observations from the emission distribution |
... |
further arguments passed to or from other methods. |
If nsim is a single integer then a HMM of that length is
produced. If nsim is a vector of integers, then
length(nsim) sequences are generated with respective lengths.
An object of class hmmdata
x |
A vector of length |
s |
A vector of length |
N |
A vector of the length of each observation sequence (used to segment x and s) |
Jared O'Connell [email protected]
Rabiner, L. (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE, 77, 257-286.
hmmspec, link{predict.hmm}
J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dnorm.hsmm) train <- simulate(model, nsim=100, seed=1234, rand.emis=rnorm.hsmm) plot(train)J<-3 initial <- rep(1/J,J) P <- matrix(c(.8,.5,.1,0.05,.2,.5,.15,.3,.4),nrow=J) b <- list(mu=c(-3,0,2),sigma=c(2,1,.5)) model <- hmmspec(init=initial, trans=P, parms.emission=b,dens.emission=dnorm.hsmm) train <- simulate(model, nsim=100, seed=1234, rand.emis=rnorm.hsmm) plot(train)
Simulates values for a specified hidden semi-Markov model
## S3 method for class 'hsmmspec' simulate(object, nsim, seed = NULL,rand.emission=NULL,...)## S3 method for class 'hsmmspec' simulate(object, nsim, seed = NULL,rand.emission=NULL,...)
object |
A |
nsim |
An integer or vector of integers (for multiple sequences) specifying the number of states to generate per sequence |
seed |
|
rand.emission |
The function used to generate observations from the emission distribution |
... |
further arguments passed to or from other methods. |
If nsim is a single integer then a HSMM of that length is produced. If nsim is a vector of integers, then length(nsim) sequences are generated with respective lengths. Note thate length is the number of states visited, each state will have a sojourn time typically >1 so the vector will be longer than nsim
An object of class hmmdata
x |
A vector of length |
s |
A vector of length |
N |
A vector of the length of each observation sequence (used to segment x and s) |
Jared O'Connell [email protected]
Guedon, Y. (2003), Estimating hidden semi-Markov chains from discrete sequences, Journal of Computational and Graphical Statistics, Volume 12, Number 3, page 604-639 - 2003
hsmmfit,
hsmmspec,
predict.hsmm
J <- 3 init <- c(0,0,1) P <- matrix(c(0,.1,.4,.5,0,.6,.5,.9,0),nrow=J) B <- list(mu=c(10,15,20),sigma=c(2,1,1.5)) d <- list(lambda=c(10,30,60),shift=c(10,100,30),type='poisson') model <- hsmmspec(init,P,parms.emission=B,sojourn=d,dens.emission=dnorm.hsmm) train <- simulate(model,rand.emis=rnorm.hsmm,nsim=100,seed=123456) plot(train,xlim=c(0,400))J <- 3 init <- c(0,0,1) P <- matrix(c(0,.1,.4,.5,0,.6,.5,.9,0),nrow=J) B <- list(mu=c(10,15,20),sigma=c(2,1,1.5)) d <- list(lambda=c(10,30,60),shift=c(10,100,30),type='poisson') model <- hsmmspec(init,P,parms.emission=B,sojourn=d,dens.emission=dnorm.hsmm) train <- simulate(model,rand.emis=rnorm.hsmm,nsim=100,seed=123456) plot(train,xlim=c(0,400))
The smooth.discrete() function provides a simple smoothing of a time series of discrete values measured at equidistant times. Under the hood of smooth.discrete() is a hidden Markov model.
smooth.discrete(y, init = NULL, trans = NULL, parms.emission = 0.5, method = "viterbi", details = 0, ...)smooth.discrete(y, init = NULL, trans = NULL, parms.emission = 0.5, method = "viterbi", details = 0, ...)
y |
A numeric vector |
init |
Initial distribution (by default derived from data; see the vignette for details) |
trans |
Transition matrix (by default derived from data; see the vignette for details) |
parms.emission |
Matrix describing the conditional probabilities of the observed states given the latent states. (See the vignette for details). |
method |
Either "viterbi" or "smoothed". The viterbi method gives the jointly most likely sequence; the smoothed method gives the sequence of individually most likely states. |
details |
Controlling the amount of information printed. |
... |
Further arguments passed on to the "hmmfit" function. |
The parameters are estimated using the Baum-Welch algorithm (a special case of the EM-algorithm).
A list with the following components:
s |
The "smoothed" states |
model |
The underlying hmm (hidden Markov model) object |
data |
The data |
initial |
The initial parameters |
S<c3><b8>ren H<c3><b8>jsgaard <sorenh at math.aau.dk>
## Please see the vignette## Please see the vignette
Prints the estimated parameters of a hmm object
## S3 method for class 'hmm' summary(object, ...)## S3 method for class 'hmm' summary(object, ...)
object |
A |
... |
further arguments passed to or from other methods. |
An object of class 'summary.hmm'
Jared O'Connell [email protected]
Returns a summary object for a hsmm object
## S3 method for class 'hsmm' summary(object, ...)## S3 method for class 'hsmm' summary(object, ...)
object |
An object of type |
... |
further arguments passed to or from other methods. |
Jared O'Connell [email protected]