2018年7月5日木曜日

$t$分布について

標準正規分布$N(0,1)$に従う母集団から2個サンプリングする場合を考える。

2個のサンプリングから、母集団の平均値を推定する場合、たった2個なので、さすがにばらつきが大きくなる。

2個のサンプリングを何回も繰り返して、算出した標本平均をヒストグラムに表すと、次のようになる。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from numpy.random import *
from statistics import mean, median,variance,stdev

%matplotlib inline

sample_array = np.zeros(20000)
n = 2
for i in range(0, 20000):
    sample = randn(n)
    sample_array[i] = mean(sample)

plt.hist(sample_array, range=(-3,3), bins=30)

標本数を増やすと、標本平均は母集団の母平均に近づくような気がしませんか?

例えば、10個のサンプリングを同様に複数回行い、標本平均をヒストグラムにすると、次のようになります。

n = 10
for i in range(0, 20000):
    sample = randn(n)
    sample_array[i] = mean(sample)

plt.hist(sample_array, range=(-3,3), bins=30)

明らかに、分布の幅が小さくなっている。

つまり標本平均の分布の分散[~1]が小さくなるということ。(この分散は標本平均の分散のことで、母分散(母集団の分散)とは異なることに注意)

サンプル数が増えると、標本平均は母平均に近づきます。またこのときの標準偏差を標準誤差とよび、$\frac{\sigma}{\sqrt{N}}$に従うことがテキスト3-5-10(p.168)に書かれています。

これを確かめるために、サンプル数を2から15まで増やしながら、標本平均の平均mu、標準誤差se、$\frac{\sigma}{\sqrt{N}}$hseを20000回のシミュレーションによって算出します。

In [4]:
df = pd.DataFrame(index=range(2,15), columns=['mu', 'se', 'hse'])
sample_mean = np.zeros(20000)
for n in range(2, 15):
    for i in range(0, 20000):
        sample = randn(n)
        sample_mean[i] = mean(sample)
    df.mu[n] = mean(sample_mean)
    df.se[n] = np.std(sample_mean, ddof=1)
    df.hse[n] = 1/np.sqrt(n)

df
Out[4]:
mu se hse
2 0.00363491 0.712645 0.707107
3 -0.00478089 0.57523 0.57735
4 -0.00340509 0.494789 0.5
5 -0.00183572 0.449929 0.447214
6 0.00406338 0.403411 0.408248
7 -0.00645362 0.378052 0.377964
8 -0.000127742 0.353343 0.353553
9 0.00400127 0.335482 0.333333
10 -0.00253611 0.313103 0.316228
11 0.00112946 0.301144 0.301511
12 -0.000473661 0.288476 0.288675
13 -0.00236576 0.276256 0.27735
14 -0.00202092 0.267012 0.267261

標本平均の平均は母平均0に近づき、標準偏差(標準誤差)は$\frac{\sigma}{\sqrt{N}}$によって近似されることがわかります。

また、今回は、標準正規分布($\mu=0, \sigma^2=1$)を母集団としていましたが、この性質は平均、分散の値にかかわらず、成り立ちます。

たとえば、平均170、標準偏差10の正規分布からのサンプルについて、同様に行います。

In [11]:
for n in range(2, 15):
    for i in range(0, 20000):
        sample = normal(170,10,n)
        sample_mean[i] = mean(sample)
    df.mu[n] = mean(sample_mean)
    df.se[n] = np.std(sample_mean)
    df.hse[n] = 10/np.sqrt(n)  #ここでは母分散を用いています

df
Out[11]:
mu se hse
2 169.933 7.0887 7.07107
3 170.042 5.81997 5.7735
4 170.015 4.97683 5
5 169.944 4.45906 4.47214
6 169.984 4.04427 4.08248
7 169.968 3.77021 3.77964
8 169.978 3.51379 3.53553
9 170 3.35437 3.33333
10 170.037 3.18029 3.16228
11 170.012 3.01981 3.01511
12 170.002 2.87617 2.88675
13 170.016 2.77772 2.7735
14 170.003 2.66788 2.67261

