Why is computing $\log p(x)$ difficult, but not the ELBO? This question is in the context where we have some observed data $x$ and some latent variables $z$ which may be used to 'explain' the data. Let's say we have some likelihood model $p(x \vert z)$ and some prior over latent variables $p(z)$. Why is approximating the marginal distribution hard or computationally expensive? Is there something inherently wrong with approximating $p(x)$ using samples from the prior as follows (this seems like it'd be an unbiased estimator as well):
$$ \log p(x) = \log \int dz \, p(x \vert z) p(z) \approx \log \sum_i p(x \vert z_i), \quad z_i \sim p(z) $$
This looks quite similar to how we would compute the evidence lower bound (ELBO), which is used to approximate the marginal log-likelihood by introducing a variational distribution $q(z;x)$:
$$ \log p(x) \geq \textrm{ELBO}(x) = \int dz \, q(z;x) \left[\log p(x \vert z) p(z) - \log q(z;x)\right]$$
The ELBO is usually approximated using multiple samples from $q$:
$$\textrm{ELBO}(x) \approx \log p(x \vert z_i) p(z_i) - \log q(z_i; x), \quad z_i \sim q(z;x)$$
I understand why computing the marginal likelihood exactly may be 'intractable' in common cases - i.e. no analytical form exists if $x$ depends nonlinearly on the $z$s or numerical evaluation of the integrals could take time exponential in the dimension of $z$, but don't understand why the above sample approximation for the ELBO is widely used (in some cases only a single sample is used!), while approximating the marginal likelihood similarly, $p(x) \approx \sum_i p(x \vert z_i)$, is not acceptable/widely used.
Edit: Thinking about it a bit more (and prompted by a comment), I guess $q(z;x)$ is being used to form an importance sampling estimate of the marginal likelihood:
$$ p(x) = \int dz \, p(x \vert z)p(z) = \int dz \, q(z;x) \frac{p(x \vert z) p(z)}{q(z;x)} $$
From which the ELBO is derived by taking logs and using Jensen's inequality, but I don't understand why this should result in a lower variance estimate than the native $p(x) \approx \sum_i p(x \vert z_i)$ estimate.
 A: The simulation method you are looking at can be generalised by using importance sampling.  As a general rule, if you take $Z \sim g$ using some density $g$ then you have:
$$\log p(x) 
= \log \bigg( \int \frac{p(x|z) p(z)}{g(z)} \ g(z) \ dz \bigg)
= \log \bigg( \mathbb{E} \bigg( \frac{p(x|z) p(z)}{g(Z)} \bigg) \bigg).$$
Simulating $Z_1,...,Z_M \sim \text{IID } g$ using some large number of simulations $M$ then gives:
$$\log p(x) 
\approx \log \bigg( \sum_{i=1}^M \frac{p(x|z_i) p(z_i)}{g(z_i)} \bigg) - \log M.$$
This method approximates the true log-mean shown above by taking the log-sample-mean for a simulated sample of values.  As in all applications of importance sampling, the method is most efficient if you choose $g$ so that it is as close as possible to the joint density $p(z,x)$.  The closer this is, the lower the variance of the approximation of the sample mean to the true mean, and so the less simulations are needed to get a good approximation.
This method requires you to compute sums of simulated terms (some of which will be very small), and then take the logarithm of the sum.  Computationally, this means that you will usually want to work in log-space (see e.g., here and here).  It should be possible to do this and get reasonable results without too much loss of accuracy, so I don't see any general need to use crude approximations like ELBO.  The latter is essentially giving you a lower bound by using Jensen's inequality, but I would recommend avoiding this unless you run in to major difficulties computing the log-sum directly from the simulations.
