1 from scipy import *
   2 from pylab import *
   3 
   4 import scipy.stats as stats
   5 import scipy.integrate as integrate
   6 
   7 
   8 rc('text',usetex=True)
   9 
  10 def stats_gamma_pdf(p,x):
  11     return stats.gamma.pdf(x,p[0],loc=p[1],scale=p[2])
  12 
  13 def my_gamma_pdf(x,a):
  14     return 1.0/my_gamma_int(a)*x**(a-1)*exp(-x)
  15 
  16 
  17 def my_gamma_int(x):
  18     return integrate.quad(lambda t:(t**(x-1)*exp(-t)),0,integrate.Inf)[0]
  19 
  20 
  21 # Gamma function test
  22 for i in range(5):
  23     print i,my_gamma_int(i+1),factorial(i)
  24 
  25 x=linspace(0.01,10,101)
  26 
  27 a=1.0
  28 g1=my_gamma_pdf(x,a)
  29 g2=stats_gamma_pdf([a,0.0,1.0],x)
  30 
  31 plot(x,g1,linewidth=3)
  32 
  33 a=2.0
  34 g1=my_gamma_pdf(x,a)
  35 g2=stats_gamma_pdf([a,0.0,1.0],x)
  36 
  37 plot(x,g1,linewidth=3)
  38 #plot(x,g2,'y:',linewidth=3)
  39 
  40 a=2.0
  41 g1=my_gamma_pdf(x-1,a)
  42 g2=stats_gamma_pdf([a,1.0,1.0],x)
  43 
  44 plot(x,g1,linewidth=3)
  45 #plot(x,g2,'y:',linewidth=3)
  46 axis([0,10,0,0.4])
  47 
  48 a=2.0
  49 g1=my_gamma_pdf(x/2,a)/2.0
  50 g2=stats_gamma_pdf([a,0.0,2.0],x)
  51 
  52 plot(x,g1,linewidth=3)
  53 #plot(x,g2,'y:',linewidth=3)
  54 axis([0,10,0,0.7])
  55 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')
  56 
  57 show()