目次
Python で統計学
\(t\) 分布
推定や仮説検定で重要になる \(t\) 分布の密度関数は， \begin{eqnarray*} f_n(t) = \frac{1}{\sqrt{n}B\left(\frac{n}{2}, \frac{1}{2}\right)} \left( 1 + \frac{t^2}{n} \right)^{-\frac{n+1}{2}} \end{eqnarray*} である．ここで \(n~(=1,2,\cdots)\) は自由度である．つまりデータサイズが \(n+1\) のときに，自由度が \(n\) になる．さらに，\(B(\alpha,\beta)\) はベータ関数であり， \begin{eqnarray*} B(\alpha, \beta) = \int_{0}^{1}x^{\alpha - 1} (1 - x)^{\beta - 1}dx, \quad \alpha > 0, ~ \beta > 0 \end{eqnarray*} で定義される．
正規分布 \(\mathcal{N}\left(\mu, \sigma^2\right)\) に従う母集団から \(n\) 個のサンプルを抽出したとき，変数 \(t = \frac{\bar{x}-\mu}{\sqrt{U^2/n}}\) は自由度 \(n-1\) の \(t\) 分布に従う．ここで，\(\bar{x}\) はサンプルから計算できる標本平均で，\(U^2\) もサンプルから計算できる不偏分散である． 自由度が \(n \geq 2\) のとき，\(t\) 分布の平均は \(0\) となる．自由度が \(n \geq 3\) のときに，分散は \(n/(n-2)\) となる．
ライブラリのインポート
Python で \(t\) 分布を利用するためには Scipy ライブラリをインポートするとよい．
import numpy as np
from scipy.stats import t
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')
確率密度関数
確率密度関数の値を取得するには t.pdf(x, df) を使用すると良い．例えば，自由度 5 の \(t\) 分布について，いくつか密度関数の値を計算してみよう．
t.pdf(x = 0, df = 5)
0.3796066898224944
t.pdf(x = 1, df = 5)
0.2196797973509806
t.pdf(x = 2, df = 5)
0.06509031032621646
確率密度関数をプロットしてみよう．
# 横軸用のnumpy配列を準備
x_step = 0.05
x_min = -5
x_max = x_step + 5
x = np.arange(x_min, x_max, x_step)
# 自由度5のt分布を計算
y5 = t.pdf(x, 5)
# グラフを描く
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(x, y5)
ax.set_xlim(-4.5, 4.5)
ax.set_ylim(0, 0.45)
ax.set_xticks(np.arange(-4, 5, 1))
ax.set_yticks(np.arange(0.0, 0.5, 0.1))
ax.grid()
# plt.savefig('stat-t-5.png', dpi=300, facecolor='white')
plt.show()
\(t\) 分布の密度関数は標準正規分布に近く，自由度が \(\infty\) になれば標準正規分布に一致することが知られている．このことを実際にプロットして確認しよう．
# 自由度1,2のt分布を計算
y1 = t.pdf(x, 1)
y2 = t.pdf(x, 2)
# 標準正規分布を計算
from scipy.stats import norm
yn = norm.pdf(x)
# グラフを描く
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(x, y1, label="n=1")
ax.plot(x, y2, label="n=2")
ax.plot(x, y5, label="n=5")
ax.plot(x, yn, label="normal")
ax.legend() # 凡例を表示
ax.set_xlim(-4.5, 4.5)
ax.set_ylim(0, 0.45)
ax.set_xticks(np.arange(-4, 5, 1))
ax.set_yticks(np.arange(0.0, 0.5, 0.1))
ax.grid()
# plt.savefig('stat-t-125n.png', dpi=300, facecolor='white')
plt.show()
上側 \(\alpha\) 点
統計学テキストの付録によく掲載されている自由度 \(n\) の \(t\) 分布の上側 \(\alpha\) 点の値を求めていみよう．この 上側 \(\alpha\) 点の値は下図の面積がちょうど \(\alpha\) になるような横軸の値 \(t_n(\alpha)\) を意味する．
上側 \(\alpha\) 点は t.isf(q, df) で求めることができる．例えばテキストの表にある「\(\alpha = 0.05\) の列」と「\(n = 5\) の行」にある「2.01505」という値は次のように求めることができる．
t.isf(q = 0.050, df = 5)
2.0150483726691575