Modelling covariance matrices by the trigonometric separation strategy with application to hidden Markov models

Abstract

Bayesian inference on the covariance matrix is usually performed after placing an inverse-Wishart or a multivariate Jeffreys as a prior density, but both of them, for different reasons, present some drawbacks. As an alternative, the covariance matrix can be modelled by separating out the standard deviations and the correlations. This separation strategy takes advantage of the fact that usually it is more straightforward and flexible to set priors on the standard deviations and the correlations rather than on the covariance matrix. On the other hand, the priors must preserve the positive definiteness of the correlation matrix. This can be obtained by considering the Cholesky decomposition of the correlation matrix, whose entries are reparameterized using trigonometric functions. The efficiency of the trigonometric separation strategy (TSS) is shown through an application to hidden Markov models (HMMs), with conditional distributions multivariate normal. In the case of an unknown number of hidden states, estimation is conducted using a reversible jump Markov chain Monte Carlo algorithm based on the split-and-combine and birth-and-death moves whose design is straightforward because of the use of the TSS. Finally, an example in remote sensing is described, where a HMM containing the TSS is used for the segmentation of a multi-colour satellite image.

Appendices

Appendix A

The vector $\left( \mu ,S,\alpha ,\Omega ,m\right) ^{\prime }$ is estimated in a single RJMCMC run only, by using the values of $\mu$, S, $\alpha$, and $\Omega$ collected in those sweeps in which the number of hidden states is equal to the most frequent value of m. Nevertheless, to make the understanding of the algorithm easier, we first present the details of the MCMC algorithm for a given number m of states and then describe the RJMCMC when m is unknown.

Known number of states

In any sweep of the MCMC algorithm, the parameters $\mu$, S, $\alpha$, $\Omega$ are accepted or rejected, after generating their elements by random walk moves. Random walks can update real parameters; thus, positive quantities $\omega _{j,i}$, standard deviations $\sigma _{i,h}$, and angles $\alpha _{i,k,l}$ are mapped onto the real line through the logarithmic transformation. As the angles belong to the interval $( A_{i,k,l};B_{i,k,l})$, we consider the natural logarithm of $[ \alpha _{i,k,l}-A_{i,k,l}] /[ B_{i,k,l}-\alpha _{i,k,l}]$. The values $\mu _{i,h}$, $\ln \sigma _{i,h}$, $\ln [ \alpha _{i,k,l}-A_{i,k,l}] /[ B_{i,k,l}-\alpha _{i,k,l}]$, and $\ln \omega _{j,i}$, for any $i;j=1,\ldots ,m$, $h=1,\ldots ,p$, $k=1,\ldots ,p-1$, for any $l=k+1,\ldots ,p$, belong to the interval $\left( -\,\infty ;+\,\infty \right)$ and can be generated by the following random walk proposals:

$$\begin{aligned}&\mu _{i,h}=\mu _{i,h}^{(old)}+U_{M} \\&\ln \sigma _{i,h}=\ln \sigma _{i,h}^{(old)}+U_{\Sigma } \\&\ln \left[ \alpha _{i,k,l}-A_{i,k,l}\right] /\left[ B_{i,k,l}-\alpha _{i,k,l} \right] =\ln \left[ \alpha _{i,k,l}^{(old)}-A_{i,k,l}\right] /\left[ B_{i,k,l}-\alpha _{i,k,l}^{(old)}\right] +U_{A} \\&\ln \omega _{j,i}=\ln \omega _{j,i}^{(old)}+U_{\Omega }, \end{aligned}$$

where $U_{\Psi }\sim \mathcal {N}\left( 0;\sigma _{\Psi }^{2}\right)$, with $\Psi \in \left\{ M;\Sigma ;A;\Omega \right\}$.

The proposals $\mu$, S, $\alpha$, $\Omega$ are accepted if $u_{\Psi }\le \min \left\{ 1;A_{\Psi }\right\}$, where $u_{\Psi }$ is a random number generated from the uniform distribution $\mathcal {U}\left( 0;1\right)$ and $A_{\Psi }$ are the acceptance ratios ($\Psi \in \left\{ M;\Sigma ;A;\Omega \right\}$), i.e.

$$\begin{aligned} A_{\Psi }=\text {likelihood ratio }\times \text { prior ratio }\times \text { ratio of the products of the Jacobians,} \end{aligned}$$

