\documentclass[11pt]{article} \usepackage{../../Math545} \usepackage{enumitem} \usepackage{bm} \usepackage{amsmath,amssymb} \pagestyle{plain} \renewcommand{\labelenumi}{\arabic{enumi}.} \renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}} \renewcommand{\labelenumiii}{\arabic{enumi}.\arabic{enumii}.\arabic{enumiii}} \setlength{\topmargin}{-0.75 in} \setlength{\textheight}{9.5 in} \parindent0in \newenvironment{packed_enum}{ \begin{enumerate} \setlength{\itemsep}{12pt} \setlength{\parskip}{0pt} \setlength{\parsep}{0pt} \setlength{\topsep}{0pt} }{\end{enumerate}} \def\myGt{\mathbf{G}_t} \def\myFt{\mathbf{F}_t} \def\myHonet{\mathbf{H}_{1t}} \def\myHtwot{\mathbf{H}_{2t}} \def\myZ{\mathbf{Z}} \def\evec{\mathbf{e}} \def\bP{\mathbf{P}} \def\bM{\mathbf{M}} \def\bK{\mathbf{K}} \def\bL{\mathbf{L}} \def\bv{\mathbf{v}} \def\ba{\mathbf{a}} \def\bOmega{{\boldsymbol{\Omega}}} \def\bSigma{{\boldsymbol{\Sigma}}} \def\beps{{\boldsymbol{\epsilon}}} \def\ctitle{MATH 545: Introduction to Time Series Analysis} \def\ttitle{Winter 2025} \def\stitle{Hidden Markov Models} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) @ \maketitle In a form of state space model, the state equation can be replaced by a Markov chain model: if $\{X_t\}$ can take values on $\{1,2,\ldots,S\}$ say, then for $t > 1$ \[ p(x_t|x_{1:(t-1)}) \equiv p(x_t|x_{t-1}) \] where \[ \Pr [X_t = l | X_{t-1} = k ] = p_{kl} \] The $S \times S$ transition matrix $\bP$, along with the (marginal) distribution of $X_1$ ($\Pr[X_1 = k], k=1,\ldots,S$) characterize the nature of the state evolution. Note that $\bP$ is a stochastic matrix whose rows sum to 1. Under mild conditions, the chain has stationary (or equilibrium) distribution, $\bpi$, defined as the solution of \[ \bpi = \bpi \bP \] where $\bpi$ represents the `long-run' behaviour of $\{X_t\}$. \medskip All of the previous types of calculation \begin{itemize} \item filtering $p(x_t|y_{1:t})$ \item prediction $p(x_t|y_{1:(t-1)})$ \item smoothing $p(x_t|y_{1:n})$ \item parameter estimation \end{itemize} are required. The main difference is that the filtering, prediction and smoothing distributions are \textbf{discrete}, and therefore somewhat more straightforward to compute. The observation variable $y_t$ can be either discrete or continuous. \bigskip \textbf{Example:} Suppose for $t=1,2,\ldots,n$ \[ Y_t = \mu(x_t) + Z_t \] and $\{X_t\}$ follows a Markov chain with transition matrix \[ \bP = \left[ \begin{array}{cc} 0.3 & 0.7\\ 0.4 & 0.6 \end{array} \right] \] and the distribution of $X_1$ is the stationarydistribution, $\bpi$, of the chain defined as the solution of \[ \bpi = \bpi \bP \] with $\bpi$ computed as \[ [1,1]^\transpose \: \bpi = \lim_{t \longrightarrow \infty} \bP^t = \left[ \begin{array}{cc} 4/11 & 7/11\\ 4/11 & 7/11 \end{array} \right] \] so $\bpi = (4/11,7/11)$. Suppose \[ \mu(x) = \left\{ \begin{array}{ll} \mu_1 \qquad & x = 1 \\[6pt] \mu_2 & x = 2 \end{array} \right. \] and $\{Z_t \} \sim \calN (0,1)$, IID. <>= n<-500 S<-2 P<-matrix(c(3,4,7,6)/10,2,2) mu<-c(2,-1) Ptmp<-P for(i in 1:100){Ptmp <- Ptmp %*% P} Ptmp Pi<-Ptmp[1,] c(4/11,7/11) set.seed(324) X<-rep(0,n) X[1]<-sample(c(1:S),size=1,prob=Pi) for(i in 2:n){X[i]<-sample(c(1:S),size=1,prob=P[X[i-1],])} X[1:20] Y<-rnorm(n,mu[X],1) par(mar=c(4,4,2,2)) cvec<-c('red','blue') plot(Y,pch=19,cex=0.5,col=cvec[X]) @ Observations from state 1 with $\mu=2$ plotted in red, observations from state 2 with $\mu=-1$ plotted in blue. \pagebreak \textbf{Key calculations:} \begin{itemize}[itemsep=0.05in] \item \emph{Filtering:} By Bayes theorem \begin{align*} p(x_t | y_{1:t}) & = p(x_t | y_t, y_{1:(t-1)}) \\[6pt] &= \dfrac{p(y_t | x_t,y_{1:(t-1)}) p(x_t |y_{1:(t-1)})}{p(y_{t}|y_{1:(t-1)})}\\[6pt] &= \dfrac{p(y_t | x_t) p(x_t |y_{1:(t-1)})}{\sum\limits_{x=1}^S p(y_t | x) p(x|y_{1:(t-1)})}\\[6pt] & \propto \ p(y_t | x_t) p(x_t |y_{1:(t-1)}) \end{align*} This is a discrete distribution on $\{1,\ldots,S\}$ that we can compute exactly (at least numerically) \medskip Note also that the denominator is the prediction distribution for $y_t$ given $y_{1:(t-1)}$ \[ p(y_{t}|y_{1:(t-1)}) = \sum\limits_{x=1}^S p(y_t | x) p(x|y_{1:(t-1)}) \] which can be used to compute the likelihood \[ p(y_1,\ldots,y_n) = p(y_1) p(y_2|y_1) p(y_3|y_1,y_2) \ldots p(y_n|y_1,y_2,\ldots,y_{n-1}) \] \item \emph{Prediction:} By the theorem of total probability \[ p(x_{t+1} | y_{1:t}) = \sum_{x=1}^S p(x_{t+1}|x) p(x|y_{1:t}) \] \begin{itemize} \item $p(x_{t+1}|x)$ is determined the $x$th row of $\bP$, for each $x=1,2,\ldots,S$ \item $p(x|y_{1:t})$ is the previous filtering distribution. \end{itemize} \item \emph{Smoothing:} By the theorem of total probability, we can compute \emph{backwards} $t=n,n-1,\ldots,1$, \[ p(x_t | y_{1:n}) = \sum_{x_{t+1}=1}^S p(x_{t}|x_{t+1},y_{1:n}) p(x_{t+1}|y_{1:n}) \] $p(x_{t+1}|y_{1:n})$ is the previous smoothing distribution. For the first term \begin{align*} p(x_{t}|x_{t+1},y_{1:n}) & \equiv p(x_{t}|x_{t+1},y_{1:t})\\[6pt] & = \frac{p(x_t|y_{1:t}) p(x_{t+1}|x_t,y_{1:t})}{p(x_{t+1}|y_{1:t})}\\[6pt] & = \frac{p(x_t|y_{1:t}) p(x_{t+1}|x_t)}{p(x_{t+1}|y_{1:t})} \end{align*} which can be computed from the filtering and prediction distributions. \end{itemize} <>= sig<-1 pX.filter<-pX.predict<-pX.smooth<-matrix(0,nrow=n,ncol=S) pX.predict[1,]<-Pi pX.filter[1,]<-dnorm(Y[1],mu,sig)*Pi Norm<-sum(pX.filter[1,]) pX.filter[1,]<-pX.filter[1,]/Norm LogLike<-log(Norm) Xpred<-Xsmooth<-rep(0,n) Xpred[1]<-which.max(pX.filter[1,]) for(i in 2:n){ pX.predict[i,]<- colSums(P*pX.filter[i-1,]) pX.filter[i,]<-dnorm(Y[i],mu,sig)*pX.predict[i,] Norm<-sum(pX.filter[i,]) pX.filter[i,]<-pX.filter[i,]/Norm LogLike<-LogLike+log(Norm) Xpred[i]<-which.max(pX.filter[i,]) } pX.smooth[n,]<-pX.filter[n,] Xsmooth[n]<-which.max(pX.smooth[n,]) for(i in (n-1):1){ pX.smooth[i,]<-pX.filter[i,] * (P %*% (pX.smooth[i+1,]/pX.predict[i+1,])) Xsmooth[i]<-which.max(pX.smooth[i,]) } #P[X=1] pmat<-cbind(Y[1:10],pX.filter[1:10,1],pX.smooth[1:10,1]) rownames(pmat)<-paste('t=',1:10) colnames(pmat)<-c('Y','Filter','Smooth') pmat table(X,Xpred) table(X,Xsmooth) @ The HMM structure is commonly implemented when the observation is also \emph{discrete}. Then the model is specified by \begin{itemize} \item \emph{transition probabilities} \[ p_{kl}=\Pr \left[ X_{t}=l|X_{t-1}=k\right] \] with initial distribution $\Pr[X_1=k]$, $k=1,\ldots,S$. \item \emph{emission probabilities} \[ e_{k} \left( m \right) =\Pr \left[ Y_{t}=m | X_{t}=k\right] \] \end{itemize} $e_k(m)$ can be specified for a finite set $(m \in \{1,2,\ldots,M\})$, or on the integers. \bigskip \textbf{The Viterbi Algorithm:} We seek the \emph{most probable} state path, that is, the most probable sequence of states in light of the data: \begin{eqnarray*} x_{1:n}^{\ast } &=&\arg \max_{x_{1:n}} p\left( x_{1:n},y_{1:n}\right) \\ && \\ &=&\arg \max_{x_{1:n}} p\left( x_{1:n}|y_{1:n}\right) \end{eqnarray*} \begin{itemize} \item $p\left( x_{1:n},y_{1:n} \right) $ is the joint probability or joint likelihood \item $p\left(x_{1:n}|y_{1:n}\right) $ is the posterior probability \end{itemize} The calculation is carried out using the \textit{Viterbi} algorithm. The recursion proceeds as follows \begin{enumerate}[label=(\roman*)] \item Initialize vector $v_{k}(0)=\Pr[X_1 = k]$, $k=1,\ldots,S$. \item Recursion: for $t=1,...,n$, compute \begin{align*} v_{l}\left( t\right) &=e_{l}\left( y_{t}\right) \max_{k}\left( v_{k}\left( t-1\right) p_{kl}\right) \\ point_{t}\left( l\right) &=\arg \max_{k}\left( v_{k}\left( t-1\right) p_{kl}\right) \end{align*} for each $l = 1,2,\ldots,S$ \item Termination \begin{align*} p\left( y_{1:n},x_{1:n}^{\ast }\right) &= \max_{k}\left( v_{k}(n)\right) \\ x_{n}^{\ast } &=\arg \max_{k}\left( v_{k} ( n )\right) \end{align*} \item Trace back the optimal path \[ x_{t-1}^{\ast }=point_{t}\left( x_{t}^{\ast }\right) \qquad t=n,n-1,\ldots,2. \] \end{enumerate} <>= V<-Point<-matrix(0,ncol=n,nrow=S) v0<-rep(1,S)/S for(l in 1:S){ V[l,1]<-dnorm(Y[1],mu[l],1)*max(P[l,] * v0) Point[l,1]<-which.max(P[l,] * v0) } Const<-sum(V[,1]) V[,1]<-V[,1]/Const Const.run<-log(Const) for(ti in 2:n){ vp<-V[,(ti-1)] for(l in 1:S){ V[l,ti]<-dnorm(Y[ti],mu[l],1)*max(P[l,] * vp) Point[l,ti]<-which.max(P[l,] * vp) } Const<-sum(V[,ti]) V[,ti]<-V[,ti]/Const Const.run<-Const.run+log(Const) } Xstar<-rep(0,n) Xstar[n]<-which.max(V[,n]) for(ti in n:2){ Xstar[ti-1]<-Point[Xstar[ti],ti] } par(mar=c(4,4,2,2)) plot(Y[1:50],pch=19,cex=0.5,xlab="t",ylab='Y',col=cvec[X[1:50]]) lines(c(1:50),mu[Xstar[1:50]],type="s") @ \pagebreak \textbf{The Forward Algorithm:} The probability of an observed sequence that yields the likelihood is given by \[ p\left( y_{1:n} \right) =\sum_{x_{1:n}} p\left( y_{1:n},x_{1:n}\right) \] However, the direct computation of this sum is not in general possible, so a recursive calculation is needed. Note first that for $1 \leq t \leq n$ \[ p\left( y_{1:t} \right) = \sum_{x_t=1}^S p(x_t,y_{1:t}) \] and we may focus on computing the terms in the sum \[ f_{x_t}(t) = p(x_t,y_{1:t}) \] for $x_t=1,2,\ldots,S$. Secondly, notice that \begin{align*} p(x_t,y_{1:t}) & = \sum_{x_{t-1}=1}^S p(x_{t-1},x_t,y_{1:t}) & \text{total probability} \\[6pt] & = \sum_{x_{t-1}=1}^S p(y_t|x_t,x_{t-1},y_{1:(t-1)}) p(x_t|x_{t-1},y_{1:(t-1)}) p(x_{t-1},y_{1:(t-1)}) & \text{chain rule}\\[6pt] & = \sum_{x_{t-1}=1}^S p(y_t|x_t) p(x_t|x_{t-1}) p(x_{t-1},y_{1:(t-1)}) & \text{conditional independence}\\[6pt] & = p(y_t|x_t) \sum_{x_{t-1}=1}^S p(x_t|x_{t-1}) p(x_{t-1},y_{1:(t-1)}) \end{align*} In this expression, $p(y_t|x_t)$ is the observation/emission probability, $p(x_t|x_{t-1})$ is the transition probability, and $p(x_{t-1},y_{1:(t-1)})$ is the joint probability from the previous recursion step which we denoted $f_{x_{t-1}}(t-1)$. Therefore, we simply need to track the values of $f_{x}(t)$ for all $x$ and $t$, with updates using the recursive formula. To set up the recursion, we need \[ f_{x_1}(1) = p(x_1,y_1) = p(y_1|x_1) p(x_1) \] which we can compute from the observation/emission probability and the probability distribution for $X_1$, which we have assumed to be $\bpi$. \begin{enumerate}[label=(\roman*)] \item Initialize vector $f_{k}(1)=e_k(y_1) \pi_k$ for $k=1,2\ldots,S$ \item Recursion: for $t=2,...,n$, and $l=1,2,\ldots,S$, \[ f_{l}\left( t\right) =e_{l}\left( y_{t}\right) \sum_{k=1}^S f_{k}\left( t-1\right) p_{kl} \] \item Final calculation: \[ p\left( y_{1:n} \right) =\sum_{k=1}^S f_{k} \left( n\right). \] \end{enumerate} Note: in this calculation, when $n$ is large, $p(y_{1:n})$ is likely to be very small, so it is preferable to work on the log scale -- however, we must compute the forward quantities on the original scale. We may scale by (for example) \[ m_t = \max_k f_k(t) \] at each step to prevent numerical underflow. <>= fmat<-matrix(0,nrow=n,ncol=S) logm<-rep(0,n) fmat[1,]<-dnorm(Y[1],mu,sig)*Pi logm[1]<-log(max(fmat[1,])) fmat[1,]<-fmat[1,]/max(fmat[1,]) for(i in 2:n){ fmat[i,]<-dnorm(Y[i],mu,sig)*(P %*% fmat[i-1,]) logm[i]<-log(max(fmat[i,])) fmat[i,]<-fmat[i,]/max(fmat[i,]) } log(sum(fmat[n,]))+sum(logm) LogLike @ We also seek the posterior (smoothing) probabilities \begin{align*} p\left( X_{t}=k|y_{1:n}\right) &=\frac{p\left( X_{t}=k,y_{1:n}\right) }{ p\left( y_{1:n}\right) }\\[6pt] &=\frac{p\left( y_{1},...,y_{t},X_{t}=k\right) p\left( y_{t+1},...,y_{n}|X_{t}=k\right) }{p\left( y_{1:n} \right) } \end{align*} for $k=1,2,\ldots,S$. Direct computation is not possible, so a recursive calculation is needed, after completing the Forward Algorithm. \medskip \textbf{The Backward Algorithm:} The backward algorithm works in reverse order through the data: \begin{enumerate}[label=(\roman*)] \item Initialize element $b_{k}\left( n\right) = f_{k}(n)$ for $k=1,\ldots,S$. \item Backwards recursion: for $t=n-1,...,1$ compute \[ b_{k}\left( t \right) =\sum_{l=1}^S p_{kl}e_{l}\left( y_{t+1}\right) b_{l}\left( t+1\right) \] for $k=1,2,\ldots,S$. \end{enumerate} <>= bmat<-matrix(0,nrow=n,ncol=S) logbm<-rep(0,n) bmat[n,]<-fmat[n,] logbm[n]<-log(max(bmat[n,])) bmat[n,]<-bmat[n,]/max(bmat[n,]) for(i in (n-1):1){ bmat[i,]<-dnorm(Y[i],mu,sig)*(t(P) %*% bmat[i+1,]) logbm[i]<-log(max(bmat[i,])) bmat[i,]<-bmat[i,]/max(bmat[i,]) } @ The forward and backward probabilities, $f_{k}(t)$ and $b_{k}(t),$ are in fact the probabilities \begin{align*} f_{k}(t) &= p\left( y_{1},...,y_{t},X_{t}=k\right) \\ b_{k}(t) &= p\left( y_{t+1},...,y_{n}|X_{t}=k\right) \end{align*} Thus we have \[ p\left( X_{t}=k,y_{1:n} \right) = f_{k}(t)b_{k}(t) \] and hence, for $k=1,\ldots,S$, \[ p\left( X_{t}=k|y_{1:n}\right) =\frac{f_{k}(t)b_{k}(t)}{p\left( y_{1:n} \right) } \] \textbf{Estimation:} The algorithms above are all expressed in terms of known transition and emission probabilities. These probabilities can be pre-specified, estimated from observed sequences, or estimated during the state sequence reconstruction. \bigskip If the state sequence is \emph{known}, let $n_{kl}$ be the number of transitions from state $k$ to state $l$, and let $r_{k}\left(b\right) $ be the number of times character $b$ is emitted from state $k$. Then we may estimate the transition and emission probabilities as follows \begin{equation} \widehat{p}_{kl}=\dfrac{n_{kl}}{\sum\limits_{l^{\prime }}n_{kl^{\prime }}} \qquad \widehat{e}_{k}\left( b \right) =\dfrac{r_{k}\left( b\right) }{\sum\limits_{b^{\prime }}r_{k}\left( b^{\prime }\right) } \label{abest} \end{equation} If the state sequence is \emph{not known}, we must reconstruct the state sequence and perform estimation concurrently. This is achieved using the \emph{Baum-Welch} algorithm, which is a special case of the \emph{Expectation-Maximization} (EM) algorithm. \medskip \textbf{The Baum-Welch Algorithm:} The Baum-Welch algorithm estimates $n_{kl}$ and $r_{k}\left(b\right) $ as the \emph{expected} number of transitions and emissions, which are computed using the forward and backward algorithms. The probability of a transition $k\rightarrow l$ $\ $at time $t$ is \[ p\left( X_{t}=k,X_{t+1}=l|y_{1:n}\right) =\frac{f_{k}\left( t\right) p_{kl}e_{l} \left( y_{t+1}\right) b_{l}\left( t+1\right) }{p\left( y_{1:n} \right) } \] Suppose that $N$ replicate observed sequences \[ y_{1:n}^1,\ldots,y_{1:n}^N \] are available. Then the expected number of times this transition will be observed (in light of the data) is \[ n_{kl} = \sum_{i=1}^N \left\{ \dfrac{ \sum\limits_{t}f_{k}^{i}\left( t\right) a_{kl}e_{l} ( y_{t+1}^{i} ) b_{l}^{i}\left( t+1\right) }{p ( y_{1:n}^i ) }\right\} \] \noindent In this formula \begin{itemize} \item $y_{1:n}^i$ is the $i^{th}$ observed sequence \item $f_{k}^{i}\left( t\right) $ is the forward probability from the forward algorithm implemented for sequence $i$ \item $b_{l}^{i}\left( t\right) $ is the backward probability from the backward algorithm implemented for sequence $i$ \end{itemize} \bigskip \noindent Similarly, \[ r_{k}(b)=\sum_{i=1}^N \left\{ \dfrac{\sum\limits_{t \: \in \: \mathcal{B}_{i}} f_{k}^{i}\left( t\right) b_{l}^{i} ( t )}{p\left( y_{1:n}^i \right) } \right\} \] where $\mathcal{B}_{i}$ is the set of indices of locations in replicate sequence\ $i$ where character $b$ is emitted. \begin{enumerate}[label=(\roman*),itemsep=0.05in] \item Initialize parameter values to some arbitrary values \item Recursion: Set $n_{kl}^{(0)},$ $r_{k}^{(0)}(b)$ to be some small count value, then for $i=1,...,N$: calculate $f_{k}^{i}(t)$ using the Forward Algorithm, and $b_{k}^{i}\left( t\right) $ using the Backward Algorithm, then set {\small \begin{align*} n_{kl}^{(i)} &= n_{kl}^{(i-1)}+\frac{1}{p\left( y_{1:n}^i \right) }\left\{ \sum_{t=1}^n f_{k}^{i}\left( t\right) p_{kl} e_{l}\left( y_{t+1}^{i}\right) b_{l}^{i}\left( t+1\right) \right\}\\[6pt] r_{k}^{(i)} ( b ) &= r_{k}^{(i-1)} (b) +\frac{1}{p\left( y_{1:n}^i \right) }\left\{ \sum_{t\in \mathcal{B}_{i}}f_{k}^{i} (t) b_{k}^{i}( t) \right\} \end{align*}} \item Estimate $p_{kl}$ and $e_{k}(b) $ using (\ref{abest}), and compute the likelihood \[ \prod\limits_{i=1}^N p(y_{1:n}^i) \] \item Terminate the recursion if the computed likelihood value is virtually unchanged compared to the previous iteration. \end{enumerate} Note that most of the above calculations carry through to the continuous emission case, that is if the emission probabilities \[ e_{k} ( m ) =\Pr [ Y_{t}= m |X_{t}=k ] \] are replaced by densities \[ e_{k} ( y ) = p ( y | X_{t}=k ) \] The forward and backward quantities become unnormalized joint densities \medskip To avoid underflow in the computation, the computations can either be computed on the log scale, or with continual rescaling of the calculation. For example, in the Viterbi algorithm, the computation of the $v_k(t)$ quantities will underflow quickly without rescaling, but rescaling (for example) by \[ c_t = \sum_{k=1}^S v_k(t) \] will prevent this. The constants $c_t$ must be stored for future use; if the forward and backward probabilities are computed after rescaling, then this rescaling must be undone when, for example, $p(y_{1:n})$ is computed. \bigskip \textbf{Direct estimation using filtering:} As for the Kalman Filter, we can use the filtering approach to compute the likelihood, and perform maximum likelihood estimation. <>= ml.HMM<-function(parv,Yv,Sv){ P<-matrix(parv[1:(Sv*(Sv-1))],nrow=Sv,byrow=T) P<-cbind(P,0) P<-exp(P) P<-t(apply(P,1,function(x){return(x/sum(x))})) mu<-parv[(Sv*(Sv-1)+1):(Sv*Sv)] sig<-exp(parv[(Sv*Sv)+1]) n<-length(Yv) Ptmp<-P;for(i in 1:100){Ptmp<-Ptmp%*%P} Pi<-Ptmp[1,] pXf<-pXp<-matrix(0,nrow=n,ncol=Sv) pXp[1,]<-Pi pXf[1,]<-dnorm(Yv[1],mu,sig)*Pi Norm<-sum(pXf[1,]) pXf[1,]<-pXf[1,]/Norm LogLike<-log(Norm) for(i in 2:n){ pXp[i,]<- colSums(P*pXf[i-1,]) pXf[i,]<-dnorm(Yv[i],mu,sig)*pXp[i,] Norm<-sum(pXf[i,]) pXf[i,]<-pXf[i,]/Norm LogLike<-LogLike+log(Norm) } return(LogLike) } S^2+1 xp<-c(0,0,2,0,0) xfit<-optim(xp,ml.HMM,Yv=Y,Sv=S,control=list(fnscale=-1),method='BFGS') par.max<-xfit$par P.hat<-matrix(par.max[1:(S*(S-1))],nrow=S,byrow=T) P.hat<-cbind(P.hat,0) P.hat<-exp(P.hat) P.hat<-t(apply(P.hat,1,function(x){return(x/sum(x))})) mu.hat<-par.max[(S*(S-1)+1):(S*S)] sig.hat<-exp(par.max[(S*S)+1]) P.hat mu.hat sig.hat @ \end{document}