For a simple LTI state space system, \begin{equation} \begin{aligned} x_{k+1}&=a x_k + w_k\\ y_k&=x_k +v_k \end{aligned}, \quad k=1,\dots,M \end{equation} where noises $w_k\sim N(0,Q)$ and $v_k\sim N(0,R)$, $a=0.5$, $Q=0.01$ and $R=0.005$. I would like to estimate parameter $a$ and noise covariances $Q$ and $R$ using EM method. I let EM method iterates 50 times and define parameter set \begin{equation} \mathscr{O}_i:=\{\hat{a}_i, \hat{Q}_i, \hat{R}_i\},\quad i=1,\dots,50. \end{equation}
The inital guesses of $a$, $Q$ and $R$ are $\hat{a}_1=0.9$, $\hat{Q}_1=1$ and $\hat{R}_1=1$, respectively. I then use Kalman smoother to determine smoothed estimates $(\hat{w}_{k|M})_{k=1}^M$ and $(\hat{v}_{k|M})_{k=1}^M$. Hence the EM algorithm can be written as \begin{align} &\min_{a,Q,R}\mathcal{Q}(\mathscr{O}|\mathscr{O}_i)=\min_{a,Q,R} (M-1)\log|Q|+M\log|R|+\frac{1}{2}tr\left(R^{-1} \sum_{k=1}^M \hat{v}_{k|M}\hat{v}_{k|M}^\top + P_\infty\right)\\ &+\frac{1}{2}tr\left(Q^{-1} \sum_{k=1}^{M-1} \hat{w}_{k|M}\hat{w}_{k|M}^\top + P_{k+1|M} - aU_\infty P_{k+1|M} - P_{k+1|M}U_\infty^\top a^\top +aP_{k|M}a^\top\right) \end{align} where $P_{\infty+1|\infty}$ is the steady state state error covariance, \begin{equation} P_\infty:=(1-L_\infty)P_{\infty+1|\infty} \quad \text{and} \quad P_{k|M}=P_\infty + U_\infty(P_{k+1|M}-P_{\infty+1|\infty})U_\infty^\top \end{equation} $L_\infty$ and $U_\infty$ are the steady state Kalman filter and smooth gain, respectively.
Therefore, the optimal solution of $a$, $Q$ and $R$ are \begin{align} \hat{a}_{i+1}&= \left(\sum_{k=1}^{M-1} P_{k+1|M} U_\infty\right)\left(\sum_{k=1}^{M-1} P_{k|M}\right)^{-1}\\ \hat{Q}_{i+1}&=\left(\sum_{k=1}^{M-1} \hat{w}_{k|M}\hat{w}_{k|M}^\top + P_{k+1|M} - \hat{a}_{i+1}U_\infty P_{k+1|M} - P_{k+1|M}U_\infty^\top \hat{a}_{i+1}^\top +\hat{a}_{i+1}P_{k|M}\hat{a}_{i+1}^\top\right)/(M-1)\\ \hat{R}_{i+1}&=\left(\sum_{k=1}^M \hat{v}_{k|M}\hat{v}_{k|M}^\top + P_\infty\right)/M \end{align} where $(P_{k|M})_{k=1}^M$, $P_\infty$, $U_\infty$, $(\hat{w}_{k|M})_{k=1}^M$ and $(\hat{v}_{k|M})_{k=1}^M$ are function of known parameter set $\mathscr{O}_i$
I tried to program it in MATLAB, however, as $i$ goes to 50, $a$ always converge to $0$, I have checked it several hours, is there any problem of my understanding of EM method?
P.S. If I assume $a$ is given, then both $\hat{Q}$ and $\hat{R}$ can converge to their true values.
Regards
Stephen.G
I finally found the mistake, when I was calculating $(\hat{w}_{k|M})_{k=1}^M$, I used \begin{equation} \hat{w}_{k|M}=\hat{x}_{k+1|M}-\hat{a}_i\hat{x}_{k|M}, \end{equation} which should be \begin{equation} \hat{w}_{k|M}=\hat{x}_{k+1|M}-a\hat{x}_{k|M}, \end{equation} where $a$ is the decision variable for the EM optimization.