Likelihood of a hidden Markov model does not increase in case of multidimensional data

I am implementing a HMM (Hidden Markov Model) for time series data to assess (define) state of the model at each time period (with continuous observations). To maximize the likelihood of the model (given an observed data, the probability of hidden state the model is in at each data point) I used expectation maximization algorithm (in case of hmm - Baum-Welch algorithm).

The problem is that in case of multidimensional data (the observation at each time is a vector), defined likelihood does not increase after each iteration. In the case of one-dimensional data when my X’s are numbers at each time point and mean and the sigma are also numbers I used cdf of Gaussian distribution (in ε-vicinity of data where ε is very small positive number). At the end, the algorithm works and at each step probabilities increase.

In my model I assumed that the conditional probability distribution of data is sampled either from an autoregressive process of order 1 or sampled from a Gaussian white noise distribution depending on which state is the model.

The steps of the algorithm are as follows:

1. Initialization of Emission, Transition and Pi matrices,

2. Calculation of forward and backward probabilities,

• Initialization $\alpha_i(1) = \pi_i \cdot p_i(\text{o}_1)$
• Induction $\alpha_i(t+1) = \big[ \sum_{j=1}^N \alpha_j(t)a_{ji} \big]p_i(\text{o}_{t+1})$
• Initialization $\beta_i(T)=1$
• Induction $\beta_i(t) = \sum_{j=1}^Na_{ij}p_j(\text{o}_{t+1})\beta_j(t+1) $

Using the forward and backward probabilities, the state probabilities are calculated as

2. Re-estimate Model Parameters

• Updating the transition probability by

• Updating the mean vectors and covariance matrices of the Gaussian by

I re-calculated the state probabilities using updated model parameters and re-estimated the model parameters using the new probabilities. I repeated these steps until the maximization of the likelihood has reached a desired convergence (i.e., when the improvement of the likelihood is negligible from one step to the next).

For one dimensional case I used scipy.stats.norm for Emission matrix calculation .

For multi-dimensional data I used scipy.stats.multivariate.normal with mean and covariance matrices.

Here is the calculation of Emission matrix:

from scipy.stats import multivariate_normal as mvn

distr = mvn(mean=mean, cov=covariance)
prob_t = abs(distr.cdf(obs + epsilion) - distr.cdf(obs - epsilion))

where epsilion is arbitrary very small number, obs is my data point (which is a vector) at time t and prob_t is conditional probability of my data being in particular state at time t.

Any idea why for multi-dimensional case my model does not work?