正規分布を標準正規分布にするための標準化は $$Z=\frac{X-\mu}{\sigma}$$ としました。つまり、正規分布に従っていれば、平均$\mu$と分散$\sigma^2$がどのような値であっても、標準化によって分布を標準正規分布$N(0, 1)$に変換できる、ということです。

標本平均についても標準化と同様に $$T=\frac{\hat{\mu}-\mu}{\sigma/\sqrt{N}}=\frac{\bar{X}-\mu}{se}$$ とすると、$T$の分布は母集団の平均と分散の値にかかわらず、同じ分布となります。

さきほどと同様の平均170、標準偏差10の分布からのサンプリングであれば、 $$T=\frac{\hat{\mu}-\mu}{\sigma/\sqrt{N}}=\frac{\bar{X}-170}{10/\sqrt{N}}$$ によって算出した$T$の分布は$t$分布に従うこととなります。

ただし、実際には母平均と母分散は分からないです。母分散は標本から算出した標本分散$\hat{\sigma}^2$で近似できますが,母平均を標本平均$\hat{\mu}$で近似すると,$T$の値はすべて0になってしまいます.

ぼくは,統計を勉強するときに,この母集団と標本集団,そして標本集団から算出した統計量(確率変数)の関係がよくわからなくなることが多いです.

話が冗長になりますが,ここで$t$分布の動機について,少し話したいと思います.

$t$分布の動機

Wikipediaには$t$分布について

正規分布する母集団の平均と分散が未知で標本サイズが小さい場合に平均を推定する問題に利用される。

と説明されています1

最初にも説明しましたが,男子高校生の身長の平均を推定するときに,サンプル数が母集団に近い300万人のときの標本平均と2人のときの標本平均では,たとえ値が同じだったとしても,データとして同じ価値(信用度?)がある,とは考えずらいですよね.このときの,小数標本に対する不信感は,平均や分散がどういう値をとるか,という点よりも,標本数が少ないという点に対して,特に抱かれるものだと思います.

$t$分布を用いた$t$推定や$t$検定は,少ないサンプルから算出された推定量が,どの程度信用できるものなのか,を教えてくれます.

ポアンカレさんとパン屋さんの例

ここで,ひとつ寓話を2

ポアンカレさんはパン屋さんで毎朝お買い物をします.毎朝,1000gのパンを買いますが,最近いつものパンが少し軽いような気がします.

そこで,毎朝,パンを美味しくいただく前に,パンを計量することにしました.その結果を

$$X=(x_1, x_2, \ldots, x_n)=(965, 1002, \ldots, 989)$$

として記録しました.

$t$推定によるパンの重さの平均の予想

平均を予想する場合は$\hat{\mu}=\bar{X}=\frac{1}{n} \sum^{n}_{i=1}x_i$です. この平均がどの程度信頼できるかを調べるために,信頼区間を算出する話も出てきました.

$$t=\frac{\hat{\mu}-\mu}{\sigma/\sqrt{n}}$$

が$t$分布に従うことより,信頼区間を算出します.$t$の値は,母集団の平均と分散がどのような値であっても,自由度$n-1$の$t$分布に従いますので,母分散を標本分散$\hat{\sigma}^2$で近似し,95%の信用区間を求めるには,下の図の赤い領域の合計が5%になるような$-t, t$を求め, $$-t<\frac{\hat{\mu}-\mu}{\hat{\sigma}/\sqrt{n}}<t$$ を変形して,母平均$\mu$の区間を標本平均$\hat{\mu}$,標本標準偏差$\hat{\sigma}$およびサンプルサイズ$n$で表すことができます.(移行して,符号を変えるだけなので,もし時間があったら式変形して,テキストと見比べてみてください)

