【統計入門】t分布・カイ二乗分布・F分布ってどんな分布なの？使い方も少しまとめた【python】

統計学習の備忘録として、記事を書きました。今回はt分布・カイ二乗分布・F分布についてまとめました。関連する信頼区間の推定とF検定もまとめました。

# t分布
母分散が未知でサイズの小さな標本から母平均などの推定・検定を行う場合に**標準正規分布の代わりに用いる**連続型の確率分布です。
t分布を説明する前にまず、標本平均$\bar x$に対してのz値についてみてます。

$$z_\bar x=\frac{\bar x -\mu}{\sqrt{\frac{\sigma^2}{n}}}$$
zは標準正規分布に従いますが、母分散$\sigma^2$がわかっていないと計算できません。現実的に母分散が明らかな場合は少ないです。そこで、この母分散$\sigma^2$を不偏分散${\hat{\sigma}}^2$で代用します。不偏分散とは、標本分散$s^2$の期待値が母分散に一致するように$\frac{n}{n-1}$をかけて調整した統計量です。$\frac{n}{n-1}$をかけるのは、nが十分に大きくないときは、標本分散の期待値は母分散よりも小さくなるためです。

$$t_\bar x=\frac{\bar x -\mu}{\sqrt{\frac{\hat\sigma^2}{n}}}=\frac{\bar x -\mu}{\sqrt{\frac{s^2}{n-1}}}$$
これがt値です。サイトなどによっては不偏分散で置き換えるのではなく、標本分散で置き換えているものもあります(下式)。

$$t_\bar x=\frac{\bar x -\mu}{\sqrt{\frac{s^2}{n}}}=\frac{\bar x -\mu}{\frac{s}{\sqrt{n}}}$$
ちなみに、$\frac{s}{\sqrt{n}}$は標本平均の標準誤差と言います。
そして、このt値は下記のような確率密度関数に従います。

$$f(t)=\frac{\Gamma((\nu+1) / 2)}{\sqrt{\nu \pi} \Gamma(\nu / 2)}\left(1+t^{2} / \nu\right)^{-(\nu+1) / 2}$$

$\nu$は自由度($n$-1)です。$\Gamma(\bullet)$はガンマ関数で、複素数を使った階乗を表す関数です。t分布は $\nu$によって形が変化します。

ここで、分布を実際に見てみましょう。

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

x = np.linspace(-6, 6, 1000)
fig, ax = plt.subplots(1,1, figsize=(10,7))

linestyles = [':', '--', '-.', '-']
deg_of_freedom = [1, 5, 30]
for df, ls in zip(deg_of_freedom, linestyles):
    ax.plot(x, stats.t.pdf(x, df), linestyle=ls, label=f'df={df}')

ax.plot(x, stats.norm.pdf(x, 0, 1), linestyle='-', label='Standard Normal Distribution')

plt.xlim(-6, 6)
plt.ylim(0, 0.4)
plt.title('t-distribution')
plt.legend()
plt.savefig('t-distribution.png')
plt.show()
```

![t-distribution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187535/d4a69808-8794-dc2a-6217-1fefc5a7ecb3.png)

自由度によって分布の形が変化していることがわかります。そして、標本サイズが大きいほど標準正規分布に近づいていることがわかります。

## 平均値の信頼区間の推定

この分布を使って母平均$\mu$が入ることを期待する信頼区間を推定してみます。

下記のようなデータが取れたとします。

```python
data = [np.random.randn() for _ in range(10)]
print(data)
# >>> [-0.14917153222917484, 0.7951720064790415, 0.662152983830839, 0.430521357874449, -2.48235088848113, 0.6166315938744059, 1.055076432212844, 0.7400193126962409, 0.90477126838906, -0.10509107744284621]
print(f"標本平均:{np.mean(data)}")
# >>> 標本平均:0.24677314572037296
print(f"標本分散:{np.var(data)}")
# >>> 標本分散:0.9702146524752354
print(f"標本標準偏差:{np.sqrt(np.var(data))}")
# >>> 標本標準偏差:0.9849947474353533
```

このデータは標準正規分布に基づいて生成されているので、真の平均$\mu$は0になります。自由度が9、95%信頼区間の推定は、下記のような式変形で求めることができます。

$$-2.262≤t_{\bar x}≤2.262$$

$$-2.262≤\frac{\bar x-\mu}{\frac{s}{\sqrt{n-1}}}≤2.262$$

$$\bar x - 2.262\frac{s}{\sqrt{n-1}}≤\mu≤\bar x + 2.262\frac{s}{\sqrt{n-1}}$$

今回はデータが10と少ないので、不偏分散置換型で計算しました。また、2.262という数はt分布表から自由度9、2.5%の交差点の値を取ってきます。なので、標本平均が上記のデータは95%の確率で母平均が

```python
bottom = np.mean(data) - 2.262*(np.sqrt(np.var(data))/(np.sqrt(len(data)-1)))
up = np.mean(data) + 2.262*(np.sqrt(np.var(data))/(np.sqrt(len(data)-1)))
print(f'{bottom} ≤ μ ≤ {up}')
# >>> -0.4959128938458835 ≤ μ ≤ 0.9894591852866295
```

ライブラリを使うと

```python
bottom, up = stats.t.interval(alpha=0.95, loc=np.mean(data), scale=np.sqrt(np.var(data)/(len(data)-1)), df=len(data)-1)
print(f'{bottom} ≤ μ ≤ {up}')
# >>> -0.49596449533733994 ≤ μ ≤ 0.9895107867780859
```
のようになると推定できます。


# カイ二乗分布
標本分散の分布など、**複数の変数を一度に扱える**のがカイ二乗分布です。

まずカイ二乗値の定義を下記に示します。

$$\chi^2_{(n)}\equiv\sum^n_{i=1}z^2_i=\frac{\sum^n_{i=1}(x_i-\mu)^2}{\sigma^2}$$

$\chi^2$の下付きの$n$は自由度を表し母平均を基準としているため$自由度=n$としています。自由度が大きくなるとカイ二乗値も大きくなりやすくなっています。

そして、カイ二乗分布の確率密度関数は以下のようになります。

$$f(x ; k)=\frac{1}{2^{k / 2} \Gamma(k / 2)} x^{k / 2-1} e^{-x / 2}$$

分布は以下のようになります。

```python
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
x = np.linspace(0, 10, 1000)
fig,ax = plt.subplots(1,1, figsize=(10,7))
 
