一个MATLAB 的程序学习

$$
\begin{align}
\mathcal{L}(\bm{\theta}, \bm{\beta}; \bm{y}) &= \int g^{(L)}(\bm{u}^{(L)};\bm{\theta}^{(L)})\prod_{j}{ \bigg[\mathcal{L}^{(L)}_{j}(\bm{\theta}, \bm{\beta} ; \bm{y}| \bm{u}^{(L)})\bigg]^{w^{(L)}_{j}}} d\bm{u}^{(L)}
\end{align}
$$

$$
\begin{align}
f_{\bm{y}, \bm{\upsilon}} \left(\bm{y}, \bm{\upsilon}, \bm{\beta}, \bm{\theta}, \sigma^2 \right)&= \frac{1}{\left(2\pi\sigma^2\right)^{\frac{n_x}{2}}} \exp{\left[\frac{-r^2(\bm{\theta}, \bm{\beta}, \bm{\upsilon}) + ||\upsilon||^2}{2\sigma^2} \right]}  \cdot \frac{1}{\left(2\pi\sigma^2\right)^{\frac{n_z}{2}}} \exp{\left[ \frac{-||\bm{\upsilon}||^2}{2\sigma^2}\right] }  \\
&= \frac{1}{\left(2\pi\sigma^2\right)^{\frac{n_x+n_z}{2}}} \exp{\left[\frac{-r^2(\bm{\theta}, \bm{\beta}, \bm{\upsilon}) }{2\sigma^2}\right] }
\end{align}
$$

$$
\begin{align}
\mathcal{L}\left( \bm{\beta}, \bm{\theta}, \sigma^2; \bm{y} \right)&=\int f_{\bm{y}, \bm{\upsilon}} \left(\bm{y}, \bm{\upsilon}, \bm{\beta}, \bm{\theta}\right) d\bm{\upsilon}\\
&=\int \frac{1}{\left(2\pi\sigma^2\right)^{\frac{n_x+n_z}{2}}} \exp{\frac{-r^2(\bm{\theta}, \bm{\beta}, \bm{\upsilon}) }{2\sigma^2} }  d\bm{\upsilon}\\
&=\frac{1}{\left(2\pi\sigma^2\right)^{\frac{n_x+n_z}{2}}} \exp{\frac{-r^2(\bm{\theta}, \bm{\hat{\beta}}, \bm{\hat{\upsilon}}) - \left|\left|  \bm{R}_{22}  (\bm{\beta} - \hat{\bm{\beta}})  \right| \right|^2 }{2\sigma^2} } \int \exp{\frac{- \left|\left|  \bm{R}_{11} (\bm{\upsilon} - \hat{\bm{\upsilon}}) + \bm{R}_{12}  (\bm{\beta} - \hat{\bm{\beta}})  \right| \right|^2 }{2\sigma^2} }  d\bm{\upsilon} \label{eq:lnlInt}
\end{align}
$$

$$
\begin{align}
\mathcal{L}\left( \bm{\beta}, \bm{\theta}, \sigma^2; \bm{y} \right)&=\frac{1}{\left(2\pi\sigma^2\right)^{\frac{n_x+n_z}{2}}} \exp{\frac{-r^2(\bm{\theta}, \bm{\hat{\beta}}, \bm{\hat{\upsilon}}) - \left|\left|  \bm{R}_{22}  (\bm{\beta} - \hat{\bm{\beta}})  \right| \right|^2 }{2\sigma^2} } \int \exp{\frac{- \left|\left|  \gamma  \right| \right|^2 }{2\sigma^2} } | {\rm det} (\bm{R}_{11}) |^{-1} d\bm{\gamma} \\
&=\frac{1}{|{\rm det}(\bm{R}_{11})| \left(2\pi\sigma^2\right)^{ \frac{n_x}{2}}} \exp{\frac{-r^2(\bm{\theta}, \bm{\hat{\beta}}, \bm{\hat{\upsilon}}) - \left|\left|  \bm{R}_{22}  (\bm{\beta} - \hat{\bm{\beta}})  \right| \right|^2 }{2\sigma^2} }  \left\{ \frac{1}{\left(2\pi\sigma^2\right)^{ \frac{n_z}{2}} } \int \exp{\frac{- \left|\left|  \gamma  \right| \right|^2 }{2\sigma^2} }  d\bm{\gamma} \right\}
\end{align}
$$

$$
\begin{align}
D \left( \bm{\beta}, \bm{\theta}, \sigma^2; \bm{y} \right)&= 2\,{\rm ln} {|{\rm det}(\bm{R}_{11})|} + n_x {\rm ln} \left(2\pi\sigma^2\right) + \frac{r^2(\bm{\theta}, \bm{\hat{\beta}}, \bm{\hat{\upsilon}}) + \left|\left|  \bm{R}_{22}  (\bm{\beta} - \hat{\bm{\beta}})  \right| \right|^2 }{\sigma^2}
\end{align}
$$

