and $\pi(z)$ is approximated by the empirical distribution of a large sample of $Z$. Therefore, by evaluating the quality of the $Z$ sample, we can study the accuracy of the learned $h_{\theta}$ function.

In our experiment, we consider three true latent distributions and generate the corresponding data sets: (a) $\pi=N(1,0.5^{2})$; (b) an exponential distribution of mean 2; (c) a mixture of normal distributions, $\pi=0.4\cdot N(0,0.5^{2})+0.6\cdot N(3,0.5^{2})$. For each case, we first train a VAE model with 5000 iterations, and then fine-tune the neural network parameter $\theta$ by running the following four training algorithms for additional 1000 iterations: (a) VAE; (b) the importance weighted autoencoders (IWAE, [@burda2015importance]) with $k=50$ importance samples; (c) Hamiltonian Monte Carlo (HMC) to approximate the true gradient; (d) SAGD. In HMC we use the same step size and chain length as SAGD, and run $L=5$ leapfrog steps to get each proposal. After training is finished, we simulate random variates of $Z$, and compare its empirical distribution $\hat{\pi}$ with the true latent distribution $\pi$. The Kolmogorov–Smirnov distance and 1-Wasserstein distance between $\hat{\pi}$ and $\pi$ are computed. Figure \[fig:latent\_estimation\] shows the data distribution, true latent distribution $\pi$, and the estimated $\hat{\pi}$ in each setting, based on one simulated data 