linestyles = [':', '--', '-.', '-']
deg_of_freedom = [1, 3, 5, 10]
for df, ls in zip(deg_of_freedom, linestyles):
    ax.plot(x, stats.chi2.pdf(x, df), linestyle=ls, label=f'df={df}')

plt.xlim(0, 10)
plt.ylim(0, 1.0)
 
plt.title('chi 2 distribution')
plt.legend()
plt.savefig('chi2distribution.png')
plt.show()
```

![chi2distribution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187535/5ea3f002-0843-8e83-2732-1b1caac45a73.png)

## 母分散の区間推定
上記の式は母平均$\mu$だったので、標本平均$\bar x$で置き換えてみます。

$$\chi^2_{(n-1)}=\frac{\sum^n_{i=1}(x_i-\bar x)^2}{\sigma^2}$$

自由度はn-1しました。ここで、不偏分散$\hat \sigma^2$を使って表すと、

$$\chi^2_{(n-1)}=\frac{(n-1)\hat \sigma^2}{\sigma^2}$$

となり、$\chi^2$値が不偏分散と比例関係にあることがわかります。なので、母分散について解くと、

$$\sigma^2=\frac{(n-1)\hat \sigma^2}{\chi^2_{(n-1)}}$$

となるので、母分散の信頼区間を推定することができます。

データはt分布の時のものを使います。

```python
data = [np.random.randn() for _ in range(10)]
print(data)
# >>> [-0.14917153222917484, 0.7951720064790415, 0.662152983830839, 0.430521357874449, -2.48235088848113, 0.6166315938744059, 1.055076432212844, 0.7400193126962409, 0.90477126838906, -0.10509107744284621]
print(f"標本平均:{np.mean(data)}")
# >>> 標本平均:0.24677314572037296
print(f"標本分散:{np.var(data)}")
# >>> 標本分散:0.9702146524752354
print(f"標本標準偏差:{np.sqrt(np.var(data))}")
# >>> 標本標準偏差:0.9849947474353533
```

自由度が9の時、信頼区間は以下のように計算できます。

$$\frac{(n-1)\hat \sigma^2}{\chi^2_{(n-1,\alpha/2)}}≤\sigma^2≤\frac{(n-1)\hat \sigma^2}{\chi^2_{(n-1,1-\alpha/2)}}$$

$$\frac{(n-1)\hat \sigma^2}{19.02}≤\sigma^2≤\frac{(n-1)\hat \sigma^2}{2.7}$$

よって、

```python
bottom = ((len(data)-1)*np.var(data, ddof=1))/19.02
up = ((len(data)-1)*np.var(data, ddof=1))/2.7
print(f'{bottom} ≤ σ^2 ≤ {up}')
# >>> 0.5101023409438672 ≤ σ^2 ≤ 3.593387601760131
```

`np.var(data, ddof=1))`の`ddof=1`とすることで不偏分散を求めています。
ライブラリを使うと、カイ二乗分布表の値を取得できます。

```python
chi2_025, chi2_975 = stats.chi2.interval(alpha=0.95, df=len(data)-1)
bottom = ((len(data)-1)*np.var(data, ddof=1))/chi2_975
up = ((len(data)-1)*np.var(data, ddof=1))/chi2_025
print(f'{bottom} ≤ σ^2 ≤ {up}')
# >>> 0.5100281214306344 ≤ σ^2 ≤ 3.5928692971228506
```



# F分布
F分布は、2つの母集団から無作為抽出した**2つの標本に基づいた統計量が従う分布**です。この性質を利用して、抽出元である2つの母集団の分散が同じかどうかを検定する際などに用いられます。
さて、F値に関してですが、F値は、**正規分布に従う2つの母集団から無作為抽出してきた2つの$\chi^2$値の比**です。正規分布に従うという点は重要です。

$$F_{(\nu_1, \nu_2)}=\frac{\chi^2_{(\nu_1)}/\nu_1}{\chi^2_{(\nu_2)}/\nu_2}$$

また、確率密度関数は以下のようになります。

$$f\left(x ; k_{1}, k_{2}\right)=\frac{\Gamma\left(\frac{k_{1}+k_{2}}{2}\right) x^{\frac{k_{1}-2}{2}}}{\Gamma\left(\frac{k_{1}}{2}\right) \Gamma\left(\frac{k_{2}}{2}\right)\left(1+\frac{k_{1}}{k_{2}} x\right)^{\frac{k_{1}+k_{2}}{2}}}\left(\frac{k_{1}}{k_{2}}\right)^{\frac{k_{1}}{2}}$$

分布は以下のようになります。

![fdistribution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187535/da2492c0-afdd-c436-c47b-01e8e2772edb.png)

## F検定
F検定は2グループの分散に差があることを検定します。2グループの母分散($\sigma^2_1, \sigma^2_2$)が同じ場合F値は以下のようになります。

$$F=\frac{\chi^2_{(\nu_1)}/\nu_1}{\chi^2_{(\nu_2)}/\nu_2}=\frac{\frac{\nu_1\hat \sigma^2_1}{\sigma^2_1}/\nu_1}{\frac{\nu_2\hat \sigma^2_2}{\sigma^2_2}/\nu_2}=\frac{\hat \sigma^2_1}{\hat \sigma^2_2}$$

母分散が同じところから標本を抽出しているので、分散に差がなければ1に近づくはずです。逆にF値が1より大きい場合(普通分散が大きい方を分子とする)は、母分散は異なる可能性が高いです。

実際に検定してみます。
データは下記のようなデータを生成してみます。

```python
# 日本人男性10人をサンプル
np.random.seed(1)
Japan = np.round([np.random.normal(64, 9, 10)],1).reshape(10)
jp_var = np.var(Japan, ddof=1)
# アメリカ人男性10人をサンプル
np.random.seed(1)
US = np.round([np.random.normal(87, 12, 10)],1).reshape(10)
us_var = np.var(US, ddof=1)
print(f'日本人の不偏分散:{jp_var}')
# >>> 日本人の標本分散:127.71955555555557
print(f'アメリカ人の不偏分散:{us_var}')
# >>> アメリカ人の標本分散:226.57377777777785
print(f'F値:{us_var/jp_var}')
# >>> F値:1.7739944113665709
```

F検定は
帰無仮説：2グループの分散は等しい
対立仮説：2グループの分散に差がある
とします。自由度(9,9)のF分布で、下記のグラフの黄色い面積がp値にあたり、p値が0.05以下なら帰無仮説を棄却します。

```python
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
x = np.linspace(0.000001, 8, 1000)
fig,ax = plt.subplots(1,1, figsize=(10,7))
 
