Inverse gamma function? - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-22T02:48:49Z http://mathoverflow.net/feeds/question/12828 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/12828/inverse-gamma-function Inverse gamma function? jeremy 2010-01-24T09:10:14Z 2012-06-03T13:50:06Z <p>This is an analysis question I remember thinking about in high school. Reading some of the other topics here reminded me of this, and I'd like to hear other people's solutions to this.</p> <p>We have the gamma function, which has a fairly elementary form as we all know,</p> <p>$\Gamma(z) = \int_0^\infty e^{-t} t^{z-1} dt = \int_0^1 \left[ \ln(t^{-1}) \right]^{z-1}$</p> <p>Which satisfies of course, $\Gamma(n) = (n-1)!$, $n\in \mathbb{N}$, and the various recurrence relations and other identities that we can all look up on wikipedia or mathwolrd or wherever. We note that the gamma function is increasing on the interval $[a,\infty]$ where $a\approx 1.46163$.</p> <p>The question is--can we come up with an explicit inverse function to the gamma function on this interval which looks similarly simple?</p> <p>My techniques at the time were to write down a differential equation that the inverse would satisfy, and solve it, which I could do in terms of a power series expansion (being in high school, ignoring the issues of convergence) to get an approximate solution. But I was never able to get a very nice looking or exact solution. I have a few more sophisticated tricks now to do this, but I would be interested to see how people with more experience with these kinds of questions would go about answering this.</p> <p>The gamma function also satisfies a reasonable number of somewhat interesting looking functional relations like $\Gamma(z)\Gamma(1-z)=\pi/\sin(\pi z)$. Does the inverse function satisfy any similar relations?</p> http://mathoverflow.net/questions/12828/inverse-gamma-function/12858#12858 Answer by Kristal Cantwell for Inverse gamma function? Kristal Cantwell 2010-01-24T19:03:19Z 2010-01-24T19:17:00Z <p>Mathematica has an inverse gamma function. It is on the web page on <a href="http://reference.wolfram.com/mathematica/tutorial/SpecialFunctions.html" rel="nofollow">special functions</a>. This would suggest that the problem is at least simple enough for computer implementation. </p> <p>I have just found more material on the inverse of the regularized incomplete gamma function from Mathematica. There are downloads on the <a href="http://functions.wolfram.com/GammaBetaErf/InverseGammaRegularized/" rel="nofollow">site</a> with information as well. Some of them are Mathematica notebooks and need the player which is free to be opened.</p> <p>The information includes differential equations, representations through equivalent functions and series representations among other things. </p> http://mathoverflow.net/questions/12828/inverse-gamma-function/28977#28977 Answer by BlueRaja for Inverse gamma function? BlueRaja 2010-06-21T18:56:27Z 2010-06-21T19:34:19Z <p>David Cantrell gives a good approximation of $\Gamma^{-1}(n)$ on <a href="http://mathforum.org/kb/message.jspa?messageID=342551&amp;tstart=0" rel="nofollow">this page</a>.</p> <p>I'll copy the result here in case that page ever goes down:</p> <p><code> k = the positive zero of the digamma function, approximately 1.461632<br> c = $\sqrt{2\pi}/e - \Gamma(k)$, approximately 0.036534<br> L(x) = $\ln((x+c)/\sqrt{2\pi})$<br> W(x) = <a href="http://en.wikipedia.org/wiki/Lambert_W_function" rel="nofollow">Lambert W function</a><br> ApproxInvGamma(x) = L(x) / W(L(x) / e) + 1/2 </code></p> http://mathoverflow.net/questions/12828/inverse-gamma-function/98267#98267 Answer by Uri Cohen for Inverse gamma function? Uri Cohen 2012-05-29T11:54:13Z 2012-05-29T11:54:13Z <p>For the benefit of generations to come I add here the python code I wrote after reading the above answers.</p> <pre><code>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&gt;=-em1, 'LambertW.py: bad argument %g, exiting.'%z if 0.0==z: return 0.0 if z&lt;-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&lt;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&gt;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)&lt;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&gt;=k, 'gamma(x) is strictly increasing for x &gt;= 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 </code></pre> http://mathoverflow.net/questions/12828/inverse-gamma-function/98301#98301 Answer by Michael Renardy for Inverse gamma function? Michael Renardy 2012-05-29T19:43:07Z 2012-05-29T19:43:07Z <p>This article addresses the question: <a href="http://www.ams.org/journals/proc/2012-140-04/S0002-9939-2011-11023-2/S0002-9939-2011-11023-2.pdf" rel="nofollow">http://www.ams.org/journals/proc/2012-140-04/S0002-9939-2011-11023-2/S0002-9939-2011-11023-2.pdf</a></p> http://mathoverflow.net/questions/12828/inverse-gamma-function/98717#98717 Answer by K K for Inverse gamma function? K K 2012-06-03T13:50:06Z 2012-06-03T13:50:06Z <p>this works for Mathematica:</p> <p>c = 0.036534 l[x_] = Log[(x + c)/Sqrt[2*Pi]] aig[x_] = l[x]/(ProductLog[l[x]/E]) + 1/2</p>