where the likelihood ratio is

$$\begin{aligned} p(y^{T}\mid \mu ,S,\alpha ,\Omega ,m) /p(y^{T}\mid \mu ^{(old)},S ^{(old)},\alpha ^{(old)},\Omega ^{(old)},m) , \end{aligned}$$

the prior ratios, respectively, are

$$\begin{aligned} \begin{array}{cccc} p\left( \mu \right) /p\left( \mu ^{(old)}\right) ;&p\left( S\right) /p\left( S^{(old)}\right) ;&p\left( \alpha \right) /p\left( \alpha ^{(old)}\right) ;&p\left( \Omega \right) /p\left( \Omega ^{(old)}\right) , \end{array} \end{aligned}$$

and the ratio of the products of the Jacobian of the logarithmic transformations of the $\sigma _{i,h}$, the $\alpha _{i,k,l}$, and the $\omega _{j,i}$, respectively, are

$$\begin{aligned} \overset{m}{\underset{i=1}{\prod }}\overset{p}{\underset{h=1}{\prod }} \sigma _{i,h}\bigg / \overset{m}{\underset{i=1}{\prod }}\overset{p}{ \underset{h=1}{\prod }}\sigma _{i,h}^{(old)} , \frac{{{\prod }_{i=1}^m}{{\prod }_{k=1}^p} {{\prod }_{l=1}^p}\left( B_{i,k,l}-\alpha _{i,k,l}\right) \left( \alpha _{i,k,l}-A_{i,k,l}\right) }{{{\prod }_{i=1}^m} {{\prod }_{k=1}^p}{{\prod }_{l=1}^p}\left( B_{i,k,l}-\alpha _{i,k,l}^{(old)}\right) \left( \alpha _{i,k,l}^{(old)}-A_{i,k,l}\right) }, \end{aligned}$$

and

$$\begin{aligned} \overset{m}{\underset{i=1}{\prod }}\overset{m}{\underset{j=1}{\prod }}\omega _{j,i}\bigg / \overset{m}{\underset{i=1}{\prod }}\overset{m}{\underset{j=1}{ \prod }}\omega _{j,i}^{(old)}. \end{aligned}$$

At the end of each iteration, the MCMC sample is post-processed as in Marin et al. (2005). Assume that at any iteration k of the MCMC algorithm, we store the values $\{\mu ^{(k)},S^{(k)}, \alpha ^{(k)},\Omega ^{(k)}\}$. Let H be the class of the m! permutations $\eta _{j}$ of the labels ($\eta _{j}\in H$, for any $j=1,\ldots ,m!$), so that $\eta _{j}\ (\mu ^{(k)},S ^{(k)},\alpha ^{(k)},\Omega ^{(k)})$ be some permutation of the parameters obtained at the k-th iteration, by which the means, the standard deviations, the angles, and the rows and the columns of the transition matrix assume a new order.

After the burn-in, if a sample of size N ($k=1,\ldots ,N$) is simulated, the post-processing algorithm works as follows:

(i):

compute the posterior mode $\left\{ \mu ^{*}, S^{*},\alpha ^{*},\Omega ^{*}\right\}$, such that

$$\begin{aligned} \left\{ \mu ^{*},S^{*},\alpha ^{*}, \Omega ^{*}\right\} =\arg \underset{k=1,\ldots ,N}{\max } p(\mu ^{(k)},S^{(k)},\alpha ^{(k)},\Omega ^{(k)}\mid y^{T},m) \end{aligned}$$

(ii):

for any $k=1,\ldots ,N$, compute $\eta ^{*}$ such that

$$\begin{aligned} \eta ^{*}=\arg \underset{\eta _{j}\in H}{\min }\left\| \eta _{j}(\mu ^{(k)},S^{(k)},\alpha ^{(k)}, \Omega ^{(k)}) -(\mu ^{*},S^{*},\alpha ^{*},\Omega ^{*}) \right\| \end{aligned}$$

and place

$$\begin{aligned} (\mu ^{(k)},S^{(k)},\alpha ^{(k)},\Omega ^{(k)}) =\eta ^{*}(\mu ^{(k)},S ^{(k)},\alpha ^{(k)},\Omega ^{(k)}) . \end{aligned}$$

