expectation of log of expectation by Monte Carlo When considering the approximation by Monte Carlo of an expectation of the form$$\mathfrak{I}=\mathbb{E}^X[\log\{\mathbb{E}^{Y|X}[h(X,Y)|X]\}]$$using a resolution of the form$$\hat{\mathfrak{I}}=\frac{1}{N}\sum_{n=1}^N\log\left\{\frac{1}{M}\sum_ {m=1}^M h(x_n,y_{nm})\right\}$$(when the $x_n$'s are generated from the correct marginal on $X$ and the $y_{nm}$'s from the proper conditional of $Y$ given $X=x_n$) produces a convergent but biased estimator. I wonder if there is a genuine way to produce an unbiased estimator in this case, other than using the trick of McLeish (see equation (1) on page 2).
 A: The methods presented in this work (https://people.maths.ox.ac.uk/gilesm/files/SLOAN80-056.pdf) concern Multi-Level Monte Carlo (MLMC) methods for expectations of this form. MLMC is typically not designed to provide unbiased estimators per se, but can usually be modified to do so using the trick of McLeish.
Broadly, if you are interested in 
$$\mathfrak{I}=\mathbb{E}^X \left[ g \left\{\mathbb{E}^{Y|X}[h(X,Y)|X] \right\} \right]$$
then the standard trick is to write
$$\hat{\mathfrak{I}}_{\ell_1, \ell_2} = \frac{1}{N_{\ell_1}} \sum_{i = 1}^{N_{\ell_1}} g \left( \frac{1}{N_{\ell_2}} \sum_{j = 1}^{N_{\ell_2}} h \left(X^i, Y^{i, j} \right) \right)$$
where $N_{\ell_1}, N_{\ell_2}$ are two sequences of positive integers increasing to infinity. One can then form a sort of double-telescoping sum representation of 
$$\lim_{\min \left(N_{\ell_1}, N_{\ell_2} \right) \to \infty} \hat{\mathfrak{I}}_{\ell_1, \ell_2}  = \mathfrak{I}$$
and the debiasing trick of McLeish applies. The work of Crisan et al. (https://arxiv.org/abs/1702.03057) is also relevant here.
