#acl AdminGroup:read,write EditorGroup:read All:read
#format wiki
#language en
#pragma section-numbers off


{{{#!latex
The Gamma distribution is a one parameter (shape $\alpha$) function defined as
\begin{equation}
f(x;\alpha)=\frac{1}{\Gamma(\alpha)} x^{\alpha-1}e^{-x}
\end{equation}
with the Gamma function $\Gamma(\alpha)$ which is $\Gamma(x+1)=x!$ for natural numbers and defined as
\begin{equation}
\Gamma(x+1)=\int_0^{\infty} t^x e^{-t} dt
\end{equation}
for $x>0$. The Gamma distribution in the SciPy stats subpackage provides a location and scaling parameter 
in addition to the shape parameter. The following code demonstrates the use of the function and includes 
an alternative implementation of the function.
}}}
{{attachment:gamma.png}}

{{{#!python
from scipy import *
from pylab import *

import scipy.stats as stats
import scipy.integrate as integrate


rc('text',usetex=True)

def stats_gamma_pdf(p,x):
    return stats.gamma.pdf(x,p[0],loc=p[1],scale=p[2])

def my_gamma_pdf(x,a):
    return 1.0/my_gamma_int(a)*x**(a-1)*exp(-x)


def my_gamma_int(x):
    return integrate.quad(lambda t:(t**(x-1)*exp(-t)),0,integrate.Inf)[0]


# Gamma function test
for i in range(5):
    print i,my_gamma_int(i+1),factorial(i)

x=linspace(0.01,10,101)

a=1.0
g1=my_gamma_pdf(x,a)
g2=stats_gamma_pdf([a,0.0,1.0],x)

plot(x,g1,linewidth=3)

a=2.0
g1=my_gamma_pdf(x,a)
g2=stats_gamma_pdf([a,0.0,1.0],x)

plot(x,g1,linewidth=3)
#plot(x,g2,'y:',linewidth=3)

a=2.0
g1=my_gamma_pdf(x-1,a)
g2=stats_gamma_pdf([a,1.0,1.0],x)

plot(x,g1,linewidth=3)
#plot(x,g2,'y:',linewidth=3)
axis([0,10,0,0.4])

a=2.0
g1=my_gamma_pdf(x/2,a)/2.0
g2=stats_gamma_pdf([a,0.0,2.0],x)

plot(x,g1,linewidth=3)
#plot(x,g2,'y:',linewidth=3)
axis([0,10,0,0.7])
legend(('gamma(x;alpha=1,loc=0,scale=1)','gamma(x;alpha=2,loc=0,scale=1)','gamma(x;alpha=2,loc=1,scale=1)','gamma(x;alpha=2,loc=0,scale=2)'),'upper right')

show()
}}}