In step (ii), for any entry of the MCMC sample, we first compute the Euclidean norm between any permuted vector of parameters and the posterior mode; then, we select that special reordered vector which is the nearest to the posterior mode. Therefore, the label switching problem is circumvented without selecting any artificial identifiability constraint.

Finally, the hidden sequence of the states is reconstructed. Each state is the maximizer of the current smoothed probabilities: after obtaining parameter estimates, we can compute backwards the smoothed probabilities of the states (Kim 1993), that is the probabilities of any state, at any time, given all observations and the estimates of the parameters, i.e. $\widehat{ \mu }$, $\widehat{\Sigma }$, $\widehat{\Gamma }$:

$$\begin{aligned} x_{t}=\arg \max _{j}P( X_{t}=j\mid y^{T},\widehat{\mu } ,\widehat{\Sigma },\widehat{\Gamma },m) , \end{aligned}$$

with $t=1,\ldots ,T$, where

$$\begin{aligned}&P( X_{t}=j\mid y^{T},\widehat{\mu },\widehat{ \Sigma },\widehat{\Gamma },m) \nonumber \\&\quad =P( X_{t}=j\mid y ^{t},\widehat{\mu },\widehat{\Sigma },\widehat{ \Gamma },m) \overset{m}{\underset{i=1}{\sum }}\frac{\gamma _{j,i} \text { }P( X_{t+1}=i\mid y^{T},\widehat{\mu }, \widehat{\Sigma },\widehat{\Gamma },m) }{P( X_{t+1}=i\mid y^{t},\widehat{\mu },\widehat{\Sigma },\widehat{\Gamma },m) }, \end{aligned}$$

(3)

for any $t=T-1,\ldots ,1$ and any $j=1,\ldots ,m$, starting from $P( X_{T}\mid y^{T},\widehat{\mu },\widehat{\Sigma }, \widehat{\Gamma },m)$, which can be obtained by means of the filtered probabilities (2).

Unknown number of states

Our RJMCMC algorithm is based on three main moves, which allow changes in the number of hidden states:

[i]:

update the parameters as described in the previous subsection;

[ii]:

split one state of the MVN-HMM into two or merge two states into one;

[iii]:

give birth or death to a state.

In move [ii], the split is randomly chosen with probability $b_{m}= \mathbb {I}(m=1)+0.5\cdot \mathbb {I}(2\le m<m_{\max })$, whereas the combine is randomly chosen with probability $d_{m}=1-b_{m}$.

In the combine move, two adjacent states, e.g. $i_{1}$ and $i_{2}=i_{1}+1$, are randomly selected and combined in state $i^{*}$, reducing by one the number of hidden states; the corresponding parameters are combined as follows:

$$\begin{aligned} \begin{array}{ll} \mu _{i^{*},h}=( \mu _{i_{1},h}+\mu _{i_{2},h}) /2 &{}\quad \text { for any }h=1,\ldots ,p \\ \sigma ^{2}_{i^{*},h}=( \sigma ^{2}_{i_{1},h}\cdot \sigma ^{2} _{i_{2},h}) ^{1/2} &{}\quad \text {for any }h=1,\ldots ,p \\ \alpha _{i^{*},k,l}=\alpha _{i_{1},k,l}+\alpha _{i_{2},k,l} &{}\quad \text {for any }k=1,\ldots ,p-1\text { and any }l=k+1,\ldots ,p \\ \omega _{i,i^{*}}=\omega _{i,i_{1}}+\omega _{i,i_{2}} &{}\quad \text {for any } i\ne i^{*} \\ \omega _{i^{*},j}=( \omega _{i_{1},j}\cdot \omega _{i_{2},j}) ^{1/2} &{}\quad \text {for any }j\ne i^{*} \\ \omega _{i^{*},i^{*}}=( \omega _{i_{1},i_{1}}\cdot \omega _{i_{2},i_{1}}) ^{1/2}+( \omega _{i_{1},i_{2}}\cdot \omega _{i_{2},i_{2}}) ^{1/2} &{} \end{array} \end{aligned}$$

(4)

In the split move, a state $i^{*}$ is picked at random and split in the two adjacent states $i_{1}$ and $i_{2}$; the corresponding parameters are split as follows, respecting the six equalities in (4). First, we generate the following $p(p+3)/2+2m+1$ random values:

