Trouble getting correct inverse-gamma samples from jax.random.gamma

[PaulScemama]

Hi!

I'm trying to get samples from an Inverse Gamma from jax.random.gamma. I do so by inverting the Gamma samples (i.e. taking $\frac{1}{\text{gamma samples}}$). Here is some code I ran:

shape = 2
scale = 2

### Jax Gamma
key = jax.random.PRNGKey(123)
gamma_samples = jax.random.gamma(key, shape, (100000,)) * scale
print(f"Gamma Sample mean: {gamma_samples.mean()}")
print(f"Gamma Analytical mean: {shape * scale}")
# Gamma Sample mean: 4.003810405731201
# Gamma Analytical mean: 4

## Jax InvGamma
invgamma_samples = 1 / gamma_samples
print(f"InvGamma Sample mean: {invgamma_samples.mean()}")
print(f"InvGamma Analytical mean: {scale / (shape - 1)}")
# InvGamma Sample mean: 0.4992263913154602
# InvGamma Analytical mean: 2.0

I'm using the analytical means for both gamma and inverse-gamma distributions where the latter requires $shape > 1$. As you can see, I am not getting near the correct mean for the Inverse-Gamma while I do for the Gamma. I'm sure I'm doing something wrong, but I can't seem to find the mistake.

Thanks in advance!

[jakevdp]

I think you're using different forms of the scale parameter in each expression. Try it this way instead:

key = jax.random.PRNGKey(123)
gamma_samples = jax.random.gamma(key, shape, (100000,)) / scale
print(f"Gamma Sample mean: {gamma_samples.mean()}")
print(f"Gamma Analytical mean: {shape / scale}")
# Gamma Sample mean: 1.0009526014328003
# Gamma Analytical mean: 1.0

## Jax InvGamma
invgamma_samples = 1 / gamma_samples
print(f"InvGamma Sample mean: {invgamma_samples.mean()}")
print(f"InvGamma Analytical mean: {scale / (shape - 1)}")
# InvGamma Sample mean: 1.9969054460525513
# InvGamma Analytical mean: 2.0

[jakevdp]

Sure, I was just using the name scale for the variable because you had used that name. You're probably right that I should have renamed the variable in my answer.

[jakevdp]

But I think it should be this, no?

print(f"InvGamma Analytical mean: {1 / (shape - 1) / scale}")
``

[PaulScemama]

Sure, I was just using the name scale for the variable because you had used that name. You're probably right that I should have renamed the variable in my answer.

Yes you're definitely right it would have been confusing to use a new variable name in your reply. I was more just thinking out loud :).

But I think it should be this, no?

print(f"InvGamma Analytical mean: {1 / (shape - 1) / scale}")
``

Wouldn't that just end up being the same as the "rate" parameterized mean? I.e.

1 / (shape-1) / scale = scale / (shape - 1)

which we discussed to be the "rate" parameterized mean although we use scale just as a poorly named variable (see the initial code block I posted).

[jakevdp]

I have to admit I'm not totally clear on the "scale" vs "rate" nomenclature convention, but whatever you name it, if you want the code as written to print the correct values for shape != 2, you need my version.

[PaulScemama]

Yeah you're definitely correct. I think I'm all mixed up with the terminology ha. Thank you for the correct way to get the inverse gamma samples though!

I think this is what I was confused about: from wolfram:

"the inverse gamma distribution with shape parameter $\alpha$ and scale parameter $\beta$ is the distribution followed by the inverse of a gamma distribution with shape parameter $\alpha$ and scale parameter $\color{red} 1 / \beta$."

I thought, instead of $\text{Gamma} [\alpha, \color{red} 1 / \beta]$ it was just $\text{Gamma} [\alpha, \beta]$. I think this is what you were trying to tell me initially.