Inverse Gamma Distribution with Newton's method

Question

I want to generate Gamma random variables using the inverse transform method. For this purpose I want to derive the inverse of the CDF of Gamma using the Newton's method. This method may be not so efficient. But can anyone help with finding some references? or explain it to me here. Thank you in advance

Answer 1

If you write the CDF as $$ F(x) = \frac{1}{\Gamma(k)}\gamma\left(k, \frac{x}{\theta}\right) \tag{1} $$ you can find the inverse by calculating the inverse of the incomplete gamma function, and there are plenty of resources to do that numerically, e.g. in python you can use scipy.special.gammaincinv

EDIT

If you need to use Newton's method to find $x$ in

$$ \mu = \frac{1}{\Gamma(k)}\gamma\left(k, \frac{x}{\theta}\right) \tag{2} $$

for a fixed $\mu$, define the function

$$ f(x) = \frac{1}{\Gamma(k)}\gamma\left(k, \frac{x}{\theta}\right) - \mu \tag{3} $$

and note that solving the problem $f(x) = 0$ is exactly the same a solving the problem (2). All you need to is to find

$$ f'(x) = \frac{1}{\Gamma(k)} \frac{{\rm d}}{{\rm d}x}\gamma\left(k, \frac{x}{\theta}\right) = \frac{1}{\Gamma(k)} \frac{{\rm d}}{{\rm d}x}\int_0^{x/\theta}{\rm d}t~ t^{k-1} e^{-t} = \frac{1}{\theta \Gamma(k)} \left(\frac{x}{\theta}\right)^{k-1} e^{-x/\theta} $$

EDIT 2

And after this, you. can use Newton's method as

$$ x_{n + 1} = x_n - \frac{f(x_n)}{f'(x_n)} $$

the resulting value of $x$ after convergence is the solution to equation (2), ieg. the inverse of $F$