• $u_{1,h}$ from $\mathcal {N}\left( 0;0.5\right)$, for any $h=1,\ldots ,p$;

• $u_{2,h}$ from $\mathcal {G}\left( 1;5\right)$, for any $h=1,\ldots ,p$;

• $u_{3,k,l}$ from $\mathcal {U}\left( 0;1\right)$, for any $k=1,\ldots ,p-1$ and any $l=k+1,\ldots ,p$;

• $v_{i}$ from $\mathcal {U}\left( 0;1\right)$, for any $i\ne i^{*}$;

• $w_{j}$ from $\mathcal {G}\left( 1;5\right)$, for any $j\ne i^{*}$;

• $\rho$ from $\mathcal {U}\left( 0;1\right)$;

• $\tau _{1}$ and $\tau _{2}$ from $\mathcal {G}\left( 1;5\right)$.

Then, we set:

$$\begin{aligned} \begin{array}{ll} \mu _{i_{1},h}=\mu _{i^{*},h}-\sigma _{i^{*},h,h}\cdot u_{1,h} &{}\quad \mu _{i_{2},h}=\mu _{i^{*},h}+\sigma _{i^{*},h,h}\cdot u_{1,h} \\ \sigma ^{2}_{i_{1},h}=\sigma ^{2}_{i^{*},h}\cdot u_{2,h} &{}\quad \sigma ^{2}_{i_{2},h}=\sigma ^{2}_{i^{*},h}/u_{2,h} \\ \alpha _{i_{1},k,l}=\alpha _{i^{*},k,l}\cdot u_{3,k,l} &{}\quad \alpha _{i_{2},k,l}=\alpha _{i^{*},k,l}\cdot \left( 1-u_{3,k,l}\right) \\ \omega _{i,i_{1}}=\omega _{i,i^{*}}\cdot v_{i} &{}\quad \omega _{i,i_{2}}=\omega _{i,i^{*}}\cdot \left( 1-v_{i}\right) \\ \omega _{i_{1},j}=\omega _{i^{*},j}\cdot w_{j} &{}\quad \omega _{i_{2},j}=\omega _{i^{*},j}/w_{j} \\ \omega _{i_{1},i_{1}}=\omega _{i^{*},i^{*}}\cdot \rho \cdot \tau _{1} &{}\quad \omega _{i_{2},i_{1}}=\omega _{i^{*},i^{*}}\cdot \rho /\tau _{1} \\ \omega _{i_{1},i_{2}}=\omega _{i^{*},i^{*}}\cdot \left( 1-\rho \right) \cdot \tau _{2} &{}\quad \omega _{i_{2},i_{2}}=\omega _{i^{*},i^{*}}\cdot \left( 1-\rho \right) /\tau _{2} \end{array} \end{aligned}$$

The split move is accepted with probability $\min \left\{ 1;A\right\}$, while the combine move is accepted with probability $\min \left\{ 1;A^{-1}\right\}$. Let the tilde mark the parameters in the model with $m+1$ states, with respect to those entering in the model with m states; the analytic expression of A is

$$\begin{aligned}&\frac{p\left( y^{T}\mid \widetilde{\mu },\widetilde{S}, \widetilde{\alpha },\widetilde{\Omega },m+1\right) }{ p\left( y^{T}\mid \mu ,S,\alpha ,\Omega ,m\right) }\cdot \frac{p(m+1)}{p(m)}\cdot \frac{p\left( \mu \right) \cdot p\left( \widetilde{S}\right) \cdot p\left( \widetilde{ \alpha }\right) \cdot p\left( \widetilde{\Omega }\right) }{p\left( \widetilde{\mu }\right) \cdot p\left( S\right) \cdot p\left( \alpha \right) \cdot p\left( \Omega \right) }\cdot \frac{d_{m+1}/m}{b_{m}/m}\cdot \nonumber \\&\quad \cdot \frac{m+1}{{{\prod }_{h=1}^p}p( u_{1,h}) \cdot {{\prod }_{h=1}^p}p( u_{2,h}) \cdot {{\prod }_{p-1}^{k=1}}{{\prod }_{l=k+1}^p} p( u_{3,k,l}) \cdot {{\prod }_{i\ne i^*}}p( v_{i}) \cdot {{\prod }_{j\ne i^*}}p( w_{j}) \cdot p( \rho ) \cdot p( \tau _{1}) \cdot p( \tau _{2}) }\cdot |J|, \nonumber \\ \end{aligned}$$

