¿Función gamma inversa?

Question

Esta es una pregunta de análisis que recuerdo haber pensado en el instituto. Leer algunos de los otros temas aquí me recordó esto, y me gustaría escuchar las soluciones de otras personas a esto.

Tenemos la función gamma, que tiene una forma bastante elemental como todos sabemos,

$\Gamma(z) = \int_0^\infty e^{-t} t^{z-1} dt = \int_0^1 \left[ \ln(t^{-1}) \right]^{z-1}$

Lo que satisface por supuesto, $\Gamma(n) = (n-1)!$ , $n\in \mathbb{N}$ y las diversas relaciones de recurrencia y otras identidades que todos podemos buscar en wikipedia o mathwolrd o donde sea. Observamos que la función gamma es creciente en el intervalo $[a,\infty]$ donde $a\approx 1.46163$ .

La pregunta es: ¿podemos encontrar una función inversa explícita a la función gamma en este intervalo que tenga un aspecto igualmente sencillo?

Mis técnicas en ese momento consistían en escribir una ecuación diferencial que la inversa satisficiera, y resolverla, lo que podía hacer en términos de una expansión de series de potencias (estando en la escuela secundaria, ignorando los problemas de convergencia) para obtener una solución aproximada. Pero nunca fui capaz de obtener una solución muy bonita o exacta. Ahora tengo algunos trucos más sofisticados para hacer esto, pero me interesaría ver cómo la gente con más experiencia con este tipo de cuestiones iría a responder a esto.

La función gamma también satisface un número razonable de relaciones funcionales de aspecto algo interesante como $\Gamma(z)\Gamma(1-z)=\pi/\sin(\pi z)$ . ¿Cumple la función inversa alguna relación similar?

Answer 1

David Cantrell da una buena aproximación a $\Gamma^{-1}(n)$ en esta página .

Copiaré aquí el resultado por si alguna vez se cae esa página:

$k$ = el cero positivo de la función digamma, aproximadamente $1.461632$
$c$ = $\sqrt{2\pi}/e - \Gamma(k)$ , aproximadamente $0.036534$
$L(x)$ = $\ln(\frac{x+c}{\sqrt{2\pi}})$
$W(x)$ = Función Lambert W
$ApproxInvGamma(x)$ = $L(x)/W(\frac{L(x)}{e}) + \frac{1}{2}$

Answer 2

Para beneficio de las generaciones venideras añado aquí el código python que escribí después de leer las respuestas anteriores.

import numpy as np
import math
import scipy.special

def _lambert_w(z):
  """
  Lambert W function, principal branch.
  See http://en.wikipedia.org/wiki/Lambert_W_function
  Code taken from http://keithbriggs.info/software.html
  """
  eps=4.0e-16
  em1=0.3678794411714423215955237701614608
  assert z>=-em1, 'LambertW.py: bad argument %g, exiting.'%z
  if 0.0==z: 
      return 0.0
  if z<-em1+1e-4:
      q=z+em1
      r=math.sqrt(q)
      q2=q*q
      q3=q2*q
      return\
       -1.0\
       +2.331643981597124203363536062168*r\
       -1.812187885639363490240191647568*q\
       +1.936631114492359755363277457668*r*q\
       -2.353551201881614516821543561516*q2\
       +3.066858901050631912893148922704*r*q2\
       -4.175335600258177138854984177460*q3\
       +5.858023729874774148815053846119*r*q3\
       -8.401032217523977370984161688514*q3*q
  if z<1.0:
      p=math.sqrt(2.0*(2.7182818284590452353602874713526625*z+1.0))
      w=-1.0+p*(1.0+p*(-0.333333333333333333333+p*0.152777777777777777777777))
  else:
      w=math.log(z)
  if z>3.0: 
      w-=math.log(w)
  for i in xrange(10):
      e=math.exp(w)
      t=w*e-z
      p=w+1.0
      t/=e*p-0.5*(p+1.0)*t/p
      w-=t
      if abs(t)<eps*(1.0+abs(w)): 
          return w
  raise AssertionError, 'Unhandled value %1.2f'%z

def _gamma_inverse(x):
  """
  Inverse the gamma function.
  http://mathoverflow.net/questions/12828/inverse-gamma-function
  """
  k=1.461632 # the positive zero of the digamma function, scipy.special.psi
  assert x>=k, 'gamma(x) is strictly increasing for x >= k, k=%1.2f, x=%1.2f' % (k, x)
  C=math.sqrt(2*np.pi)/np.e - scipy.special.gamma(k) # approximately 0.036534
  L=np.log((x+C)/np.sqrt(2*np.pi))
  gamma_inv = 0.5+L/_lambert_w(L/np.e)
  return gamma_inv

Answer 3

Este artículo aborda la cuestión: http://www.ams.org/journals/proc/2012-140-04/S0002-9939-2011-11023-2/S0002-9939-2011-11023-2.pdf

Answer 4

Se puede obtener un mejor valor numérico del cero positivo de la función digamma utilizando la función incorporada de Mathematica FindInstance :

$k \approx 1.46163214496836234126265954233$

Entonces un mejor valor numérico de $c$ se encuentra y se utiliza:

$c=\sqrt{2\pi}/e - \Gamma(k)$

$c \approx 0.0365338144849004166003358471688$

Entonces se puede utilizar la relación mencionada en otras respuestas:

$\Gamma_{approx}^{-1}(n)=L(x)/W(\frac{L(x)}{e}) + \frac{1}{2}$

Sin embargo, un valor numérico real más preciso de la función gamma inversa se encuentra utilizando directamente FindInstance :

invgamnum[n_] = 
  FindInstance[Gamma[x] == n, x, Reals, WorkingPrecision -> 30]

El valor numérico obtenido se fijó en 30 dígitos decimales; esto puede cambiarse a una precisión numérica arbitraria.

Si Reals no se añade en el código anterior, FindInstance suele dar una solución de valor complejo.

FindInstance se puede utilizar para encontrar más de una solución (si existen).

Por ejemplo, para $\Gamma(x) = e$ con Mathematica:

FindInstance[Gamma[x] == E, x, Reals, 2, WorkingPrecision -> 30]

Los dos valores numéricos obtenidos son:

$x=\Gamma^{-1}(e) \approx 0.328713113325854368520566502893$

$x=\Gamma^{-1}(e) \approx 3.31244088253916037027609875238$

Nótese que si se utilizara la relación que implica la función W de Lambert, se encontraría un valor numérico único menos preciso:

$\Gamma^{-1}(e) \approx 3.3111260003863245271459515483$

Answer 5

Mathematica tiene una función gamma inversa. Está en la página web de funciones especiales . Esto sugiere que el problema es, al menos, lo suficientemente sencillo como para ser implementado por ordenador.

Acabo de encontrar más material sobre la inversa de la función gamma incompleta regularizada de Mathematica. Hay descargas en el sitio con información también. Algunos de ellos son cuadernos de Mathematica y necesitan el reproductor que es gratuito para ser abiertos.

La información incluye ecuaciones diferenciales, representaciones mediante funciones equivalentes y representaciones en serie, entre otras cosas.