% \VignetteIndexEntry{Smoothing discrete data (II)}
% \VignetteDepends{mhsmm}
% \VignetteKeyword{Hidden Markov Model}

\documentclass{article}
\usepackage[top=1in,left=1in,bottom=1in,right=1in]{geometry}
%\usepackage{a4wide}
\usepackage[T1]{fontenc}

\title{Smoothing discrete data (II)\\ -- using a hidden Markov model as
  implemented in the \texttt{mhsmm} package}

\author{S\o{}ren H\o{}jsgaard and Jared O'Connell}

\begin{document}


\SweaveOpts{prefix.string=fig/dhmm}
\setkeys{Gin}{width=0.5\textwidth} 
\renewenvironment{Schunk}{\linespread{.85}\small}{}

\maketitle

\parindent0pt\parskip5pt

@ 
<<print=F,echo=F>>=
options("width"=80)
@ %def 






We apply the \verb'mhsmm' package for a simple smoothing task. The
data pertain to classification of a cows eating behaviour over
time. The true eating status $S_t$ is in the vector \verb'cls0' where
'1' denotes not eating and '2' denotes eating. The time resolution is
one minute.  The observed variables $x_t$ in \verb'cls1' are actually
not observations per se but the result of a classification obtained by
a neural network (using \verb'nnet()' from the \verb'nnet'
package). 

See also the vignette ``Smoothing discrete data (I)'' for an
alternative approach to smoothing the data.


@ 
<<>>=
load("clsX.RData")
length(cls0)
length(cls1)
library(mhsmm)
@ %def 

@ 
<<clsfig1, fig=T, include=F>>=
plot(cls1[1:2000], type='l', ylim=c(.8,2))
addStates(cls0[1:2000])
@ %def 

\begin{figure}[h]
  \centering
  \includegraphics[width=\textwidth,height=6cm]{fig/dhmm-clsfig1}
  \caption{Observed and true eating states}
  \label{fig:clsfig1}
\end{figure}

A simple 'smoothing' of the observed states can be obtained as follows:

The density function for the emission distribution is
@ 
<<>>=
dpmf <- function(x,j,model) model$parms.emission$pmf[j,x]
@ %def 


An initial setting of the parameters is as follows:

@ 
<<>>=
J <- 2
init <- c(.5,.5)
P <- matrix(c(.9,.1,.1,.9),nrow=J)
B <- list(pmf=matrix(.1,ncol=J,nrow=J))
diag(B$pmf) <- .9
init.spec <- hmmspec(init,trans=P,parms.emission=B,dens.emission=dpmf)
init.spec
@ %def 

To fit the model we need to provide the function for the M--step of
the EM--algorithm:

@ 
<<>>=
mstep.pmf <- function(x,wt) {
  ans <- matrix(ncol=ncol(wt),nrow=ncol(wt))
  for(i in 1:ncol(wt))
    for(j in 1:ncol(wt))
      ans[i,j] <- sum(wt[which(x==j),i])/sum(wt[,i])
  list(pmf=ans)
}
@ %def 




For training the model we use the first 1000 cases
@ 
<<>>=
samp <- 1:2640
train <- list(s=cls0[samp], x=cls1[samp], N=length(cls0[samp]))
valid <- list(x=cls1[-samp], N=length(cls1[-samp]))
@ %def  

We fit the model with
@ 
<<>>=
hmm.obj <- hmmfit(train, init.spec,mstep=mstep.pmf)
summary(hmm.obj)
@ %def 


Two types of predictions can be made: Default is to use the Viterbi
algorithm for producing the jointly most likely sequence of states
given the observed data:

@ 
<<>>=
vit <- predict(hmm.obj, valid)
@ %def 

Alternatively we can get the individually most likely state sequence as:
@ 
<<>>=
smo <- predict(hmm.obj, valid, method="smoothed")
@ %def 

The prediction results are quite similar:
@ 
<<>>=
normtab <- function(tt) round(sweep(tt,1,rowSums(tt),"/"),2)
SS <- cls0[-samp]
XX <- cls1[-samp]
cls2.vit <- vit$s 
cls2.smo <- smo$s 
normtab(table(SS,XX))
normtab(table(SS,cls2.vit))
normtab(table(SS,cls2.smo))
@ %def 

@ 
<<clsfig2, fig=T, include=F>>=
show <- 1:2000
cls0b <- cls0[-samp]
cls1b <- cls1[-samp]
c(length(SS),length(cls2.vit),length(cls2.smo), length(cls0b), length(cls1b))
plot(cls1b[show], type='l', ylim=c(.8,2))
addStates(list(cls0b[show],cls2.vit[show], cls2.smo[show]))
@ %def 

\begin{figure}[h]
  \centering
  \includegraphics[width=\textwidth,height=6cm]{fig/dhmm-clsfig2}
  \caption{Observed and true eating states}
  \label{fig:clsfig1}
\end{figure}





% --------------------------------------------------------



% @ 
% <<>>=
% library(mhsmm)

% J <- 3
% init <- c(0,0,1)
% P <- matrix(c(0,.1,.4,.5,0,.6,.5,.9,0),nrow=J)
% d <- list(lambda=c(10,30,60),shift=c(10,100,30),type='poisson')
% @ %def 

% our emission distribution is just specified with a matrix where
% row is hidden state and column is probability of observing that state
% @ 
% <<>>=
% B <- list(pmf=matrix(.1,ncol=3,nrow=3))
% diag(B$pmf) <- .8
% @ %def 


% simulation function
% @ 
% <<>>=
% rnorm.pmf <- function(j,model) sample.int(nrow(model$emission$pmf),1,prob=model$emission$pmf[j,])
% @ %def 

% @ 
% <<>>=
% model <- hsmmspec(init,P,emission=B,sojourn=d,r=rnorm.pmf)
% train <- simulate(model,nsim=100,seed=123456)
% @ %def 



% @ 
% <<>>=
% mmm <- hmmspec(init=init, trans=P, emis=B,r=rnorm.pmf)
% ttt <- simulate(mmm,nsim=100,seed=123456)
% @ %def 

% @ 
% <<fig=T>>=
% plot(ttt)
% @ %def 


% plots look terrible...
% @ 
% <<>>=
% plot(train,xlim=c(0,600))
% @ %def 

% density function
% @ 
% <<>>=
% dpmf <- function(x,j,model) model$emission$pmf[j,x]
% @ %def 

% reestimation is straightforward

% @ 
% <<>>=
% mstep.pmf <- function(x,wt) {
%   ans <- matrix(ncol=ncol(wt),nrow=ncol(wt))
%   for(i in 1:ncol(wt))
%     for(j in 1:ncol(wt))
%       ans[i,j] <- sum(wt[which(x==j),i])/sum(wt[,i])
%   list(pmf=ans)
% }
% @ %def 

% @ 
% <<>>=
% hhh <- hmmfit(ttt, mmm,f=dpmf,mstep=mstep.pmf)
% yhat <- predict(hhh,ttt)
% @ %def 




% @ 
% <<>>=
% h.pmf <- hsmmfit(train,model,f=dpmf,mstep=mstep.pmf,M=300)
% yhat <- predict(h.pmf,train)
% mean(train$s!=yhat$s) #error rate
% table(train$s,yhat$s)
% @ %def 

% jjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjjj




% @ 
% <<>>=
% J <- 2
% init <- c(0,1)
% P <- matrix(c(.9,.1,.1,.9),nrow=J)
% B <- list(pmf=matrix(.1,ncol=J,nrow=J))
% diag(B$pmf) <- .8
% init.spec <- hmmspec(init,trans=P,emission=B)
% @ %def 

% @ 
% <<>>=
% tset <- 1:100
% train <- list(x=cls0[tset], N=length(cls0[tset]))
% valid <- list(x=cls0[-tset], N=length(cls0[-tset]))
% hhh <- hmmfit(train, init.spec,f=dpmf,mstep=mstep.pmf)

% yhat <- predict(hhh, valid)

% oo <- cls0[-tset]
% ee <- yhat$s 
% table(oo,ee) 
% @ %def  






\end{document}