StanとPythonでベイズ統計モデリング　その3　Chapter6

アヒル本(StanとRでベイズ統計モデリング)のChapter6にPythonで取り組んでいきます。この章は丁寧に分布を解説していくものなので、内容の復習は飛ばします。おざなりにされそうな章ですが、自分でパラメータをいじって分布からサンプリングしてみると新しい発見もありますので挑戦してみるのをおすすめします。

今回はRどころかStanも出てきません(笑)

練習問題

練習問題(1)

ベルヌーイ分布 or カテゴリカル分布に従う乱数を生成せよ…ということですが、numpyには両方とも存在しません。二項分布からのサンプリングを一回だけ行うものがベルヌーイ分布であり、多項分布からのサンプリングを一回だけ行うものがカテゴリカル分布であるため、それで代用しました。

カテゴリカル分布についてはRの結果と対応させるために少し手間を加えます。

import numpy as np

random_bernoulli = np.random.binomial(1, 0.5, 100)

random_categorical = np.random.multinomial(1, [0.2, 0.3, 0.5], 100)
np.where(random_categorical == 1)[1]

練習問題(2)

テキストで推奨されている、ベータ分布、ディリクレ分布、ガンマ分布、2変量正規分布、コーシー分布を生成していきます。ついでに可視化します。

ベータ分布

alphaが大きいほど右側に偏り、全体的に値が大きいほど真ん中(0.5辺り)に集合していくことが分かりました。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

alpha = 0.5
beta = 0.5
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

alpha = 1 
beta = 0.5
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

alpha = 1 
beta = 1 
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

alpha = 5 
beta = 5 
r_beta = np.random.beta(alpha, beta, size=1000)
g = sns.distplot(r_beta, kde=False, label="alpha={0}, beta={1}".format(alpha, beta))

g.legend()
plt.show(g)

f:id:ajhjhaf:20170531172702p:plain

ディリクレ分布

トピックモデルで使われる分布です。ベータ分布の多次元拡張なので、alphaに対する挙動はベータ分布と同じです。

今回は3次元で行いました(可視化の限界のためです…)。可視化の方法は先行事例を参考にしています。

import matplotlib.tri as tri
import numpy as np

def plot_points(X, barycentric=True, border=True, **kwargs):
    '''
    Copy from https://gist.github.com/tboggs/8778945
    '''
    corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])
    triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
    if barycentric is True:
        X = X.dot(corners)
    plt.plot(X[:, 0], X[:, 1], 'k.', **kwargs)
    plt.axis('equal')
    plt.axis('off')
    plt.xlim(0, 1)
    plt.ylim(0, 0.75**0.5)
    if border is True:
        plt.triplot(triangle, linewidth=1)
        
fig = plt.figure()
ax = fig.add_subplot(2,2,1)
alpha=[3,3,3]
r_diri = np.random.dirichlet(alpha, size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)

ax = fig.add_subplot(2,2,2)
alpha=[0.3,0.3,0.3]
r_diri = np.random.dirichlet([0.3,0.3,0.3], size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)

ax = fig.add_subplot(2,2,3)
alpha=[3,0.3,0.3]
r_diri = np.random.dirichlet([3,0.3,0.3], size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)

ax = fig.add_subplot(2,2,4)
alpha = [1,1,1]
r_diri = np.random.dirichlet([1,1,1], size=1000)
ax.set_title("alpha={0}".format(alpha))
plot_points(r_diri)
plt.show()

f:id:ajhjhaf:20170531172713p:plain

ガンマ分布

ガンマ分布は離散値を取る分布です。numpyの場合はshapeとscaleの２つのパラメータをとるのですが、scaleは教科書の $\beta$ の逆数になります。というわけで

平均値は $shape * scale$
標準偏差は $\sqrt{shape}*scale$

になります。ご注意下さい。

shape = 5
scale = 5
r_gamma = np.random.gamma(shape, scale, size=1000)
g = sns.distplot(r_gamma, label="shape={0}, scale={1}".format(shape, scale), kde=False)

shape = 5
scale = 20 
r_gamma = np.random.gamma(shape, scale, size=1000)
g = sns.distplot(r_gamma, label="shape={0}, scale={1}".format(shape, scale), kde=False)

shape = 40 
scale = 5
r_gamma = np.random.gamma(shape, scale, size=1000)
g = sns.distplot(r_gamma, label="shape={0}, scale={1}".format(shape, scale), kde=False)

g.legend()
plt.show()

f:id:ajhjhaf:20170531174245p:plain

2変量正規分布

平均と分散にどう2変量が関わるか示す相関係数を入れ込みます。共分散行列の書き方はもっと綺麗なものがありそうですが。。

import seaborn as sns
import matplotlib.pyploy as plt
import numpy as np

mu1 = 5
mu2 = 5
sigma1 = 10
sigma2 = 40
rho = 0.6
mu = np.array([mu1, mu2])
cov = np.array([[sigma1**2, sigma1*sigma2*rho], [sigma1*sigma2*rho, sigma2**2]])
r_mnorm = np.random.multivariate_normal(mu, cov, size=1000)

g = sns.jointplot(r_mnorm[:, 0], r_mnorm[:, 1], kind="kde", n_levels=50)
plt.show(g)

f:id:ajhjhaf:20170531172724p:plain

コーシー分布

numpyではパラメータを選べない標準コーシー分布しか使えません。そのためscipyを利用しました。