$$
\begin{align}
r^2(\bm{\theta}, \bm{\beta}, \bm{\upsilon}) &= \left| \left| \left[ \begin{matrix} \bm{\Omega}^{\frac{1}{2}} \bm{y} \\ \bm{0} \end{matrix} \right] - \left[ \begin{matrix} \bm{\Omega}^{\frac{1}{2}} \bm{Z\Lambda}(\bm{\theta}) & \bm{\Omega}^{\frac{1}{2}}\bm{X} \\ \bm{\Psi}^{\frac{1}{2}} & \bm{0} \end{matrix} \right] \left[ \begin{matrix} \bm{\upsilon} \\ \bm{\beta} \end{matrix} \right] \right|\right|^2 \label{eq:WeMixWr2} \\
&= \left|\left| \bm{\Omega}^{\frac{1}{2}} \left( \bm{y} - \bm{Z\Lambda}(\bm{\theta}) \bm{\upsilon} - \bm{X} \right) \right| \right|^2 + \left| \left| \bm{\Psi}^{\frac{1}{2}} \bm{\upsilon} \right| \right|^2
\end{align}
$$

so, it like this

where $\bm{Z}$ is now a block matrix that incorporates all of the $\bm{Z}$ matrices for the various levels:

The likelihood of $\bm{y}$, conditional on $\bm{\upsilon}$ is then
$$
\begin{align}
f_{\bm{y}|\bm{\upsilon}=\bm{\upsilon}}(\bm{y}, \bm{\upsilon}) &= \prod_{i=1}^{n_x} \left[ {\frac{1}{\left(2\pi\sigma^2\right)^{\frac{1}{2}}} \exp\left[-\frac{\left|\left|\bm{y}_i-\bm{X}_i\bm{\beta}-\bm{Z}_i\bm{\Lambda}(\bm{\theta})\bm{\upsilon} \right|\right|^2}{2\sigma^2}\right) } \right]^{\bm{\Omega}_{ii}} \\
&= \prod_{i=1}^{n_x}  {\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\bm{\Omega}_{ii}}{2}}} \exp \left(- \frac{ \left|\left| \bm{\Omega}_{ii}^{\frac{1}{2}} \left(\bm{y}_i-\bm{X}_i\bm{\beta}-\bm{Z}_i\bm{\Lambda}(\bm{\theta})\bm{\upsilon} \right) \right|\right|^2}{2\sigma^2}\right) } \\
&= {\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii}}{2}}} \exp \left( - \frac{ \left|\left|\bm{\Omega}^{\frac{1}{2}} \left(\bm{y}-\bm{X}\bm{\beta}-\bm{Z}\bm{\Lambda}(\bm{\theta})\bm{\upsilon} \right) \right|\right|^2}{2\sigma^2}\right) }
\end{align}$$
And the unconditional density of $\bm{\upsilon}$ is
$$
\begin{align}
f_{\bm{\upsilon}}(\bm{\upsilon}) &=\prod_{j=1}^{n_z} \left[ \frac{1}{\left(2\pi\sigma^2\right)^{\frac{1}{2}}} \exp \left[- \frac{ \left|\left|  \bm{\upsilon}_j \right|\right|^2}{2\sigma^2}\right]  \right]^{\bm{\Psi}_{jj}} \\
&=\prod_{j=1}^{n_z} \frac{1}{\left(2\pi\sigma^2\right)^{\frac{\bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ \left|\left| \bm{\Psi}_{jj}^{\frac{1}{2}}  \bm{\upsilon}_j \right|\right|^2}{2\sigma^2}\right] \\
&= \frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_j \bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ \left|\left|  \bm{\Psi}^{\frac{1}{2}}  \bm{\upsilon} \right|\right|^2}{2\sigma^2}\right] \label{eq:WeMixWfu}
\end{align}$$
where $\sum \bm{\Psi}_{jj}$ is the sum of the terms in the diagonal matrix $\bm{\Psi}$.

