Unbiased Estimator for $\log\left[\int p(x\mid z)p(z) \, dz\right]$ The naive Monte Carlo estimator is an unbiased estimator for $\int p(x\mid z)p(z) \, dz$, is there a convenient unbiased estimator for $\log \left[\int p(x\mid z)p(z)\,dz \right]$
 A: Path sampling is a way to evaluate the log integral by an unbiased estimator. Let us introduce a temperature index $0\le t\le 1$ and a sequence of conditional functions $p_t(x|z)$ such that
$$p_0(x|z)=1\qquad\qquad\text{and}\qquad\qquad p_1(x|z)=p(x|z)$$
Then, if $$\mathfrak{Z}_t(x)=\int p_t(x|z)\,q(z)\text{d}z$$
\begin{align*}\frac{\text{d}}{\text{d}t} \log(\mathfrak{Z}_t(x)) 
&= \frac{1}{\mathfrak{Z}_t(x)}\frac{\text{d}}{\text{d}t}\mathfrak{Z}_t(x)\\
&=\frac{1}{\mathfrak{Z}_t(x)}\int\frac{\text{d}}{\text{d}t}p_t(x|z)\,q(z)\text{d}z\\
&=\frac{1}{\mathfrak{Z}_t(x)}\int\frac{\text{d}}{\text{d}t}p_t(x|z)\,q(z)\text{d}z\\
&=\frac{1}{\mathfrak{Z}_t(x)}\int\frac{\text{d}}{\text{d}t}\log(p_t(x|z))\,p_t(x|z)q(z)\text{d}z\\
&=\mathbb{E}\left[\left.\frac{\text{d}}{\text{d}t}\log(p_t(x|Z))\right| X=x\right]\end{align*}
Therefore, since
$$\log(\mathfrak{Z}_1(x)/\mathfrak{Z}_0(x)) = \int_0^1 \frac{\text{d}}{\text{d}t} \log(\mathfrak{Z}_t(x))\,\text{d}t$$
and $\mathfrak{Z}_0(x)=1$, a remarkable identity is that
$$\log(\mathfrak{Z}_1(x))=\log\left\{\int p(x|z) q(z)\text{d}z\right\}=\int_0^1  \mathbb{E}\left[\left.\frac{\text{d}}{\text{d}t}\log(p_t(x|Z))\right| X=x\right]\,\text{d}t$$
Hence, assuming simulations $z_i=\zeta(x,t,\epsilon_i)$ ($i=1,\ldots,I$) from $$q(z|x)=\frac{q(z) p(x|z)}{\mathfrak{Z}_1(x))}$$ are available (e.g. by MCMC), an unbiased estimator of $\log(\mathfrak{Z}_1(x))$ is
$$\int_0^1 \frac{1}{I}\sum_{i=1}^I\frac{\text{d}}{\text{d}t}\log(p_t(x|\zeta(x,t,\epsilon_i)))\,\text{d}t$$
A: It can be shown that by using sufficiently large sample sizes for the MC approximation, a lower bound to the marginal log-likelihood is tightened. While this estimator is hence biased, in can serve the same practical purposes. It was shown in [1] that
$$
\log p(x) \ge \mathcal{L}_{K+1} \ge \mathcal{L}_{K} 
$$
for $\mathcal{L}_K = \mathbb E_{z_1, \ldots, z_K \sim q^K} \left[\log \frac{1}{K} \sum_{k=1}^K \frac{p(x, z_k)}{q(z_k)}\right]$, where $q^K \equiv \underbrace{q \otimes q \otimes \cdots \otimes q}_{K times}$. 
This also holds for $q(z) = p(z)$.
[1] Burda, Yuri, Roger Grosse, and Ruslan Salakhutdinov. "Importance weighted autoencoders." arXiv preprint arXiv:1509.00519 (2015).
A: The Monte Carlo estimator using $n$ samples of
$$l = \int p(x \mid z) g(z) dz = E(p(x \mid Z))$$
with $Z\sim g$ is
$$\hat l_n = \frac 1n \sum_{i = 1}p(x\mid Z_i) \qquad Z_i\sim g$$
Following Durbin and Koopman
(1997), the error of the log
marginal likelihood is given by
$$\begin{align*}
\log\hat l_n - \log l
  &= \log\left(1 + \frac{\hat l_n - l}{l}\right) \\
  &= \frac{\hat l_n - l}{l}
    - \frac 12\left(\frac{\hat l_n - l}{l}\right)^2 +
    O\left(\left(\frac{\hat l_n - l}{l}\right)^3\right)