t分布3

$t$検定によるパンの重さの検証

$t$推定によって,ポアンカレさんは,「パンの重さの平均は950g」と予想しました.

でも,ポアンカレさんはまだパン屋を疑うわけにはいかない,と考えます.もっと慎重に考えましょう.

例えば,$H_0$:「パンの重さの平均(母平均)は1000gだ(1000g以上)」という仮説を立てて,この仮説がどうにも標本の分布と相反している場合に初めて$H_1$:「パンの重さの平均(母平均)は1000gではない(1000g以下)」という仮説を支持し,パン屋の悪事を警察に通報しようと決めます.

このときの$H_0$が帰無仮説で,$H_1$が対立仮説です.帰無仮説を棄却し,対立仮説を支持したいので,導きたい結論の反対の内容を帰無仮説にします.

帰無仮説の条件のもとでは,母集団を1000としていますので,$\mu$は未知の値ではないので,

$$t_i=\frac{x_i-\mu}{\hat{\sigma}/\sqrt{n}}$$

として,$T=\{t_i\}$の値を求めます.この$T$が自由度$n-1$の$t$分布に従うため,

$$t=\frac{\hat{\mu}-\mu}{\hat{\sigma}/\sqrt{n}}$$

の値が求められます.ここで求めた値(仮に$t^*$とする)が$t$分布のどこにあるか,で帰無仮説が正しそうかどうか,評価ができます.

例えば,$t^*$が$t$分布のちょうど真ん中に位置していれば帰無仮説は正しそうですが,$t^*$が$t$分布の外れのほう(特に左の端っこ)にあれば,どうも帰無仮説は正しくないだろう,と考えられます.

これが,$t$検定の流れです.(話を端折りましたが,詳細はテキストをご参照ください)

まとめ

これでめでたくポアンカレさんは警察に訴え,パン屋さんを懲らしめることができました.ヨカッタ.

それは,さておき,ふたたびWikipediaから$t$分布の説明を引用しますと,

t分布は、連続確率分布の一つであり、正規分布する母集団の平均と分散が未知で標本サイズが小さい場合に平均を推定する問題に利用される。また、 2つの平均値の差の統計的有意性を検討するt検定で利用される。

となりますが,多少イメージできましたでしょうか.

テキストを読まれるときには,$\mu$とか$\sigma$だとかが,母集団のパラメータなのか,標本から算出されるものなか,仮説の条件から与えられるものなのか,混乱してしまってないか,確認されながら読まれると,諸々整理しながら理解できるのではないか,と思います.

また,母集団について推定したいのか,根拠を求めて仮説を検証したいのか,という動機も念頭に置きながら読まれるのがよいでしょう.

回帰になると,まずは仮定したモデルにおいて変数に対する係数を推定し,その推定した係数が妥当か,そのモデルが妥当かを検証し,改めてモデルを更新して係数を推定しなおして…みたいに,手法をとっかえひっかえしながら,行ったり来たりしますんで,焦らず,ほどよいペースで読み進めていただければ,と思います. (ときに,概要をつかむために,深い理解はさておくのも大事です!)


  1. https://ja.wikipedia.org/wiki/T%E5%88%86%E5%B8%83

  2. 寓話ですが,こういう話が流布されて,にわかに信じられてしまうほど,ポアンカレさんは一風変わった方だったのでしょうか.http://www.netlorechase.net/entry/2016/09/24/164544

  3. http://www.biwako.shiga-u.ac.jp/sensei/mnaka/ut/tdistinvtab.html から画像を引用しています

$t$分布について

標準正規分布$N(0,1)$に従う母集団から2個サンプリングする場合を考える。 2個のサンプリングから、母集団の平均値を推定する場合、たった2個なので、さすがにばらつきが大きくなる。 2個のサンプリングを何回も繰り返して、算出した標本平均をヒストグラムに表すと、次の...