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

    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))


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)

will be

t_x = t_dist.ppf(q, t_nu, t_mu, np.sqrt(t_s))
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.