details visible already in the early work by @opper. Plugging the transformation $f(\cdot)$ into gives $$\mathcal{L}(\lambda) = \int \phi(\epsilon) \log \frac{p(x,
  f(\epsilon,\lambda))|\operatorname*{det}_{J_f}(\epsilon,
  \lambda)|}{\phi(\epsilon)}\ \delta \epsilon,$$ where the integral is over the standard distribution that does not depend on $z$, and $|\operatorname*{det}_{J_f}(\epsilon, \lambda)|$ is the absolute value of the determinant of the Jacobian of $f(\epsilon, \lambda)$. Consequently, it can be replaced by a stochastic approximation $$\mathcal{L}(\lambda) = \frac{1}{M} \sum_{m=1}^M \log \frac{p(x,
  f(\epsilon_m,\lambda))|\operatorname*{det}_{J_f}(\epsilon_m,
  \lambda)|}{\phi(\epsilon_m)},$$ where $\epsilon_m$ is drawn from $\phi(\epsilon)$. We can now easily compute the gradients using the chain rule, by first differentiating $\log p(x, f(\epsilon, \lambda))$ w.r.t $z$ and then $z = f(\epsilon,
\lambda)$ w.r.t $\lambda$, resulting in

\_ () \_[m=1]{}\^M . \[eq:MCELBO\]

The combination of the standard distribution $\phi(\epsilon)$ and the transformation $f(\epsilon, \lambda)$ defines the approximation family. For example, $\phi(\epsilon) = \mathcal{N}(0, I)$ and $f(\epsilon, \lambda=\{\mu, L\}) = \mu + L \epsilon$ defines arbitrary Gaussian approximations [@opper; @dsvi], where $L$ is the Cholesky factor of the covariance. To create richer approximations we can concatenate multiple transformations; @advi combines the transformation above with a rich family of univariate transformations designed for different kinds of parameter constraints. For example, by using $f(\epsilon, \lambda=\{\mu, L\}) = \operatorname*{Softplus}(\mu + L \epsilon)$ we 