はしくれエンジニアもどきのメモ

情報・Web系技術・Englishの勉強メモ・備忘録です。

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

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

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

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

参考:

t分布の混合性

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

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

関連リンク:

cartman0.hatenablog.com

この考え方を拡張して,

  • 平均(位置)が同じ,
  • 分散がガンマ分布に従うスケールパラメータ$h$でスケールされた$\sigma^{2}/h$

正規分布をhの全範囲(無限)についてガンマ分布の確率で重みづけした混合分布がt分布になる.

式に直すと以下になる.


\begin{eqnarray}
Pr(h) = \text{Gam}_{h}\left[ \frac{\nu}{2}, \frac{\nu}{2} \right] \\\\
Pr(\vec{x} | h) = \text{Norm}_{\vec{x}}\left[ \vec{\mu}, \frac{\Sigma}{h} \right] \\\\
\end{eqnarray}

\begin{eqnarray}
Pr(\vec{x} )
= \int_{0}^{\infty} Pr(\vec{x}, h) dh
&=& \int_{0}^{\infty} Pr(\vec{x} | h)Pr(h) dh \\\\
&=& \int_{0}^{\infty} \text{Norm}_{\vec{x}}\left[ \vec{\mu}, \frac{\Sigma}{h} \right] \cdot \text{Gam}_{h}\left[ \frac{\nu}{2}, \frac{\nu}{2} \right] dh \\\\
&=& \int_{0}^{\infty} w_{h} \text{Norm}_{\vec{x}}\left[ \vec{\mu}, \frac{\Sigma}{h} \right]  dh \\\\
&=& \text{Stud}_{\vec{x}} \left[\vec{\mu}, \mathbf{\Sigma}, \nu \right]
\end{eqnarray}

この積分(周辺化)を計算するとt分布のpdfと同じになる.

$\nu$がt分布の自由度になっている.

関連リンク:

cartman0.hatenablog.com

シミュレーション

以下のシミュ―ションをして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()

出来上がる図は以下になる.

  • 左:混合対象の30000個中の30個の正規分布$N(0, 1/h)$
  • 中:混合して密度化した分布とt分布
  • 右:対数スケールに直したもの

真ん中の図で誤差が発生する理由は,以下である. ただ,誤差が小さためログスケールではほぼ一致しているように見える.

  • 混合する正規分布の数が無限でない.30000個だけ
  • 密度化の際に横軸をxを$(-\infty, \infty)$にしていない.$(-100, 100)$の範囲で近似している.
  • また,密度化の際にビン幅が大きい(分解能が低い)ため.

シミュレーション結果