なおここまでのnumpyで生成した他の分布も、rvs(random variates sampling?)メソッドを使えば同様に得ることが出来ます。ただnumpyとscipyで仕様が異なり、numpyの正規分布を由来とする乱数はnumpy.random.normalなのに、scipyはscipy.stats.normだったりします(解せぬ)。

コーシー分布は面白い分布ですね。パラメータがあるものの、平均値や分散をそこから推定できない。しかもめちゃめちゃ変なところに外れ値みたいに値が登場する性質を持っています。

import seaborn as sns
import matplotlib.pyploy as plt
from scipy.stats import cauchy

r_cauchy = cauchy.rvs(0, 1, 100)
sns.distplot(r_cauchy)
plt.show()

f:id:ajhjhaf:20170531172819p:plain

練習問題(3)

$x_1 - x_2 \sim Normal(\mu_1 - \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2})$ のような分布になります。これは以前に触れた部分ですね。

正規分布の線形結合で表される分布は、記事内のpdfを読めば全部同じように表されることが分かります。

y1 = np.random.normal(50, 20, 2000)
y2 = np.random.normal(20, 15, 2000)
diff = y1 - y2
sns.distplot(diff)
plt.show()

f:id:ajhjhaf:20170531172937p:plain

練習問題(4)

何か自分でテキストに載っていない分布を対象にすること。という問題でした。

今回は負の二項分布を取り扱ってみます。生物学の分野では、遺伝子発現量(RNA-seq or cDNA-seq)の推定モデル(DESeqやedgeR)でよく使われている分布です。負の二項分布はポアソン分布とガンマ分布の混合モデルです。ホクソエムさんの分かりやすい解説があるので、そちらをご覧ください。ポアソン分布の $\gamma$ が確率的に変動することで、過分散に対応できるようになっています。

負の二項分布の解説は英語版Wikipediaを見るのが一番いいと思います。以下そこからのざっくりとした抽出です。

「負の」の意味

負の二項分布は定義4パターン存在します。4つはそれぞれ

成功率pの試行をr回失敗(成功)するまでに、成功(失敗)した回数kの分布
成功率pの試行をr回失敗した時の、試行回数nの分布
成功率pの試行をk回成功した時の、失敗回数rの分布
成功率pの試行をk回成功した時の、試行回数nの分布

となります。4つのパターンがありますが、一番良く使われるのは1つ目のものです。

1つ目のパターンを式で表すと、

f:id:ajhjhaf:20170603170330p:plain:w300

となります。このパターンは試行が成功でも失敗でも式に大きな変化がなく、逆で紹介されている場合もあります。その時は逆側の事象が起こる確率が1-pになるように対応が変化するだけです。ところで、成功回数と失敗回数の和が試行回数です。つまり $n = k + r$ となります。これを使って整理すると

f:id:ajhjhaf:20170603171743p:plain:w300

となりました。ところで普通の二項分布の式は次のようになります。

f:id:ajhjhaf:20170603171236p:plain:w300

つまり負の二項分布だと二項分布と比べてトライアル $n$ の数が-1されていることになります。

平均と分散

負の二項分布の平均と分散は次のようになります。

$\mu = \frac{pr}{1-p}$

$var = \mu(1+\frac{\mu}{r})$

詳しい導出についてはこちらを参考にするのが適切だと思いました。日本語版Wikipediaの記事は、どこか間違っているような…。

実数への拡張

負の二項分布のパラメータ $r$ は正の実数へ拡張することが出来ます。「rが実数の『成功率pの試行をr回失敗(成功)するまでに、成功(失敗)した回数kの分布』ってなんやねん」と僕は思いました。おそらくここで肝心なのは、 $r$ を正の実数に拡張しても元々の定義へ矛盾せずに一般化できるということではないかと思います。この時ガンマ関数 $\Gamma()$ を利用して、

$NegativeBinomial(k|n, p) \sim \frac{\Gamma(n)}{k!\Gamma(n-k)}$

のように書くことが出来ます。

平均と形状パラメータによる分布の指定

負の二項分布のパラメータを指定する時に $n, p$ ではなく、平均パラメータ $\mu$ と形状パラメータ $\alpha$ を使って記述されることがあります。形状パラメータとは分散をスッキリさせるために $r$ の逆数を取って $\alpha = \frac{1}{r}$ としたものです。一応形状パラメータと書きましたが、他にも呼び方が色々あり、

shape parameter
dispersion parameter
clustering coefficient
heterogeneity
aggregation

のようにとりとめがありません。文字の使い方も適当らしく、論文によっては $\phi$ が使われているものもありました。

分布の可視化

可視化をする時にbinsを自分で指定しています。何も指定しないと最適なbinsをseabornが決めてくれるのですが、範囲が小数で分割された場合には離散値である負の二項分布の値がダブルカウントされるのか、物凄く変な分布になってしまいます。気をつけて下さい。

import numpy as np
from scipy.stats import nbinom
import seaborn as sns
import matplotlib.pyplot as plt

# r_nb = np.random.negative_binomial(10, 0.3, 10000)
r_nb = nbinom.rvs(10, 0.3, size=10000)

sns.distplot(r_nb, kde=False, bins=np.arange(0,100,5))
plt.show()

f:id:ajhjhaf:20170531193122p:plain

余談ですがedgeRの論文は事実上、RNA-seqの解析に負の二項分布をあてはめだだけの論文です。2ページしかありません。そのわりに1000回を超えれば場外ホームランとなる世界で、2500回も引用されている(Scopus調べ)ので、なんとも夢があります。

Easy to type

個人的な勉強の記録です。データ分析、可視化などをメイントピックとしています。