(5)

where $p(m+1)/p(m)$ cancels out; $b_{m}/m$ is the probability of splitting the special state $i^{*}$, while $d_{m+1}/m$ is the probability of merging one of the m pairs $\left( i_{1};i_{2}\right)$ of adjacent states; factor $\left( m+1\right)$ is the ratio $\left( m+1\right) !/m!$, in which the factorials arise from the exchangeability assumption on the states; J is the Jacobian of the transformation from $(\omega _{i,i^{*}},v_{i},\omega _{i^{*},j},w_{j},\omega _{i^{*},i^{*}},\rho ,\tau _{1},\tau _{2},\mu _{i^{*},h},u_{1,h},$$\sigma ^{2}_{i^{*},h},u_{2,h},\alpha _{i^{*},k,l},u_{3,k,l})$ to $(\widetilde{\omega } _{i,i_{1}},\widetilde{\omega }_{i,i_{2}},\widetilde{\omega }_{i_{1},j}, \widetilde{\omega }_{i_{2},j},\widetilde{\omega }_{i_{1},i_{1}},\widetilde{ \omega }_{i_{1},i_{2}},\widetilde{\omega }_{i_{2},i_{1}},\widetilde{\omega } _{i_{2},i_{2}},\widetilde{\mu }_{i_{1}},\widetilde{\mu }_{i_{2}},\widetilde{ \sigma }^{2}_{i_{1},h},\widetilde{\sigma }^{2}_{i_{2},h},\widetilde{\alpha }_{i_{1},k,l}$, $\widetilde{\alpha }_{i_{2},k,l})$.

Note that the Jacobian can be decomposed into the product of five subdeterminants, i.e. $J_{1}$ for the transformation from $\left( \omega _{i,i^{*}},v_{i}\right)$ to $( \widetilde{\omega }_{i,i_{1}}, \widetilde{\omega }_{i,i_{2}})$, $J_{2}$ for the transformation from $( \omega _{i^{*},j},w_{j})$ to $( \widetilde{\omega } _{i_{1},j},\widetilde{\omega }_{i_{2},j})$, $J_{3}$ for the transformation from $( \omega _{i^{*},i^{*}},\rho ,\tau _{1},\tau _{2})$ to $( \widetilde{\omega }_{i_{1},i_{1}}, \widetilde{\omega }_{i_{1},i_{2}},\widetilde{\omega }_{i_{2},i_{1}}, \widetilde{\omega }_{i_{2},i_{2}})$, $J_{4}$ for the transformation from $( \mu _{i^{*},h},u_{1,h},\sigma ^{2}_{i^{*},h},u_{2})$ to $( \widetilde{\mu }_{i_{1,h}},\widetilde{\mu }_{i_{2},h}, \widetilde{\sigma }^{2}_{i_{1},h},\widetilde{\sigma }^{2}_{i_{2},h})$, $J_{5}$ for the transformation from $( \alpha _{i^{*},k,l},u_{3,k,l})$ to $( \widetilde{\alpha }_{i_{1},k,l}, \widetilde{\alpha }_{i_{2},k,l})$. Hence, we have $|J|=|J_{1}\cdot J_{2}\cdot J_{3}\cdot J_{4}\cdot J_{5}|$, where

$$\begin{aligned} J_{1}= & {} \underset{i\ne i^{*}}{\prod }\left( -\,\omega _{i,i^{*}}\right) \qquad J_{2}=\left( -\,2\right) ^{m-1}\cdot \underset{j\ne i^{*} }{\prod }\frac{\omega _{i^{*},j}}{w_{j}}\qquad J_{3}=-\,2\cdot \omega _{i^{*},i^{*}}^{3}\cdot \frac{\rho \cdot \left( 1-\rho \right) }{ \tau _{1}\cdot \tau _{2}} \\ J_{4}= & {} \left( -\,4\right) ^{p}\cdot \overset{p}{\underset{h=1}{\prod }}\frac{ \sigma _{i^{*},h}^{3}}{u_{2,h}}\qquad J_{5}=\overset{p-1}{\underset{ k=1}{\prod }}\underset{l=k+1}{\overset{p}{\prod }}\left( -\,\alpha _{i^{*},k,l}\right) . \end{aligned}$$

