Solution (source code)

= Solution

Put $w_k=C_k^{-1}$, $u_i=H_i^{-1}$, $A=\sum_kw_k$, $B=\sum_i u_i$, and form the <weighted means>
$$
\bar q_w=\frac{\sum_kw_kq_k}{A},
\qquad
\bar r_u=\frac{\sum_i u_ir_i}{B}.
$$
The two equations obtained from the <score function> are
$$
A(\bar q_w-M_0)+B(\bar r_u-M_0+\theta)=0,
\qquad
-B(\bar r_u-M_0+\theta)=0.
$$
Consequently the <maximum-likelihood estimators> are
$$
\boxed{\widehat M_0=\bar q_w,
\qquad \widehat\theta=\bar q_w-\bar r_u.}
$$
The <Hessian matrix> of the log likelihood is
$$
\begin{pmatrix}-(A+B)&B\\B&-B\end{pmatrix}.
$$
Its first leading principal minor is negative and its <determinant> is $AB>0$, so it is a <negative-definite matrix>. Thus the stationary point is the unique global maximum.