How do I express sigma square follows Inverse Gamma distribution?

1

Hi All
I found a model shown below in a book on Bayesian statistics.
x \sim N(\mu, \sigma^2)
\sigma^2\sim IG(0.01, 0.01)
\mu\sim N(5, 4\sigma^2)
(IG means InverseGamma distribution)
To write this model in PyMC5, I wrote the following which didn’t work.

with pm.Model() as model:
    sigma**2 = pm.InverseGamma('sigma', 0.01, 0.01)
    mu = pm.Normal('mu', mu=5., sigma=2*sigma)

    y_pred = pm.Normal('y_pred', mu=mu, sigma=sigma, observed=x0)

    start = pm.find_MAP()    
    step = pm.step_methods.Metropolis()

I don’t know how to express “\sigma^2 follows InverseGamma distribution”.
Does anyone know how to write this?
Or is it impossible to write this model in PyMC5?

2

It’s possible but doesn’t seem necessary in your case. Just transform the variable afterwards:

with pm.Model() as model:
    sigma_sq = pm.InverseGamma('sigma_sq', 0.01, 0.01)
    mu = pm.Normal('mu', mu=5., sigma=2*pm.math.sqrt(sigma_sq))

    y_pred = pm.Normal('y_pred', mu=mu, sigma=sigma, observed=x0)
1 Like
3

That’s great. So simple.
Thank you so much!

1 Like
4

You can also code the multiplier directly within the inverse gamma to make it a clean prior that doesn’t need to get post-processed.

sigma_sq = pm.InverseGamma('sigma_sq', 0.01, 0.04)
mu = pm.Normal('mu', mu=5., sigma=pm.math.sqrt(sigma_sq))

This relies on the fact that if X \sim \textrm{invGamma}(\alpha, \beta), then k \cdot X \sim \textrm{invGamma}(\alpha, k \cdot \beta) (first line on the Wikipedia page of related distributions).

Before using these vague Gamma priors, I’d strongly suggest reading section 4.3 of Andrew Gelman’s article (comment!) from Bayesian Analysis about why they are less vague than people think:

http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf

3 Likes