Correct shape rate parameterisation of gamma and inverse gamma with scipy.stats

simgamma

from pylab import *

import scipy.stats as stats

########################################################################################
#
#       QUICK ILLUSTRATION OF USING SCIPY.STATS FOR GAMMA AND INV-GAMMA
#       HOW TO DO SHAPE-RATE (alpha, beta) PARAMETERISATION CORRECTLY
#
########################################################################################



g_mean = lambda a,b: float(a)/float(b)
g_var = lambda a,b: float(a)/(float(b)**2)

def report_gamma(a,b,N):
  X = stats.gamma.rvs(a,scale=1./b,size=N)
  print "Reporting gamma results with \t alpha=%s, beta=%s " % (a,b)
  print "Theoretical mean: ", g_mean(a,b), "\tEmpirical mean: ", mean(X)
  print "Theoretical variance: ", g_var(a,b), "\tEmpirical variance: ", var(X)



ig_mean = lambda a,b: float(b)/(a - 1.)
ig_var = lambda a,b: (float(b)**2)/((a-1)**2 * (a-2))

def report_invgamma(a,b,N):
  X = stats.invgamma.rvs(a,scale=b,size=N)
  print "Reporting inv-gamma results with \t alpha=%s, beta=%s " % (a,b)
  print "Theoretical mean: ", ig_mean(a,b), "\tEmpirical mean: ", mean(X)
  print "Theoretical variance: ", ig_var(a,b), "\tEmpirical variance: ", var(X)


N = 100000
for p in range(4,10):
  report_gamma(p,p,N)
  print "\n"
  report_invgamma(p,p,N)
  print "\n\n"

# CONCLUSION:
# FOR GAMMA USE:        scale=1./Beta   to do rate parameterisation
# FOR INVGAMMA USE:     scale=Beta      to do rate parameterisation