Given an easily-computable probability density function $f(x)$, what algorithm can we use to numerically approximate percentiles?
For instance, we might be looking for $x$ such that given $X \sim f(x)$, $p(X \le x) = 95\%$.
Here is a use case that motivates this:
We're estimating the shape $\alpha$ and inverse scale $\theta$ of a gamma distribution (from a sequence of samples). Via Wikipedia we find a conjugate prior of the following form, originally developed by Miller (1980):
$$ f(\alpha, \theta) \propto \frac{\theta^{v \alpha - 1} p^{\alpha - 1} e^{-s\theta}}{\Gamma(\alpha)^n} $$
This somewhat unwieldy function represents our belief about likely values for $\alpha$ and $\theta$. We are only interested in $\theta$, so we find its pdf by integration:
$$ g(\theta) \propto \int f(\alpha, \theta) d\alpha $$
Miller (1980) suggests that we need numerical integration to compute this function.
Now we can easily plot this $g(\theta$) function by repeatedly approximating it for many $\theta$ values. From this density plot we can visually glean the most probably region for $\theta$. So far so good.
But we also want to compute the $\theta$ values for the 5th and 95th percentiles, $G^{-1}(0.05)$ and $G^{-1}(0.95)$, where $G^{-1}(y)$ is the inverse cdf. So my question is, how do we approximate those given $g(\theta)$?
There must surely be a way to numerically approximate this by harnessing some well-known algorithm, but I cannot find a way to express $G^{-1}$ in terms of $g$. Am I missing some intermediate step to make it work? The best I can come up with is getting a lot of samples for $g$ and summing them.
References:
Miller, Robert B. "Bayesian Analysis of the Two-Parameter Gamma Distribution" Technometrics, 1980, 22(1), 65-69
I haven't done math in a while, so it's all rusty -- if there's mistakes or misphrasings in this question, please edit away!
