105 Expectation Maximization (EM) Principles
NONMEM 7 Technical Guide
Maximization-expectation methods separate the process of expectation (integration) and maximization. To find improved estimates for \(\boldsymbol{\mu}\) modeled \(\boldsymbol{\theta}_\mu\), it is convenient to first minimize the negative logarithm of \(p(\mathbf{y} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}, \boldsymbol{\Omega})\) with respect to \(\boldsymbol{\mu}\), which is equivalent to maximizing \(p(\mathbf{y} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}, \boldsymbol{\Omega})\). We can do this as follows:
\[ \frac{\partial L_i}{\partial \boldsymbol{\mu}_i} = \tag{1}\]
\[ \frac{ \partial -\log\left(\int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi}\right) }{ \partial \boldsymbol{\mu}_i } = \tag{2}\]
\[ \frac{ -\partial\left(\int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi}\right) / \partial \boldsymbol{\mu}_i }{ \int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} } = \frac{ \int_{-\infty}^{\infty} \bigl[-\partial p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) / \partial \boldsymbol{\mu}_i\bigr] d\boldsymbol{\phi} }{ \int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} } = \tag{3}\]
\[ \frac{ \int_{-\infty}^{\infty} \bigl[-\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) / \partial \boldsymbol{\mu}_i\bigr] p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} }{ \int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} } = \tag{4}\]
\[ \int_{-\infty}^{\infty} \left[\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\mu}_i }\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \tag{5}\]
\[ E_{\boldsymbol{\phi}, i}\left(\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\mu}_i }\right) = \mathbf{g}_{\boldsymbol{\mu}_i} \tag{6}\]
where \(\mathbf{g}_{\boldsymbol{\mu}_i}\) is the gradient with respect to \(\boldsymbol{\mu}_i\), and \(E_{\boldsymbol{\phi}, i}()\) represents the expected value when integrating over all \(\boldsymbol{\phi}\), and
\[ z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) = \frac{ p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) }{ \int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} } \tag{7}\]
is called the conditional density of \(\boldsymbol{\phi}\) for individual \(i\). The conditional density integrated over all possible \(\boldsymbol{\phi}\) evaluates to 1:
\[ \begin{aligned} \int_{-\infty}^{\infty} z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} &= \int_{-\infty}^{\infty} \frac{ p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega}) }{ \int_{-\infty}^{\infty} p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} } d\boldsymbol{\phi} \\ &= 1 \end{aligned} \tag{8}\]
The relationship
\[ \frac{ \partial -\log\bigl(p(\mathbf{y}_i \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\mu}_i } = E_{\boldsymbol{\phi}, i}\left(\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\mu}_i }\right) \tag{9}\]
holds for any joint density \(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\). Now, to evaluate specifically for a parameter density \(h\) that is multivariate normal:
\[ \int_{-\infty}^{\infty}\left[\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\mu}_i }\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \tag{10}\]
\[ \int_{-\infty}^{\infty}\left[\frac{ -\partial \log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\mu}_i }\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \tag{11}\]
\[ \int_{-\infty}^{\infty} -\boldsymbol{\Omega}^{-1}(\boldsymbol{\phi} - \boldsymbol{\mu}_i) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \tag{12}\]
\[ -\boldsymbol{\Omega}^{-1} \int_{-\infty}^{\infty} (\boldsymbol{\phi} - \boldsymbol{\mu}_i) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \mathbf{g}_{\boldsymbol{\mu}_i} \tag{13}\]
We can perform the above algebraic manipulation because \(\boldsymbol{\mu}\) (and therefore \(\boldsymbol{\theta}_\mu\)) appears only in the parameter density \(h\), but does not appear in the data density \(l\). We define
\[ \bar{\boldsymbol{\phi}}_i = \int_{-\infty}^{\infty} \boldsymbol{\phi} z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \tag{14}\]
as the conditional mean \(\boldsymbol{\phi}\) vector for individual \(i\), so that
\[ \begin{aligned} \frac{\partial L_i}{\partial \boldsymbol{\mu}_i} &= -\boldsymbol{\Omega}^{-1} \int_{-\infty}^{\infty} (\boldsymbol{\phi} - \boldsymbol{\mu}_i) z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \\ & = -\boldsymbol{\Omega}^{-1}\left[\int_{-\infty}^{\infty} \boldsymbol{\phi} z(\boldsymbol{\phi} \mid \cdots) d\boldsymbol{\phi} - \boldsymbol{\mu}_i \int_{-\infty}^{\infty} z(\boldsymbol{\phi} \mid \cdots) d\boldsymbol{\phi}\right] \\ &= -\boldsymbol{\Omega}^{-1}(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) = \mathbf{g}_{\boldsymbol{\mu}_i} \end{aligned} \tag{15}\]
There are several ways of determining \(\bar{\boldsymbol{\phi}}_i\) which are described later, and are called the expectation (=integrating or averaging) step.
1 Maximization
To determine the \(\boldsymbol{\mu}\) modeled theta that reduces the objective function, we must solve:
\[ \begin{aligned} \frac{\partial L}{\partial \boldsymbol{\theta}_\mu} &= \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \frac{\partial L_i}{\partial \boldsymbol{\mu}_i} \\ & = -\sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1} (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i) \\ &= \mathbf{g}_{\boldsymbol{\theta}_\mu} \end{aligned} \tag{16}\]
So that
\[ \mathbf{g}_{\boldsymbol{\theta}_\mu} = \mathbf{0} \tag{17}\]
To evaluate Equation 16 fully, an optimization algorithm is necessary which varies \(\boldsymbol{\theta}_\mu\), and evaluating \(L\) at each \(\boldsymbol{\theta}_\mu\). Keep in mind that in addition to \(\boldsymbol{\mu}_i\) varying with \(\boldsymbol{\theta}_\mu\), \(\bar{\boldsymbol{\phi}}_i\) also varies
with \(\boldsymbol{\theta}_\mu\) through the conditional density \(z\), so this minimization process can be computationally expensive. Alternatively, we can perform a limited maximization step in which \(\bar{\boldsymbol{\phi}}_i\) is kept constant, while only \(\boldsymbol{\mu}_i\) is varied with changes in \(\boldsymbol{\theta}_\mu\). This separation of the expectation step from the maximization step is characteristic of the EM algorithm.
Evaluating Equation 16 by this limited optimization is equivalent to minimizing the following surrogate objective function (keeping \(\bar{\boldsymbol{\phi}}_i\) constant):
\[ \begin{aligned} & E_{\boldsymbol{\phi}, i}\bigl(-\log(h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega}))\bigr) = \\ &\quad \frac{1}{2} E_{\boldsymbol{\phi}, i}\left(\sum_{i=1}^{m} (\boldsymbol{\phi} - \boldsymbol{\mu}_i(\boldsymbol{\theta}_\mu))^{\prime} \boldsymbol{\Omega}^{-1} (\boldsymbol{\phi} - \boldsymbol{\mu}_i(\boldsymbol{\theta}_\mu))\right) + \frac{1}{2} E_{\boldsymbol{\phi}, i}\bigl(m \log(\det(\boldsymbol{\Omega}))\bigr) = \\ &\quad \frac{1}{2} \sum_{i=1}^{m}\left[\overline{\boldsymbol{\phi}^{\prime}\boldsymbol{\Omega}^{-1}\boldsymbol{\phi}}_i - 2\boldsymbol{\mu}_i^{\prime} \boldsymbol{\Omega}^{-1} \bar{\boldsymbol{\phi}}_i + \boldsymbol{\mu}_i^{\prime} \boldsymbol{\Omega}^{-1} \boldsymbol{\mu}_i\right] + \frac{1}{2} m \log(\det(\boldsymbol{\Omega})) = \\ &\quad \frac{1}{2} \sum_{i=1}^{m}\left[\overline{\boldsymbol{\phi}^{\prime}\boldsymbol{\Omega}^{-1}\boldsymbol{\phi}}_i - 2\boldsymbol{\mu}_i^{\prime} \boldsymbol{\Omega}^{-1} \bar{\boldsymbol{\phi}}_i + \boldsymbol{\mu}_i^{\prime} \boldsymbol{\Omega}^{-1} \boldsymbol{\mu}_i\right] + \frac{1}{2} m \log(\det(\boldsymbol{\Omega})) + \frac{1}{2} \sum_{i=1}^{m}\left[\bar{\boldsymbol{\phi}}_i^{\prime} \boldsymbol{\Omega}^{-1} \bar{\boldsymbol{\phi}}_i - \bar{\boldsymbol{\phi}}_i^{\prime} \boldsymbol{\Omega}^{-1} \bar{\boldsymbol{\phi}}_i\right] \propto \\ &\quad \frac{1}{2} \sum_{i=1}^{m} (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i(\hat{\boldsymbol{\theta}}_\mu))^{\prime} \boldsymbol{\Omega}^{-1} (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i(\hat{\boldsymbol{\theta}}_\mu)) \equiv L_c \end{aligned} \tag{18}\]
The \(L_c\) is called the (negative) complete data log likelihood, and it can be shown (see (Schumitzky 1995)) that any \(\hat{\boldsymbol{\theta}}_\mu\) that reduces \(L_c\), will reduce \(L\) by an at least equivalent amount, or:
\[ L(\boldsymbol{\theta}_\mu) - L(\hat{\boldsymbol{\theta}}_\mu) \geq L_c(\boldsymbol{\theta}_\mu) - L_c(\hat{\boldsymbol{\theta}}_\mu) \tag{19}\]
where \(\hat{\boldsymbol{\theta}}_\mu\) is an improved value over the present value \(\boldsymbol{\theta}_\mu\). That is, any improvement value \(\hat{\boldsymbol{\theta}}_\mu\) that reduces \(L_c\) (where \(\bar{\boldsymbol{\phi}}_i\) was kept constant), will also reduce \(L\) (in which \(\bar{\boldsymbol{\phi}}_i\) varies with \(\boldsymbol{\theta}_\mu\)), by at least the same amount as it reduced \(L_c\).
The easiest way to minimize \(L_c\) is to perform a least squares analysis, by producing the following positive definite Hessian matrix:
\[ \begin{aligned} \mathbf{H}_{\boldsymbol{\theta}_\mu} &= E\left(\frac{\partial^2 L_c}{\partial \boldsymbol{\theta}_\mu^2}\right) \\ &= \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } E\left(\frac{\partial^2 L_{ci}}{\partial \boldsymbol{\mu}_i^2}\right) \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \\ &= \sum_{i=1}^{m} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \boldsymbol{\Omega}^{-1} \frac{ \partial \boldsymbol{\mu}_i }{ \partial \boldsymbol{\theta}_\mu } \end{aligned} \tag{20}\]
And performing the following update with a variable step size \(\alpha \leq 1\):
\[ \hat{\boldsymbol{\theta}}_\mu = \boldsymbol{\theta}_\mu + \alpha \mathbf{H}_{\boldsymbol{\theta}_\mu}^{-1} \mathbf{g}_{\boldsymbol{\theta}_\mu} \tag{21}\]
This is the maximization step of the EM algorithm. If all of the mu’s have linear relationships with respect to theta, then the step size that minimizes \(L_c\) with respect to the mu’s is \(\alpha = 1\). However, if the mu’s are not linearly related to thetas, then \(\alpha\) must be adjusted to minimize \(L_c\) with respect to mu. This can be done by selecting a value \(\alpha\), evaluate \(\hat{L}_c\) using the proposed \(\hat{\boldsymbol{\theta}}_\mu\), and if \(\hat{L}_c\) is not smaller than the present \(L\) evaluated at the present \(\boldsymbol{\theta}_\mu\), try another value of \(\alpha\), etc. In NONMEM, \(\alpha = 1\) is first selected, tested, and if necessary, \(\alpha\) is reduced by geometrical decrements of square root of 2 until an \(\hat{L}_c\) is found that is less than \(L_c\). More elaborate search algorithms (such as conjugate gradient or variable metric methods) for thetas not linearly modeled with respect to mu could be used for the expectation-maximization methods, but no real time savings occurs in doing so for population analysis problems.
In the next iteration, the updated \(\boldsymbol{\theta}_\mu\) are used to evaluate a new set of conditional means \(\bar{\boldsymbol{\phi}}_i\) in the expectation step, followed by a limited maximization step to update \(\boldsymbol{\theta}_\mu\) again. By repeatedly performing the expectation step 14, and evaluating the maximization step as expressed in equations 15 through 21, the gradient \(\mathbf{g}_{\boldsymbol{\theta}_\mu}\) becomes smaller, and estimates \(\hat{\boldsymbol{\theta}}_\mu\) that maximize the marginal density (satisfy equation 17) are eventually obtained (Schumitzky 1995).
Again, because \(\boldsymbol{\mu}\) appeared only in the parameter density \(h\) as the mean to this multivariate normal density, and does not appear in the data density \(l\), and the parameters \(\boldsymbol{\theta}_\mu\) to be estimated appear in the objective function only through \(\boldsymbol{\mu}\), this allowed us to obtain a gradient evaluation with a simple construction as given in Equation 16. For those \(\boldsymbol{\theta}\) that are not expressed in the model through \(\boldsymbol{\mu}\) the \(\boldsymbol{\theta}\) may appear anywhere in the joint density. No shortcut evaluation can be made by maximizing just the parameter density portion. Thus, to optimize the population objective function in these \(\boldsymbol{\theta}\) as well, we need to differentiate the entire joint density. Through a similar process as we showed in differentiating with respect to \(\boldsymbol{\mu}\),
\[ \frac{ \partial L_i }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } = \tag{22}\]
\[ \int_{-\infty}^{\infty}\left[\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} }\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \tag{23}\]
\[ E_{\boldsymbol{\phi}, i}\left(\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} }\right) = \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} \tag{24}\]
\[ \begin{aligned} \frac{ \partial L }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } &= \sum_{i=1}^{m} \frac{ \partial L_i }{ \partial \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} } \\ &= \sum_{i=1}^{m} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} \\ &= \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \end{aligned} \tag{25}\]
A Hessian matrix may be constructed as follows:
\[ \mathbf{H}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} = \sum_{i=1}^{m} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}, i}^{\prime} \tag{26}\]
\[ \hat{\boldsymbol{\theta}}_{\enclose{horizontalstrike}{\mu}} = \boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}} + \mathbf{H}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}}^{-1} \mathbf{g}_{\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}} \tag{27}\]
To minimize the objective function with respect to the inter-subject variance parameters, we recognize that \(\boldsymbol{\Omega}\) is symmetrical, and we must vary only the lower triangular portion of the matrix. Defining \(\mathbf{A}\) as the lower triangular matrix of \(\boldsymbol{\Omega}\), and minimizing with respect to \(\mathbf{A}\), we have
\[ \frac{ \partial -\log\bigl(p(\mathbf{y} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \mathbf{A} } = \sum_{i=1}^{m} E_{\boldsymbol{\phi}, i}\left(\frac{ -\partial \log\bigl(p(\mathbf{y}_i, \boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \mathbf{A} }\right) = \tag{28}\]
\[ \sum_{i=1}^{m} \int_{-\infty}^{\infty}\left[\frac{ -\partial \log\bigl(l(\mathbf{y}_i \mid \boldsymbol{\phi}) h(\boldsymbol{\phi} \mid \boldsymbol{\mu}_i, \boldsymbol{\Omega})\bigr) }{ \partial \mathbf{A} }\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \tag{29}\]
\[ \sum_{i=1}^{m} \int_{-\infty}^{\infty} \operatorname{Lower}\left[\mathbf{R}_i - \frac{1}{2} \operatorname{diag}(\mathbf{R}_i)\right] z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} = \mathbf{g}_\mathbf{A} \tag{30}\]
where
\[ \mathbf{R}_i = \boldsymbol{\Omega}^{-1}\bigl(\boldsymbol{\Omega} - (\boldsymbol{\phi} - \boldsymbol{\mu}_i)(\boldsymbol{\phi} - \boldsymbol{\mu}_i)^{\prime}\bigr) \boldsymbol{\Omega}^{-1} \tag{31}\]
and \(\mathbf{g}_\mathbf{A}\) is the gradient with respect to \(\mathbf{A}\). The derivation from equation 29 to 30 requires evaluating partial derivatives of matrix components, the tools for which are derived in appendix A.
We define
\[ \begin{aligned} \bar{\boldsymbol{\Omega}}_i &= \int_{-\infty}^{\infty} (\boldsymbol{\phi} - \boldsymbol{\mu}_i)(\boldsymbol{\phi} - \boldsymbol{\mu}_i)^{\prime} z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \\ &= E_{\boldsymbol{\phi}, i}\bigl[(\boldsymbol{\phi} - \boldsymbol{\mu}_i)(\boldsymbol{\phi} - \boldsymbol{\mu}_i)^{\prime}\bigr] \end{aligned} \tag{32}\]
as the contribution to the evaluated population variance from each individual \(i\). Then,
\[ E_{\boldsymbol{\phi}, i}(\mathbf{R}_i) = \boldsymbol{\Omega}^{-1}\bigl(\boldsymbol{\Omega} - \bar{\boldsymbol{\Omega}}_i\bigr) \boldsymbol{\Omega}^{-1} \tag{33}\]
and
\[ \begin{aligned} \operatorname{Lower}\left[\sum_{i=1}^{m} E_{\boldsymbol{\phi}, i}(\mathbf{R}_i) - \frac{1}{2} \operatorname{diag}\left(\sum_{i=1}^{m} E_{\boldsymbol{\phi}, i}(\mathbf{R}_i)\right)\right] &= \mathbf{g}_\mathbf{A} \\ &= \mathbf{0} \end{aligned} \tag{34}\]
is equivalent to solving for
\[ \begin{aligned} \sum_{i=1}^{m} E_{\boldsymbol{\phi}, i}(\mathbf{R}_i) &= \boldsymbol{\Omega}^{-1}\left(m\boldsymbol{\Omega} - \sum_{i=1}^{m} \bar{\boldsymbol{\Omega}}_i\right) \boldsymbol{\Omega}^{-1} \\ &= \mathbf{0} \end{aligned} \tag{35}\]
which suggests the following update for \(\boldsymbol{\Omega}\):
\[ \hat{\boldsymbol{\Omega}} = \frac{1}{m} \sum_{i=1}^{m} \bar{\boldsymbol{\Omega}}_i \tag{36}\]
Note for any given \(\boldsymbol{\Omega}\)
\[ \mathbf{g}_\mathbf{A} = \operatorname{Lower}\left[m\boldsymbol{\Omega}^{-1}(\boldsymbol{\Omega} - \hat{\boldsymbol{\Omega}})\boldsymbol{\Omega}^{-1} - \frac{m}{2} \operatorname{diag}\left(\boldsymbol{\Omega}^{-1}(\boldsymbol{\Omega} - \hat{\boldsymbol{\Omega}})\boldsymbol{\Omega}^{-1}\right)\right] \tag{37}\]
Thus, with repeatedly evaluating the expectation step 32, and utilizing the result to evaluate the next estimate of the intersubject variance (maximization step 36), when the “output” \(\hat{\boldsymbol{\Omega}}\) equals the “input” \(\boldsymbol{\Omega}\), then the gradient \(\mathbf{g}_\mathbf{A}\) is equal to \(\mathbf{0}\).
Note that equation 32 may be rearranged as follows, which will be useful later:
\[ \bar{\boldsymbol{\Omega}}_i = (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} + \int_{-\infty}^{\infty} (\boldsymbol{\phi} - \bar{\boldsymbol{\phi}}_i)(\boldsymbol{\phi} - \bar{\boldsymbol{\phi}}_i)^{\prime} z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \tag{38}\]
Defining
\[ \bar{\mathbf{B}}_i = \int_{-\infty}^{\infty} (\boldsymbol{\phi} - \bar{\boldsymbol{\phi}}_i)(\boldsymbol{\phi} - \bar{\boldsymbol{\phi}}_i)^{\prime} z(\boldsymbol{\phi} \mid \mathbf{y}_i, \boldsymbol{\mu}_i, \boldsymbol{\Omega}) d\boldsymbol{\phi} \tag{39}\]
as the conditional variance of \(\boldsymbol{\phi}\) for individual \(i\), then
\[ \bar{\boldsymbol{\Omega}}_i = (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} + \bar{\mathbf{B}}_i \tag{40}\]
so that
\[ \hat{\boldsymbol{\Omega}} = \frac{1}{m} \sum_{i=1}^{m} (\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)(\bar{\boldsymbol{\phi}}_i - \boldsymbol{\mu}_i)^{\prime} + \frac{1}{m} \sum_{i=1}^{m} \bar{\mathbf{B}}_i \tag{41}\]
Thus, the update variance inter-subject variance is evaluated as the sum of the sample variance of the conditional means and the average conditional variance. To summarize, the EM algorithm consists of an expectation step evaluating conditional means \(\bar{\boldsymbol{\phi}}_i\) and conditional variances \(\bar{\mathbf{B}}_i\), keeping \(\boldsymbol{\theta}_{\enclose{horizontalstrike}{\mu}}\) and \(\boldsymbol{\Omega}\) constant, followed by a limited maximization step to obtain updated \(\boldsymbol{\theta}_\mu\) and \(\boldsymbol{\Omega}\), keeping \(\bar{\boldsymbol{\phi}}_i\) and \(\bar{\mathbf{B}}_i\) constant.