# mismatch in sampling between t distribution and normal-inverse-gamma distribution

I am looking at equivalence of sampling between t distribution and normal-inverse-gamma (NIG) distribution in python. The results don't match, and I want to see if there's a mistake in how I am sampling or if it is a problem with scipy/numpy implementation of the same.

The results are expected from Bayesian linear regression, specifically for posterior predictive distribution.

To give you a background, under NIG prior, the posterior predictive distribution of $$y$$ can be obtained as follows. Posterior for $$\beta$$ will be a $$NIG(\beta, \sigma^2|\mu, V, a, b)$$. Also $$y$$ follows $$N(y|x\beta, \sigma^2)$$. Once you have samples of $$NIG(\beta, \sigma^2|\mu, V, a, b)$$, it is trivial to get samples of the predictions of $$y$$. This distribution of $$y$$ will be a $$t-distribution$$ with parameters $$t(y|\mu=x\beta, \frac{b}{a} (1 + xVx^T))$$.

But the samples from python are now matching.

from scipy.stats import invgamma, t as t_dist
from numpy.random import multivariate_normal

def normal_inverse_gamma(mu, cov, a, b, n=1):
"""
Generates n random samples from β where β ~ NIG(mu, V, a, b), i.e.:
β | σ^2 ~ N(mu, σ^2 * V)
σ^2 ~ IG(a, b)
"""
cov = 0.5 * (cov + cov.T)
s2 = invgamma.rvs(a, loc=0, scale=b, size=n)
beta = np.apply_along_axis(lambda x: multivariate_normal(mu, x * cov), 0, s2[None, :]).T
return s2, beta

def test_nig_samples_with_tdist():
"""
Compares the equivalence of predictions based on NIG samples, and t distribution
"""
print()

cov = np.array(
[
[2.658472e-01, -1.720448e-01, -2.694973e-02, -1.201772e-01, 4.615540e-17],
[-1.720448e-01, 1.173027e-01, 1.528254e-02, 7.190824e-02, -2.985954e-17],
[-2.694973e-02, 1.528254e-02, 6.707327e-03, 1.625543e-02, -4.724725e-18],
[-1.201772e-01, 7.190824e-02, 1.625543e-02, 6.329014e-02, -2.087379e-17],
[4.615540e-17, -2.985954e-17, -4.724725e-18, -2.087379e-17, 1.999996e-03],
]
)

cov = 0.5 * (cov + cov.T)
mu = np.array([-0.68469919, 0.02140081, 0.06961789, 0.05209994, 3.9761895])
a = 251.005
b = 368.6142568638883

x = np.array([0.2, 0.3, 0.6, 0.1, 1.0])

q = 0.2

# x at quantile q using t_dist
t_nu = 2 * a
t_mu = np.dot(x, mu)
t_s = (b / a) * (1 + np.dot(x, cov).dot(x))

t_x = t_dist.ppf(q, t_nu, t_mu, t_s)

print("based on t-distribution", t_x)

n_samples = 20000
sigma_2, beta = normal_inverse_gamma(mu, cov, a, b, n=n_samples)
t_samples = np.array(
[
np.random.normal(loc=np.dot(beta[i, :], x), scale=np.sqrt(sigma_2))
for i in range(n_samples)
]
)

print("based on NIG samples", np.percentile(t_samples, q=q * 100))

return


For q=0.2, the results are:

based on t-distribution 2.648105110347587
based on NIG samples 2.8693436388254305


so very different. The distribution based on $$NIG$$ is more peaked than $$t$$.

What's wrong here?

The problem was in my call to python for $$t-distribution$$.
t_x = t_dist.ppf(q, t_nu, t_mu, t_s)

t_x = t_dist.ppf(q, t_nu, t_mu, np.sqrt(t_s))