¿Función gamma inversa?

Question

Esta es una pregunta de análisis en la que recuerdo haber pensado en el instituto. La lectura de 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 cual 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 aumenta en el intervalo $[a, \infty ]$ donde $a \approx 1.46163$ .

La pregunta es ¿podemos llegar a una función inversa explícita a la función gamma en este intervalo que parezca similarmente simple?

Mis técnicas en ese momento eran escribir una ecuación diferencial que la inversa satisfaría, y resolverla, lo cual podía hacer en términos de una expansión de las series de potencia (estando en la escuela secundaria, ignorando las cuestiones de convergencia) para obtener una solución aproximada. Pero nunca fui capaz de obtener una solución exacta o de aspecto agradable. Tengo algunos trucos más sofisticados ahora para hacer esto, pero me interesaría ver cómo la gente con más experiencia en este tipo de preguntas podrí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)$ . ¿La función inversa satisface alguna relación similar?

Answer 1

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

Copiaré el resultado aquí en caso de que esa página se caiga:

$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 de Lambert W
$ApproxInvGamma(x)$ = $L(x)/W( \frac {L(x)}{e}) + \frac {1}{2}$

Answer 2

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 3

Para el beneficio de las generaciones venideras añado aquí el código pitón 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 4

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.

Answer 5

Esto debería funcionar para Mathematica:

c = 0.036534
l[x_] = Log[(x + c)/Sqrt[2*Pi]]
aig[x_] = l[x]/(ProductLog[l[x]/E]) + 1/2