igamma.py

from scipy.stats import rv_continuous
from scipy.special import gammaln, gammaincinv, gammainc
from numpy import log,exp

class igamma_gen(rv_continuous):
  def _pdf(self, x, a, b):
      return exp(self._logpdf(x,a,b))
  def _logpdf(self, x, a, b):
      return a*log(b) - gammaln(a) -(a+1)*log(x) - b/x
  def _cdf(self, x, a, b):
      return 1.0-gammainc(a,b/x) # why is this different than the wiki?
  def _ppf(self, q, a, b):
      return b/gammaincinv(a,1-q)
  def _munp(self, n, a, b):
      args = (a,b)
      super(igamma_gen, self)._munp(self, n, *args)

igamma = igamma_gen(a=0.0, name='invgamma', longname="An inverted gamma",
          shapes = 'a,b', extradoc="""

Inverted gamma distribution

igamma.pdf(x,a,b) = b**a*x**(-a-1)/gamma(a) * exp(-b/x)
for x > 0, a > 0, b>0.
"""
)