Este tipo de fluctuación salvaje se debe a errores de redondeo en coma flotante en los cálculos.
La función de riesgo de un $\Gamma(a,1)$ con el parámetro de forma $a$ y el parámetro de escala $1$ , es igual a
$$H(x; a) = \frac{x^{a-1}\exp(-x)}{\int_x^\infty t^{a-1} \exp(-t) dt }.$$
El máximo solicitado en la pregunta es también el valor límite ya que $x\to\infty$ porque la función de riesgo en este caso es creciente.
Tanto el numerador como el denominador son funciones diferenciables que se aproximan a cero como $x$ aumenta, por lo que se aplica la Regla de L'Hopital, que nos dice que el valor límite de la relación es el límite de la relación de las derivadas:
$$\lim_{x\to\infty} H(x;a) = \lim_{x\to\infty}\frac{\exp(-x)\left((a-1)x^{a-2} - x^{a-1}\right)}{-x^{a-1}\exp(-x)} = \lim_{x\to\infty} 1 - \frac{a-1}{x} = 1.$$
Cuando la escala se cambia a $b$ el PDF debe ser dividido por $b$ para compensar (para mantener el área total igual a la unidad), lo que implica el valor límite de la función de riesgo para una distribución Gamma con parámetro de escala $b$ es $1/b$ .
Una forma mejor de calcular esta función de riesgo para grandes $x$ es utilizar los primeros términos de su expansión de Taylor en torno a $x=\infty$ :
$$H(x;a,b) \approx x^{-a} \left(\frac{(a-1) \left(\frac{1}{b}\right)^{-a-1}}{x^2}-\frac{(a-1) \left(\frac{1}{b}\right)^{-a}}{x}+\left(\frac{1}{b}\right)^{1-a}\right ) \left(\frac{x}{b}\right)^a.$$
Esto será extremadamente preciso mucho antes de $x$ es tan grande que el cálculo de la relación por fuerza bruta se rompe. Por supuesto, no es preciso para pequeño $x$ .