t分布が平均が同じ正規分布の無限混合分布であることをシミュレーションする

厳密には，t分布は，平均が同じ，分散がガンマ分布に従うスケールパラメータ$h$でスケールされた $\sigma^{2} / h$ での正規分布で hをすべての範囲（無限）について足し合わせた（積分した）混合分布になっている．

ということで，大量の正規分布を作りpdfの和を求めてt分布に近づくかをシミュレーションしてみる．

参考：

t分布の混合性

潜在変数と混合分布の考え方を使うとt分布を導出できる．

混合正規分(MoG)では，潜在変数がcategorical分布（試行回数1の多項分布）に従うとし，重みをcategorical分布の確率とした混合分布であった．

シミュレーション

以下のシミュ―ションをしてt分布に近づくか検証する．

1次元の標準正規分布$N(0, 1)$をベースにする．
自由度$\nu=1$として，$h \sim \text{Gamm}(\frac{\nu}{2}, \frac{\nu}{2})$でスケールパラメータをサンプリング
正規分布$N(0, 1/h)$を30000個生成し，その密度を足し合わせる．
ヒストグラムの密度化（幅×密度の合計が1になるように）を使って正規化する．
t分布$t(\mu = 0, \sigma^{2}=1,\nu = 1)$に近づくかをみる．

import scipy as sp
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

%matplotlib inline

def init_plot(axs):
    for ax in axs:
        ax.set_xlim(-5, 5)
        ax.set_xlabel("$x$", fontsize=20)
        ax.tick_params(axis='both', which='major', labelsize=18)
        ax.legend(fontsize=18)


fig = plt.figure(figsize=(30, 10))
gs = gridspec.GridSpec(1, 3)

distN = 30000
df = 1
x = sp.linspace(-100, 100, 10000)

ax0 = fig.add_subplot(gs[0, 0])
ax1 = fig.add_subplot(gs[0, 1])
ax2 = fig.add_subplot(gs[0, 2])

h = gamma_df(df).rvs(size=distN)
pdf = sp.zeros_like(x)
for i in sp.arange(distN):
    N = stats.norm(loc=0, scale=sp.sqrt(1/h[i]))
    pdf += N.pdf(x)

    if i % 1000 == 0:
        ax0.plot(x, N.pdf(x))

normed_pdf = ((pdf[:-1] + pdf[1:]) / 2 ) / sp.diff(x)
normed_pdf = normed_pdf / sp.sum((pdf[:-1] + pdf[1:]) / 2 )
print(sp.sum(normed_pdf))
print(sp.sum(normed_pdf * sp.diff(x)))
ax1.plot(x, stats.norm.pdf(x), label="N(0, 1)", linestyle="--")
ax1.plot((x[:-1]+x[1:])/2, normed_pdf, label=f"sum{distN} N(0, 1/h)")
ax1.plot(x, stats.t(df=df).pdf(x), label=f"t(df={df})")
ax1.set_ylabel("density", fontsize=20)

ax2 = fig.add_subplot(gs[0, 2])
ax2.plot(x, stats.norm.pdf(x), label="N(0, 1)", linestyle="--")
ax2.plot((x[:-1]+x[1:])/2, normed_pdf, label=f"sum{distN} N(0, 1/h)", alpha=.6)
ax2.plot(x, stats.t(df=df).pdf(x), label=f"t(df={df})")
ax2.set_ylim(1e-5, 1)
ax2.set_yscale("log")
ax2.set_ylabel("density(log scale)", fontsize=20)

init_plot([ax0,ax1,ax2])
plt.tight_layout()
plt.savefig("sim_mix_Norm.png")
plt.show()