ish to investigate the relationship between the
normal-inverse-gamma distribution and the inverse gamma distribution.

>>> import numpy as np
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng(527484872345)
>>> mu, lmbda, a, b = 0, 1, 20, 20
>>> norm_inv_gamma = stats.normal_inverse_gamma(mu, lmbda, a, b)
>>> inv_gamma = stats.invgamma(a, scale=b)

One approach is to compare the distribution of the `s2` elements of
random variates against the PDF of an inverse gamma distribution.

>>> _, s2 = norm_inv_gamma.rvs(size=10000, random_state=rng)
>>> bins = np.linspace(s2.min(), s2.max(), 50)
>>> plt.hist(s2, bins=bins, density=True, label='Frequency density')
>>> s2 = np.linspace(s2.min(), s2.max(), 300)
>>> plt.plot(s2, inv_gamma.pdf(s2), label='PDF')
>>> plt.xlabel(r'$\sigma^2$')
>>> plt.ylabel('Frequency density / PMF')
>>> plt.show()

Similarly, we can compare the marginal distribution of `s2` against
an inverse gamma distribution.

>>> from scipy.integrate import quad_vec
>>> from scipy import integrate
>>> s2 = np.linspace(0.5, 3, 6)
>>> res = quad_vec(lambda x: norm_inv_gamma.pdf(x, s2), -np.inf, np.inf)[0]
>>> np.allclose(res, inv_gamma.pdf(s2))
True

The sample mean is comparable to the mean of the distribution.

>>> x, s2 = norm_inv_gamma.rvs(size=10000, random_state=rng)
>>> x.mean(), s2.mean()
(np.float64(-0.005254750127304425), np.float64(1.050438111436508))
>>> norm_inv_gamma.mean()
(np.float64(0.0), np.float64(1.0526315789473684))

Similarly, for the variance:

>>> x.var(ddof=1), s2.var(ddof=1)
(np.float64(1.0546150578185023), np.float64(0.061829865266330754))
>>> norm_inv_gamma.var()
(np.float64(1.0526315789473684), np.float64(0.061557402277623886))

Nc