I want to sample from the inverse gamma function in Python, but it is implemented in scipy invgamma [ ^ 1] was designed so that one of the parameters could not be tampered with (the parameter β of the function that appears in the next section is fixed at β = 1: please let me know if anyone knows the reason). I made it after practice (both posting to python and Qiita) (including memorandum).
The inverse gamma function is a continuous probability distribution expressed in the following form.
f(x,α,β) = \frac{β^α}{Γ(α)}x^{-α-1}e^{\frac{-β}{x}} (x>0)\\
(However, Γ(α) = \int_{0}^{∞}x^{α-1}e^{-x}Gamma function represented by dx)
I'm not sure why scipy's invgamma is fixed at β = 1 (I don't know much about anything other than my field ...), but I use MCMC (Markov Chain Monte Carlo) to priori-distribute the variance of the normal distribution. I wrote this article this time because I want to use other than β = 1 when using it. If you want to know more about MCMC, please see @ pynomi's article [^ 2].
from scipy import stats
from scipy.special import gamma
import numpy as np
###Probability density function of inverse gamma distribution###
class invgamma(stats.rv_continuous):
def _pdf(self, x,alpha,beta):
px = (beta**alpha)/gamma(alpha)*x**(-alpha-1)*np.exp(-beta/x)
return px
###Sampling from inverse gamma function###
invgamma = invgamma(name="invgamma", a=0.0)
sample_from_invgamma = invgamma.rvs(size=1, alpha = 1, beta = 1.0)
Like this. (OS: Windows10, Python3.7, Development environment: Operation confirmed with Spyder, when the random number is fixed with β = 1, the sampling value matches the invgamma of scipy, so I think it is probably correct.)
In writing this article and code, I also referred to @ physics303's article [^ 3]. Thank you very much.
Recommended Posts