In move [iii], birth and death are chosen with probability $b_{m}$ and $d_{m}$, respectively. In a death move, a state is selected at random and then suppressed along with the corresponding parameters. In a birth move a new state $i^{*}$ is added to the previous m and the new parameters are drawn from their respective priors; the position of the new state is generated at random. The birth move is accepted with probability $\min \left\{ 1;A\right\}$, while the death move is accepted with probability $\min \left\{ 1;A^{-1}\right\}$; the analytic expression of A is

$$\begin{aligned}&\frac{p\left( y^{T}\mid \widetilde{\mu },\widetilde{S}, \widetilde{\alpha },\widetilde{\Omega },m+1\right) }{ p\left( y^{T}\mid \mu ,S,\alpha ,\Omega ,m\right) }\cdot \frac{p(m+1)}{p(m)}\cdot \frac{p\left( \mu \right) \cdot p\left( \widetilde{S}\right) \cdot p\left( \widetilde{ \alpha }\right) \cdot p\left( \widetilde{\Omega }\right) }{p\left( \widetilde{\mu }\right) \cdot p\left( S\right) \cdot p\left( \alpha \right) \cdot p\left( \Omega \right) }\\&\quad \cdot \frac{d_{m+1}/\left( m+1\right) }{b_{m}/\left( m+1\right) }\cdot \\&\quad \cdot \frac{m+1}{p( \mu _{i^{*}}) \cdot p( S_{i^{*}}) \cdot p( \alpha _{i^{*}}) \cdot {{\prod }_{i\ne i^*}}p( \omega _{i,i^{*}}) \cdot {{\prod }_{j\ne i^{*}}}p( \omega _{i^{*},j}) \cdot p( \omega _{i^{*},i^{*}}) }\cdot |J|, \end{aligned}$$

where $p(m+1)/p(m)$ cancels out; the ratio of the products of the prior densities multiplied by the reciprocal of the product of the densities of the new-born parameters is equal to 1; $b_{m}/\left( m+1\right)$ is the probability of giving birth to a new state in the special position $i^{*}$, while $d_{m+1}/\left( m+1\right)$ is the probability of killing a special state; factor $\left( m+1\right)$ has the same meaning as in (5); the Jacobian J is 1.

Appendix B

Estimates (posterior means) of the transition matrix:

$$\begin{aligned} \mathbf {\Gamma }=\left[ \begin{array}{ccc} 0.960 &{}\quad 0.029 &{}\quad 0.011 \\ 0.045 &{}\quad 0.921 &{}\quad 0.034 \\ 0.017 &{}\quad 0.032 &{}\quad 0.951 \end{array} \right] \end{aligned}$$

Estimates (posterior means) of the mean vectors:

$$\begin{aligned} \mathbf {\mu }_{1}= & {} \left( -\,1.292,-\,0.026,1.222\right) \quad \mathbf {\mu } _{2}=\left( -\,0.282,-\,0.064,0.467\right) \\ \mathbf {\mu }_{3}= & {} \left( -\,1.416,0.068,1.318\right) \end{aligned}$$

Estimates (posterior means) of the covariance matrices:

$$\begin{aligned} \mathbf {\Sigma }_{1}= & {} \left[ \begin{array}{ccc} 0.144 &{}\quad 0.273 &{}\quad 0.117 \\ 0.273 &{}\quad 1.309 &{}\quad 0.260 \\ 0.117 &{}\quad 0.260 &{}\quad 0.165 \end{array} \right] \quad \mathbf {\Sigma }_{2}=\left[ \begin{array}{ccc} 0.402 &{}\quad 0.249 &{}\quad 0.170 \\ 0.249 &{}\quad 0.836 &{}\quad 0.210 \\ 0.170 &{}\quad 0.210 &{}\quad 0.254 \end{array} \right] \\ \mathbf {\Sigma }_{3}= & {} \left[ \begin{array}{ccc} 0.195 &{}\quad 0.210 &{}\quad 0.120 \\ 0.210 &{}\quad 0.789 &{}\quad 0.190 \\ 0.120 &{}\quad 0.190 &{}\quad 0.105 \end{array} \right] \end{aligned}$$