Sampling instead of delta method Let's say I estimating some parameter vector $\theta$ with some estimator $\hat \theta$. I have that asymptotically, 
$\hat \theta \text{ }\dot  \sim \text{ } N(\theta, \Sigma)$
Now let's say I want to compute a confidence interval for $f(\theta)$, where $f$ is some differentiable function. Given a consistent estimator of $\Sigma$ defined as $\hat \Sigma$, I can create an asymptotic confidence interval using the delta method, i.e.
$\hat {f(\theta)} \text{ } \dot \sim  \text{ } N\left(f(\theta),\text{ } f'(\theta)^T \Sigma \text{ } f'(\theta) \right)$
For a variety of reasons, this may be suboptimal in my case. Mostly the issue that $f$ is very non-linear about $\theta$, relative to area covered by $N(\theta, \Sigma)$.
In my head, a much better solution to this is to actually sample from 
$N(\hat \theta, \hat \Sigma)$
and then plug directly into $f$. We can then form confidence intervals in exact same manner we normally form credible intervals from MCMC samples. 
Is this a well known practice? Is there reference or a name for this procedure? I realize it's not possible that I'm the first person to think of this, but I don't know where to start looking about references for this. 
 A: The procedure you describe can be used to compute the standard errors of $f(\hat{\theta})$. Let us assume that knowledge of $\theta$ is available and that we are able to sample from:
$$
\tilde\theta_1,...,\tilde\theta_M \sim N(\theta, \Sigma).
$$
Then an estimate for the bias of $f(\hat\theta)$ is:
$$
\frac{1}{M}\sum_{m=1}^{M} f(\tilde\theta_m) - f(\theta),
$$
and an estimate of the variance is given by:
$$
\hat\Sigma_f = \frac{1}{M-1}\sum_{m=1}^{M}
\left(f(\tilde\theta_m) - \bar f \right)
\left(f(\tilde\theta_m) - \bar f \right)^{T}.
$$
If we choose to believe in the asymptotic normality of $f(\hat\theta)$ then symmetric confidence intervals can be constructed based on $\hat\Sigma_f$ in the usual way. 
If a normality assumption is undesirable, it is possible to make a weaker assumption that:
$$
f(\hat\theta) - f(\theta) \approx^D f(\tilde\theta) - f(\hat\theta)
$$
Assume w.l.o.g. that $f$ is one dimensional and denote by $D_q$ the $q$th quantile of the distribution of $f(\tilde\theta) - f(\hat\theta)$. Under the assumption above, a valid confidence intervals at a level $\alpha$ for $f(\hat\theta)$ is given by:
$$
\left(
f(\hat\theta) - D_{1 - \alpha/2},
f(\hat\theta) - D_{\alpha / 2}
\right)
$$
About the first step where we assumed knowledge of $\theta$, in practice we will sample:
$$
\hat\theta_1,...,\hat\theta_M \sim N(\hat\theta, \Sigma).
$$
However, it is important to keep in mind that if $\hat\theta$ is far from $\theta$ or not normal, then the resulting confidence intervals may be inaccurate. 
