Why use the EM Algorithm and not just marginalise the complete likelihood? On the wikipedia article for Expectation-Maximization it states

Given the statistical model which generates a set $\mathbf{X}$ of observed data, a set of unobserved latent data or missing values $\mathbf{Z}$, and a vector of unknown parameters $\boldsymbol\theta$, along with a likelihood function $L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z}) = p(\mathbf{X}, \mathbf{Z}\mid\boldsymbol\theta)$, the maximum likelihood estimate (MLE) of the unknown parameters is determined by maximizing the marginal likelihood of the observed data
$$L(\boldsymbol\theta; \mathbf{X}) = p(\mathbf{X}\mid\boldsymbol\theta) = \int  p(\mathbf{X},\mathbf{Z} \mid \boldsymbol\theta) \, d\mathbf{Z} = \int  p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol\theta) p(\mathbf{Z} \mid \boldsymbol\theta) \, d\mathbf{Z} $$
However, this quantity is often intractable since $\mathbf{Z}$ is unobserved and the distribution of $\mathbf{Z}$ is unknown before attaining $\boldsymbol\theta$.

I don't understand this last sentence. Surely $\mathbf{Z}$ being unobserved is why we integrate over $\mathbf{Z}$, and we set values of $\boldsymbol{\theta}$ during the optimisation procedure which then allows us to compute $p(\mathbf{Z} \mid \boldsymbol{\theta})$?
I initially thought it was intractable because there are too many settings of $\mathbf{Z}$ to consider, but in the expectation step of the EM algorithm we calculate
\begin{align}
  Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)}) &= \mathbb{E}_{\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)}}\left[\log L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z})\right]\\
  &= \int_Z p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)})\log L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z})d\mathbf{Z}
 \end{align}
which doesn't look any easier than the first integral. If we can approximate this expectation using some method like monte-carlo, why not just approximate the first integral instead?

So I've discussed this with a colleague.
Consider the marginalisation
$$p(\mathbf{X} \mid \boldsymbol{\theta}) = \int p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\theta})p(\mathbf{Z} \mid \boldsymbol{\theta}) d\mathbf{Z}.$$
This can be rewritten as the expectation
$$\mathbb{E}_{\mathbf{Z} \mid \boldsymbol{\theta}}\left[ p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\theta}) \right].$$
If this can be calculated exactly, we're fine. However, if the expectation is intractable, we need to compute a numerical approximation. Hence to maximise w.r.t. $\boldsymbol{\theta}$ we will need to use some iterative procedure like gradient-ascent.
Suppose we have the current estimate $\boldsymbol{\hat{\theta}^{(t)}}$ for $\boldsymbol{\theta}$. If we now try and approximate the expectation by sampling from $\mathbf{Z}$, we get
$$\frac{1}{n}\sum_{i=1}^n p(\mathbf{X} \mid \mathbf{Z}_i, \boldsymbol{\theta})$$
but the $\mathbf{Z}_i$ were sampled from the distribution $\mathbf{Z} \mid \boldsymbol{\hat{\theta}^{(t)}}$, so this actually approximates the expectation
$$\mathbb{E}_{\mathbf{Z} \mid \boldsymbol{\hat{\theta}^{(t)}}}\left[ p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\theta}) \right] \approx \frac{1}{n}\sum_{i=1}^n p(\mathbf{X} \mid \mathbf{Z}_i, \boldsymbol{\theta})$$
which is not the same as the original expectation we wanted. Hence if we maximised this new expectation w.r.t. $\boldsymbol{\theta}$ we won't be doing Maximum Likelihood Estimation.
The EM algorithm therefore takes a different approach.
