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?
