# Finding a mixture of 1st and 0'th order Markov models that is closest to an empirical distribution

I am interested in finding the distribution "$p^*$" closest to an empirical distribution $\hat{p}$ where $p^*$ is a mixture of first and zeroth order Markov models. That is, I want to find $$p^* = \arg\min_p \sum_{i,j} D\left(\,\hat{p}(j|i)\, \| \, p(j|i) \,\right)$$ subject to the following constraints

1. $p(j|i) = \alpha \cdot r(j) + (1 - \alpha) \cdot q(j|i)$
2. $r(j) \ge 0$ for all $j$
3. $\sum_j r(j) = 1$
4. $q(j|i) \ge 0$ for all $i$ and $j$
5. $\sum_j q(j|i) = 1$ for all $i$

where $\alpha$ is a mixture parameter in $[0,1]$ that is given and fixed. (The model is not identifiable if $\alpha$ is also a parameter to be learned.)

I know I can use EM here -- the missing data is the indicator variable for which Markov chain was invoked at a given time-step. But I am more interested in taking a direct optimization approach, particularly if this is a convex problem.

How might one derive the updates for for optimizing this model? Any assistance in this endeavor is greatly appreciated.

You can solve this problem by fitting only the more complex model, then transforming the resulting parameter estimates $q^*$ to recover $q$ and $r$. No iterative optimization needed!

Please forgjve me if jn thjs answer J have mjxed up $j$ and $i$. Jt's iust very djffjcult for me to keep them strajght.

Suppose your dataset consists of one observation of length $T$: $\{X_t\}_{n=0}^T$. (Below, I'll extend my answer to a case with multiple series.) I understand your problem to be that at each time $t$, there is a latent variable $z_t$ that is Bernoulli with known parameter $\alpha$. If $z_t=0$, then $x_{t}$ is drawn from $r$, and if $z_t=1$, then $x_{t}$ is drawn from $q(\cdot|x_{t-1})$. If $z$ is shared for all series, that's a very different situation, which this answer does not address. So, the contribution of epoch $t$ to the likelihood is this.

$$P(x_{t}|x_{t-1}) = \sum_z P(z_t, x_t |x_{t-1} ) = \sum_z P(z_t)P( x_t | x_{t-1}, z_t) = (1-\alpha) r(x_{t}) + \alpha q(x_{t}|x_{t-1})$$

This is much easier to look at if one defines $q^*(i|j)\equiv (1-\alpha) r(j) + \alpha q(j|i)$: it becomes just $q^*(x_{t}|x_{t-1})$. The complete-data likelihood (conditioning on $x_0$) is $\prod_{t=1}^{T} q^*(x_{t}|x_{t-1})$, and you could use the usual maximum likelihood estimates: $\hat q^*(j|i) \equiv \frac{N_{i\rightarrow j }}{N_{j}}$, where $N_{j}$ is the number of occurrences of $i$ and $N_{i\rightarrow j }$ is the number of $i$ to $j$ transitions. You can also add pseudocounts if you want to "smooth" things out; those could be interpreted as MAP estimates under suitable priors on $q^*$. (How that relates to priors on $q$ and $r$ would be interesting to investigate...)

Define the estimates for $r$ and $q$ (I'll call them $\hat r$ and $\hat q$) to be any numbers satisfying $\alpha\hat q(j|i) + (1-\alpha)\hat r(j) = \hat q^*(j|i)$. If the state space has cardinality $J$, you now face a linear system with $J^2 + J$ unknowns ($\hat q$'s and $\hat r$'s). The number of equations is tricky to count because some of the sum-to-one constraints are redundant, but my calculation is:

• J^2 equations from the definition I gave
• J equations from the sum-to-one constraints on $\hat q$
• 1 equation from the sum-to-one constraints on $\hat r$, which I am going to ignore because I suspect it's redundant.

If you vectorize $\hat q,\hat r$ as follows (dropping the hats, sorry):

$$v = [q(1|1), q(1|2),... q(1|J), j(1), q(2|1), ... j(2), ... q(J|1), ... r(J)]$$

and write the RHS as

$$w = [q^*(1|1), q^*(1|2), ...q^*(1|J), 1, ..., q^*(J|1), q^*(J|2), ...q^*(J|J), 1]$$,

then your system is $Av = w$, with $w$ known and $A$ known. $A$ is square and has a block diagonal structure with blocks of size $J+1$. Each block has the form

\begin{bmatrix} \alpha & 0 & 0 & 1-\alpha \\ 0 & \alpha & 0 & 1-\alpha \\ 0 & 0 & \alpha & 1-\alpha \\ 1 & 1 & 1 & 0 \\ \end{bmatrix} .

This is invertible, as can be seen by subtracting, from the bottom row, $\alpha^{-1}$ times each of the top three, to yield an upper triangular block with a nonzero diagonal. This implies that the parameters $q$ and $r$ are identifiable, and fittingly, it only holds when $\alpha$ is not zero or 1. (Of course, identifying $r$ and $q$ is impossible if there is one of them we don't sample from.)

Actually, those row operations also encode an algorithm to compute $r$ and $q$, which maybe I could have doped out more directly if I weren't such a matrix nut. The procedure is:

For all $j$:

• $a(j|i) \leftarrow \alpha^{-1}q^*(j|i)$ for all $i$ (rescale top $J$ rows)
• $b(j) \leftarrow 1 - \sum_i a(j|i)$ (subtract top J rows from bottom)
• $b(j) \leftarrow b(j) (\frac{J}{\alpha})^{-1}$ (rescale coeff in bottom row to $(1-\alpha)$
• $b(j) \leftarrow b(j)\alpha/J$ (rescale coeff to $(1-\alpha)$ in bottom row
• $a(j|i) \leftarrow a(j|i) - b(j)$ (wipe out top J of right-hand column)
• $q(j|i) \leftarrow a(j|i)/\alpha$ (convert most of diagonal to 1's)
• $r(j) \leftarrow b(j|i)/(1-\alpha)$ (convert bottom of diagonal to 1)

Let me know if you don't end up with positive probabilities that sum to 1. I suspect that if you specify too small of an $\alpha$, you'll have a bad time with this. You may want to start with a big $\alpha$ and inch it up slowly until something goes wrong. That way you can measure the degree to which your system can be simply described.

• Hi Eric, thanks a bunch. When you say '$\hat{j}$' I assume you mean '$\hat{r}$' ? – ted Jan 22 '18 at 4:51
• Yes, that's jight. – eric_kernfeld Jan 22 '18 at 13:53
• This solves the question (which stated $\alpha$ is fixed) only in the special case where the formulation does not constrain $\hat{p}$ at all (compared to the direct MLE of $p(i\mid j)$). OP does not explain the motivation behind the problem, but what would be the benefit of considering the mixture formulation in the first place, if $\alpha$ is then tuned so that the mixture formulation does not matter? – Juho Kokkala Jan 23 '18 at 7:59
• However, perhaps this can be developed into a full answer by considering the constrained optimization of $q^\ast$ with the constraint implied by the mixture formulation (did not think further whether this is solvable) – Juho Kokkala Jan 23 '18 at 8:00
• There are two motivations for including the zeroth order component of the mixture model. (1) Detrend the chain by identifying states which are commonly visited from all states. So "r" is a background distribution (noise), and "q" is the signal representing transitions which depend on context (the current state). (2) The second motivation is statistical. By including a zeroth order component, we can better predict transitions out of sample, since our training data contains many states which were infrequently visited. – ted Jan 23 '18 at 19:38