Plotting Normal Inverse Gamma Distirbution

Scipy stats doesnt have Normal Inverse Gamma distirbution.

We would like to incorporate Normal Inverse Gamma distirbution in “scipy.stats” package.

Learning about Normal Inverse Gamma(NIG) distribution will lead you to a plot like this from wikipedia. .

It was intruiging enough to find out how to plot this graph in python and was sure that there will be some already plots available. But to my suprise there is no blogs or docs to plot NIG in python. The closest I found was in R langugage in [1] by Frank Portman.

So I spent some time to plot NIG in python below is the snippet for it. Special thanks to Jake Vadendeplas[2] for his wonderful blogs about visualization in python.

Normal Inverse Gamma Distribution

Let the input \(x\) on which its modelled be : \[ x = [\mu, \sigma^2] \]

Probability density function (PDF)

\[ f(x | \delta, \alpha, \beta, \lambda ) = \sqrt{\left(\frac{\lambda}{2 \pi x[\sigma^2]} \right)} \frac{\beta^\alpha}{\Gamma(\alpha)} \left(\frac{1}{x[\sigma^2]} \right)^{(\alpha + 1)} \exp{ \left( - \frac{2\beta + \lambda \left(x[\mu] - \delta \right)^2 }{ 2 x[\sigma]^2} \right)} \]

from scipy.stats import rv_continuous
from scipy.stats import norm
from scipy.stats import gengamma
from scipy.special import gamma
from scipy.stats import expon
import numpy as np

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-white')

class norminvgamma():
  r"""A normal inverse gamma random variable.
    The mu (``mu``) keyword specifies the parmaeter mu.
    %(before_notes)s
    Notes
    -----
    The probability density function for `norminvgamma` is:
    .. math::
        x = [\mu, \sigma^2]
        f(x | \delta, \alpha, \beta, \lamda) = 
               \sqrt(\frac{\lamda}{2 \pi x[\sigma^2}])
               \frac{\beta^\alpha}{\gamma(\alpha)}
               \frac{1}{x[\sigma^2]}^(\alpha + 1)
               \exp(- \frac{2\beta + \lamda(x[\mu] - delta)^2}{2 x[\sigma^2] })
        
    for a real number :math:`x` and for positive number :math: `\sigma^2` > 0
    %(after_notes)s
    %(example)s
    """
  def __init__(self, delta, alpha, beta, lamda):
    self.argcheck(delta, alpha, beta, lamda)
    self.delta = delta
    self.alpha = alpha
    self.beta = beta
    self.lamda = lamda

  def argcheck(self, delta, alpha, beta, lamda):
        return (alpha > 0) 

  def rvs(self, size=1):
    sigma_2 = gengamma.rvs(self.alpha, self.beta,  size=size)
    sigma_2 = np.array(sigma_2)
    return [[norm.rvs(self.delta, s/self.lamda), s] for s in sigma_2]
  
  def pdf(self, xmu, xsigma2):
    t1 = ((self.lamda)**0.5) * ((self.beta)**self.alpha)
    t2 = (xsigma2 * (2 * 3.15)**0.5) * gamma(self.alpha)
    t3 = (1 / xsigma2**2)**(self.alpha + 1)
    t4 = expon.pdf((2*self.beta + self.lamda*(self.delta-xmu)**2)/(2*xsigma2**2))
    #print (t1, t2, t3, t4)
    return (t1/t2)*t3*t4
    

  def stats(self):
    #ToDo
    return

  def plot(self,zoom=0.9, axs=None):
    
    steps = 50
    max_sig_sq = gengamma.ppf(zoom, self.alpha, self.beta) * self.lamda
    #print(max_sig_sq)
    mu_range = np.linspace(self.delta - 1 * max_sig_sq, self.delta + 1 * max_sig_sq, num=steps)
    #print (mu_range[0], mu_range[-1])
    sigma_range = np.linspace(0.01, max_sig_sq, num=steps)
    mu_grid, sigma_grid = np.meshgrid(mu_range, sigma_range)
    pdf_mesh = self.pdf(mu_grid, sigma_grid)

    if axs:

      contours = axs.contour(mu_grid, sigma_grid, pdf_mesh,  20, cmap='RdGy');

      plt.clabel(contours, inline=True, fontsize=8)
      #extent=[mu_range[0], mu_range[-1], sigma_range[0], sigma_range[-1]]
      axs.imshow(pdf_mesh, extent=[mu_range[0], mu_range[-1], sigma_range[0], sigma_range[-1]],
                origin='lower', cmap='Blues', alpha=0.5)
      axs.axis('equal')
      axs.set_title("("+str(self.delta)+","+str(self.alpha)+"," 
                    +str(self.beta)+","+str(self.lamda)+")")
      #plt.colorbar();
    else:
      assert True, "Pass the axes to plot from matplotlib"

Varying different range of \(\alpha\)


#norminvgamma = norminvgamma_gen()
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(15,5))
fig.suptitle('Vertically alpha')
nig = norminvgamma(delta=0,alpha=1,beta=1, lamda=1)
samples = nig.rvs(size=10)
nig.plot(axs=axs[0])
nig = norminvgamma(delta=0,alpha=2,beta=1, lamda=1)
nig.plot(zoom=0.7, axs=axs[1])
nig = norminvgamma(delta=0,alpha=4,beta=1, lamda=1)

nig.plot(zoom=0.2, axs=axs[2])

Varying different range of \(\beta\)

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(15,5))
fig.suptitle('Varying beta')
nig = norminvgamma(delta=0,alpha=1,beta=1, lamda=1)
samples = nig.rvs(size=10)
nig.plot(axs=axs[0])
nig = norminvgamma(delta=0,alpha=1,beta=2, lamda=1)
nig.plot( axs=axs[1])
nig = norminvgamma(delta=0,alpha=1,beta=3, lamda=1)

nig.plot(axs=axs[2])

Varying different range of \(\lambda\)

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(15,5))
fig.suptitle('Verying lamda ')
nig = norminvgamma(delta=0,alpha=1,beta=1, lamda=1)
samples = nig.rvs(size=10)
nig.plot(axs=axs[0])
nig = norminvgamma(delta=0,alpha=1,beta=1, lamda=2)
nig.plot( axs=axs[1])
nig = norminvgamma(delta=0,alpha=1,beta=1, lamda=4)

nig.plot(axs=axs[2])

References

https://frankportman.github.io/bayesAB/reference/plotNormalInvGamma.html
https://jakevdp.github.io/PythonDataScienceHandbook/04.04-density-and-contour-plots.html