逆ガンマ関数をpythonで実装してみた

#導入
Pythonで逆ガンマ関数からサンプリングを行いたいのに、scipyに実装されているinvgamma¹はなぜか片方のパラメータがいじれない仕様になっていたので(次の節で出てくる関数のパラメータβがβ=1で固定：理由を知っている人がいたら教えてください)、練習（pythonとQiitaへの投稿の両方）がてら作ってみました（備忘録込み）。

#逆ガンマ関数とは
逆ガンマ関数とは、以下の形で表される連続な確率分布です。

f(x,α,β)　=　\frac{β^α}{Γ(α)}x^{-α-1}e^{\frac{-β}{x}} 　　(x>0)\\
（ただし　Γ(α)　= \int_{0}^{∞}x^{α-1}e^{-x}dx　　で表されるガンマ関数)

scipyのinvgammaがなんでβ=1で固定なのかよく分かりませんが(自分の分野以外のことはたいして知らない...)、私はMCMC（マルコフ連鎖モンテカルロ法）で正規分布の分散の事前分布に使用する時に、β=1以外も使いたいので今回この記事を書きました。MCMCを詳しく知りたい方は@pynomiさんの記事²を見るといいと思います。

#実際に書いてみた

from scipy import stats
from scipy.special import gamma
import numpy as np

###逆ガンマ分布の確率密度関数###
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

###逆ガンマ関数からサンプリング###
invgamma = invgamma(name="invgamma", a=0.0)    
sample_from_invgamma = invgamma.rvs(size=1, alpha = 1, beta = 1.0) 

こんな感じ。
（OS:Windows10, Python3.7, 開発環境:Spyderで動作確認済み、β=1として乱数を固定したとき、scipyのinvgammaとサンプリングの値が一致したので、多分あってると思います。）

#引用
この記事とコードを書くに当たって、@physics303さんの記事³も参考にしました。ありがとうございます。

1. scipy.stats.invgamma ↩

2. pythonでギブス・サンプリングによるベイズ線形回帰を実装してみた ↩

3. Pythonで任意の確率密度関数からサンプリングを行う方法 ↩