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.