\end{align*}$$
implying that the bias is
$$
E(\log\hat l_n) =  \log l
  - \frac 12\frac{E\left((\hat l_n - l)^2\right)}{l^2} +
    O\left(\frac{E\left((\hat l_n - l)^3\right)}{l^3}\right)
$$
Thus, the bias is approximately proportional to the variance
of the estimated log-likelihood using the delta method and
that $E(\hat l_n) =l$.
The bias corrected version from Durbin and Koopman (1997) is
therefor
$$
\log\hat l_n + \frac 12 \frac{\widehat{\text{Var}}(\hat l_n)}{\hat l_n^2}
$$
where $\widehat{\text{Var}}$ is the estimated variance. That is,
$$\widehat{\text{Var}}(\hat l_n) = \frac 1{n(n -1)} \sum_{i = 1}^n\left(p(x\mid Z_i) - \hat l_n\right)^2.$$
The adjusted estimator is convenient as we already have all the quantities we need from the standard Monte Carlo estimator.
The above can also be
extended to importance sampling. Pick an importance distribution with density $h$ with the same support as $g$. Then replace $p(x\mid Z_i)$ in $\hat l_n$ (and the variance estimate) with
$$
\frac{p(x\mid Z_i)g(Z_i)}{h(Z_i)}
$$
and sample $Z_i\sim h$.
Example
Take $g(z) = \phi(z)$ where $\phi$ is the standard normal
distribution's density function and
$$
p(x\mid z) = \Phi(0.5 + z)\Phi(-0.25 - z)
$$
where $\Phi$ is the standard normal distribution's cumulative density function. In this case $l$ is the bivariate normal CDF integrated up to 0.5 and -0.25 with covariance matrix
$$\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}$$
which can be computed to floating point precision which makes this example nice.
A small simulation study investigating the bias of the adjusted and unadjusted Monte Carlo estimator as a function of the sample size is run below.
# the true value of l
truth <- mvtnorm::pmvnorm(
  upper = c(.5, -.25), sigma = matrix(c(2, -1, -1, 2), 2))

n_reps <- 100000L # the number of replications
# the sample sizes
n_samples <- seq(log(10), log(1000), length.out = 40) |> exp() |> floor()

# compute the estimates
set.seed(1)
ests <- sapply(n_samples, \(n_samp){
  smps <- matrix(rnorm(n_samp * n_reps), n_samp)
  weight <- pnorm(.5 + smps) * pnorm(-0.25 - smps)
  rbind(Estimate = colMeans(weight), Var = apply(weight, 2L, var) / n_samp)
}, simplify = "array")

# compute the bias estimates
err_unadjusted <- log(ests["Estimate", , ]) - log(truth)
bias_ests_unadjusted <- colMeans(err_unadjusted)
err_adjusted <- 
  log(ests["Estimate", , ]) + ests["Var", , ] / (2 * ests["Estimate", , ]^2) - 
  log(truth)
bias_ests_adjusted <- colMeans(err_adjusted)

# plot the bias estimates versus the sample size
ylim <- range(bias_ests_unadjusted, bias_ests_adjusted)
par(mar = c(5, 5, 1, 1))
plot(n_samples, bias_ests_unadjusted, pch = 16, bty = "l", ylim = ylim,
     ylab = "Bias estimate (unadjusted)", xlab = "Sample size", log = "x")
abline(h = 0, lty = 2)
plot(n_samples, bias_ests_adjusted, pch = 16, bty = "l", ylim = ylim,
     ylab = "Bias estimate (adjusted)", xlab = "Sample size", log = "x")
abline(h = 0, lty = 2)



The plots show that the adjusted estimator seems to be unbiased but the standard Monte Carlo estimator is downward biased consistent with the answer by bayerj. However, the bias of the latter gets smaller as the number of samples is increased also consistent with the answer by bayerj.
A side note, the plot above shows that the variance of the log-likelihood estimator is about $2 \cdot 0.012$ in the worst case when we use $n = 10$ samples. Thus,
the standard error of the estimator is $\sqrt{2 \cdot 0.012} = 0.125$. Hence, a bias of $0.012$ may not be the main concern.
