Closes #3813

This implements reparameterized gradient for distributions.Gamma. The gradient is implemented by directly approximating the reparameterized gradient function dx/dalpha following Knowles (2015). The approximation is accurate to within 0.5% relative error for a wide range of alphas.

Derivation

First consider the beta variable. If x ~ Gamma(alpha, beta) then x / beta ~ Gamma(alpha, 1). Since division is already implemented in PyTorch, we can thus reduce our problem to computing a reparameterized gradient of a standard gamma x ~ Gamma(alpha) = Gamma(alpha, 1) wrt alpha.

This PR implements a function standard_gamma_grad(x, alpha) that directly approximates the reparameterized gradient defined (for any continuous univariate distribution) as

                d/dalpha cdf(x; alpha)     d/dalpha cdf(x; alpha)
dx / dalpha = - ---------------------- = - ----------------------
                  d/dx cdf(x; alpha)           pdf(x; alpha)

This definition is used in the unit tests in tests/test_distributions.py, which compute d/dalpha cdf(x;alpha) via finite difference of the scipy.stats.gamma.cdf() function.

The approximation is split into three regions:

• For small x we use a power series approximation of cdf(x, alpha).
  Until digamma() is implemented in PyTorch, we use a finite difference of lgamma().
• For large alpha we use the approximation

  standard_gamma_grad(x, alpha) = sqrt(x/alpha)
  

• For intermediate x,alpha we use a rational function approximation

  standard_gamma_grad(x, alpha) = exp(PQ(log(x / alpha), log(alpha)))
  

  where PQ(u,v) is a rational function of order 2 in u and 3 in v.

For complete derivation, see this Jupyter Notebook.