The joint distribution of $\bm{\upsilon}$ and $\bm{y}$ is then the product of eqs. \ref{eq:WeMixWfyu} and \ref{eq:WeMixWfu}:
$$\begin{align}
f_{\bm{y}, \bm{\upsilon}} \left(\bm{y}, \bm{\upsilon}, \bm{\beta}, \bm{\theta}, \sigma^2 \right)&= f_{\bm{y}|\bm{\upsilon}=\bm{\upsilon}}(\bm{y}, \bm{\upsilon})\cdot f_{\bm{\upsilon}}(\bm{\upsilon})\\
&={\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii}}{2}}} \exp \left[- \frac{ \left|\left|\bm{\Omega}^{\frac{1}{2}} \left(\bm{y}-\bm{X\beta}-\bm{Z}\bm{\Lambda}(\bm{\theta})\bm{\upsilon}  \right) \right|\right|^2}{2\sigma^2}\right] } \cdot \frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_j \bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ \left|\left|  \bm{\Psi}^{\frac{1}{2}}  \bm{\upsilon} \right|\right|^2}{2\sigma^2}\right] \\
&=\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii} + \sum_j \bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ \left|\left|\bm{\Omega}^{\frac{1}{2}} \left(\bm{y}-\bm{X\beta}-\bm{Z}\bm{\Lambda}(\bm{\theta})\bm{\upsilon} \right) \right|\right|^2 + \left|\left|  \bm{\Psi}^{\frac{1}{2}}  \bm{\upsilon} \right|\right|^2}{2\sigma^2} \right] \\
&=\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii} + \sum_j \bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ r^2(\bm{\theta}, \bm{\beta}, \bm{\upsilon})}{2\sigma^2} \right]
\end{align}$$
Using the same logic for the results in eq. \ref{eq:r2sub}, $r^2$ can be written as a sum of the value at the optimum ($\hat{\bm{\beta}}$ and $\hat{\bm{\upsilon}}$) and deviations from that:
$$\begin{align}
f_{\bm{y}, \bm{\upsilon}} \left(\bm{y}, \bm{\upsilon}, \bm{\beta}, \bm{\theta}, \sigma^2 \right)&=\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii} + \sum_j \bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ r^2(\bm{\theta}, \hat{\bm{\beta}}, \hat{\bm{\upsilon}}) - \left| \left| \bm{R}_{22}(\bm{\beta} - \hat{\bm{\beta}})\right| \right|^2 - \left| \left| \bm{R}_{11}(\bm{\upsilon} - \hat{\bm{\upsilon}}) + \bm{R}_{12}(\bm{\beta} - \hat{\bm{\beta}}) \right| \right|^2}{2\sigma^2} \right]
\end{align}$$
Now, finding the integral of this over $\bm{\upsilon}$,
$$\begin{align}
\mathcal{L}(\bm{\beta}, \bm{\theta}, \sigma^2; \bm{y})&=\int f_{\bm{y}, \bm{\upsilon}} \left(\bm{y}, \bm{\upsilon}, \bm{\beta}, \bm{\theta}, \sigma^2 \right)  d\bm{\upsilon} \\
&=\int \frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii} + \sum_j \bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ r^2(\bm{\theta}, \hat{\bm{\beta}}, \hat{\bm{\upsilon}}) - \left| \left| \bm{R}_{22}(\bm{\beta} - \hat{\bm{\beta}})\right| \right|^2 - \left| \left| \bm{R}_{11}(\bm{\upsilon} - \hat{\bm{\upsilon}}) + \bm{R}_{12}(\bm{\beta} - \hat{\bm{\beta}}) \right| \right|^2}{2\sigma^2} \right] d\bm{\upsilon} \\
&=\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii} + \sum_j \bm{\Psi}_{jj}}{2}}} \exp \left[- \frac{ r^2(\bm{\theta}, \hat{\bm{\beta}}, \hat{\bm{\upsilon}}) - \left| \left| \bm{R}_{22}(\bm{\beta} - \hat{\bm{\beta}})\right| \right|^2 }{2\sigma^2}\right] \int \exp \left[- \frac{\left| \left|  \bm{R}_{11}(\bm{\upsilon} - \hat{\bm{\upsilon}}) + \bm{R}_{12}(\bm{\beta} - \hat{\bm{\beta}}) \right| \right|^2}{2\sigma^2} \right] d\bm{\upsilon} \label{eq:tmpLA}
\end{align}$$

Notice that while the unweighted integral has $n_z$ dimensions, this weighted integral has $\sum_j \bm{\Psi}_{jj}$ dimensions---the number of (population) individuals values to integrate out.
$$\begin{align}
\mathcal{L}(\bm{\beta}, \bm{\theta}, \sigma^2; \bm{y})
&=\frac{1}{\left(2\pi\sigma^2\right)^{\frac{\sum_i \bm{\Omega}_{ii}}{2}}} \exp \left[- \frac{ r^2(\bm{\theta}, \hat{\bm{\beta}}, \hat{\bm{\upsilon}}) - \left| \left| \bm{R}_{22}(\bm{\beta} - \hat{\bm{\beta}})\right| \right|^2 }{2\sigma^2}\right] \\ & \left\{ \frac{1}{(2\pi\sigma^2)^{\frac{\sum_j \bm{\Psi}_{jj}}{2}}} \int \exp \left[- \frac{\left| \left|  \bm{R}_{11}(\bm{\upsilon} - \hat{\bm{\upsilon}}) + \bm{R}_{12}(\bm{\beta} - \hat{\bm{\beta}}) \right| \right|^2}{2\sigma^2} \right] d\bm{\upsilon} \right\}
\end{align}$$

$$
\left(h\frac{\partial}{\partial x}  + k\frac{\partial}{\partial y} \right) ^m
f(x_0,y_0)  =\sum_{p=0}^{m}C_m^ph^pk^{m-p}
\frac{\partial^mf}{\partial x^p\partial y^{m-p}}\Big|_{x_0-y_0}
$$

$\rightarrow$