df = (9,9)
ls = '-'
y = stats.f.pdf(x, df[0], df[1])
ax.plot(x, y, linestyle=ls, label=f'k = {df[0]}, {df[1]}')

plt.xlim(0, 8)
plt.ylim(0, 1.0)
plt.fill_between(x, y, 0, where=x>=us_var/jp_var, facecolor='y',alpha=0.5)
plt.title('F distribution')
plt.legend()
plt.savefig('fdistribution_p.png')
print(f'p-value:{stats.f.sf(us_var/jp_var, len(Japan)-1, len(US)-1)}')
# >>> p-value:0.20301975133837194
plt.show()
```

![fdistribution_p.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187535/9962097a-4182-2401-2a95-3aea36d17342.png)

p値が0.203...で0.05より大きく帰無仮説は棄却できないので、2グループの分散に差があるとは言えなそうです。

# 参考
- [自由度](https://blog.apar.jp/data-analysis/5952/)
- [標本分散と不偏分散](https://bellcurve.jp/statistics/course/8614.html)
- [標準誤差](https://bellcurve.jp/statistics/course/8616.html)
- [ガンマ関数](https://ja.wikipedia.org/wiki/ガンマ関数)
- [t分布表](https://bellcurve.jp/statistics/course/8970.html)
- [カイ二乗分布表](https://bellcurve.jp/statistics/course/9210.html)