まず重要な用語をまとめておきましょう:
手元にある標本から母集団の性質(母数)を推測することを推測統計といいます。 推測統計は推定と検定に分けられます。
「推定」は、母数の値に対し具体的な値を求める点推定、つまり一つの値で推定の結果を表すものと、区間推定、つまりある幅を持った区間で結果を表すものに分けられます。
また「検定」は母集団について述べた異なる立場の二つの主張(仮説)のうちどれを採用するかを決定するものです。例えば、日本の中学生のテスト得点が5年前と現在とで変化したという仮説と、変化していないという仮説がある場合、どちらを採択するかをきめるのが検定です。
母集団から抽出した標本を用いて、母数の点推定を行う手順についてみていきます。 ここで、標本のデータの個数のことをサンプルサイズ、 もしくは標本の大きさと言います。
例えば、10人の17歳男性の身長データがあるとします。これは日本人の17歳男性全体という『母集団』から抽出された『標本』とみなせます。そしてこのとき標本のデータの個数が10ですから、n = 10と書きます。ここでnはサンプルサイズを表す記号です。
import numpy as np
Height=np.array([165.2, 175.9, 161.7, 174.2, 172.1, 163.3, 170.6, 168.4, 171.3])
この標本から、「17歳の日本人男性全体の身長」(母集団)の平均(母数の一つ)を点推定してみましょう。 それには、次のようにします。ここでmeanは、平均値を求めるRの関数でしたね。
print(Height.mean())
まず用語を確認しましょう:
</A>
</CENTER>
ここで、標本統計量からどのように母数を推定するかを表にしておきましょう:
母数 | 推定量 | 関係する関数 | オプション |
---|---|---|---|
母平均 | 標本平均 | np.mean | |
母分散 | 不偏分散 | np.var | ddof=1 |
母標準偏差 | 不偏分散の正の平方根 | np.std | ddof=1 |
母相関係数 | 標本相関係数 | np.corrcoef | 対応する要素 |
# kadai 3-1 計算・解答用
標本から推定値を得る手順よりも大事なこと:
推定値を計算する元になる、標本に含まれる「個々のデータの値」の決定方法について見ていきましょう。
[参考: サイコロの出目の確率分布] 理論的には、サイコロの出目の確率分布はつぎのようになります:
サイコロの出る目 | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
確率 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 |
それに対し、シミュレーションにより、サイコロの出る目の度数分布は次のようになりました: (1) 6回サイコロを振った場合:
サイコロの目 | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
度数 | 0 | 2 | 1 | 1 | 0 | 2 |
相対度数 | 0 | 2/6 | 1/6 | 1/6 | 0 | 2/6 |
サイコロの目 | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|
度数 | 10,097 | 9,918 | 9,936 | 9,960 | 10,050 | 10,039 |
相対度数 | 1.0097 | 0.9918 | 0.9936 | 0.9960 | 1.0050 | 1.0039 |
正規分布(normal distribution, ガウス分布 Gaussian distributionともいう)は統計学で最もよく用いられる確率分布なので、ここで復習しておこう
正規分布の表現方法: N (μ,σ2) --- μは平均値、σ2は分散 (σは標準偏差)
</A>
</CENTER>
![]() |
![]() |
正規分布 N (0,1)のグラフ | 正規分布 N (1,1) のグラフを重ね書き |
![]() |
![]() |
![]() |
μ-σからμ+σ の範囲 | μ-2σからμ+2σ の範囲 | μ-3σからμ+3σの範囲 |
from scipy import stats as st
import matplotlib.pyplot as plt
x = np.linspace(-4, 4, 10000)
plt.plot(x, st.norm.pdf(x,0,1), label=r'N(%i,%i)$' % (0,1))
plt.axis([-4,4,-0.02,0.5])
plt.axhline(0,c='g',ls='-.')
plt.axvline(-2,c='b',ls='--')
plt.axvline(2,c='b',ls='--')
x = np.linspace(-2, 2, 10000)
plt.fill_between(x, stats.norm.pdf(x,0,1))
plt.show()
正規母集団から単純無作為抽出を行う場合、特にサンプルサイズが小さいと、 確率変数が正規分布に従っているかどうかは標本データのヒストグラムを見てもわからない(左図)。 しかし大きなサイズになれば、ほぼ正規分布に従うことわかる(右図).
import numpy.random as random
help(random.normal)
samplesSmall = random.normal(50, 10.0, 10) # N(50,10^2)に従う乱数を10個生成
samplesSmall
plt.hist(samplesSmall)
plt.show()
samplesLarge = random.normal(50, 10.0, 10000) # N(50,10^2)に従う乱数を10000個生成
plt.hist(samplesLarge)
plt.show()
標本分布とは、標本統計量(標本平均、標本分散など)に関する確率分布のこと
標本分布は標本における個々のデータの実現値を表した度数分布ではなく、標本統計量の確率分布
標本分布は母集団分布と標本統計量の種類、そしてサンプルサイズが決まると理論的(数学的)に導かれる。
実際のデータから作成されるわけではない
標本分布を調べることで、推定値の性質が分かる可能性がある
標本分布を調べるときの観点:
Pythonを用いて「経験的」に標本分布を求める方法---「経験的」とは「現実に得られたデータに基づく」という意味、 理論的な標本分布を「近似的」に再現
母集団からサンプルサイズnの標本を何度も繰り返し抽出し、そのたびに標本統計量の実現値を計算して記録する
これを実行するときの問題
仮定:母集団分布は$N(50,10^2)$、 サンプルサイズ n =10
samples = random.normal(50, 10.0, 10) # N(50,10^2)に従う乱数を10個生成
samples # データの確認
np.mean(samples) # 平均
この結果、母平均の推定値は 55.5
[疑問: 55.4716... なのに答は 55.5?] 注意: 55.4716... と表示されたとしても、答として、これをこのまま答えないこと。
有効数字の桁数を考慮して、55.5 と答える(有効数字は2桁、計算途中なので3桁とする)。
母平均は50なので、推定値との差は 5.5 --- この差が標本誤差(誤差)
一回きりの標本抽出で推定値がどれくらいの誤差をもつかは、いろいろ
先の操作を10,000回繰り返してみよう:
import numpy as np
import numpy.random as random
meanArray = np.array([ np.mean(random.normal(50, 10.0, 10)) for _ in range(10000)])
こうして得られた標本平均のヒストグラムと母集団の正規分布とを、次のコードにより重ねあわせてみると、かなり近い関係にあることが分かります:
from scipy import stats as st
import matplotlib.pyplot as plt
%matplotlib inline
plt.hist(meanArray,normed=True) # 頻度ではなく割合で表示
x = np.linspace(35, 65, 10000)
plt.plot(x, st.norm.pdf(x,50,10.0**0.5), ls='--',label=r'N(%i,%i)$' % (0,1))
plt.grid()
plt.show()
平均$\mu$、分散$\sigma^2$の正規分布 $N(\mu, \sigma^2)$に従う母集団からサンプルサイズ$n$の標本を抽出したとき、標本平均の標本分布は、理論的に$N(\mu, \sigma^2/n)$ に従います。
誤差の絶対値が5以下、すなわち母平均の推定値が45以上、55以下の範囲で得られた回数を調べてみましょう:
# 45以上55以下のものを1それ以外を0としたリストを作り、そこで0でない(つまり1)の個数を数える
np.count_nonzero([1 if 45 <= x <= 55 else 0 for x in meanArray])
全体の89%は本当の母数±5に収まっていることがわかります。理論的にはこのデータは$N(50,10^2)$に従うので、次章の知識を先取りしscipy.statsモジュールの関数cdfを用いた理論的な値と比較してみましょう:
from scipy.stats import norm
help(norm.cdf)
norm.cdf(55, loc=50, scale=10.0**0.5) - norm.cdf(45,loc=50,scale=10.0**0.5)
このように、 上位二桁まで一致していることが確かめられました。
平均$\mu$、分散$\sigma^2$の正規分布 $N(\mu,\sigma^2)$ に従う母集団からサンプルサイズ$n$の標本を抽出したとき、標本平均の標本分布は、理論的に$N(\mu,\sigma^2/n)$に従うことは、次のようにして実験的に確かめることができます。
import numpy.random as random
for k in [50,60,70,80,90,100,150,200,250,300]: # サンプルサイズを指定
meanArray = np.array( [np.mean(random.normal(50, 10.0, k))
for _ in range(10000)]) # N(50,10^2)に従うk個の乱数の平均を10,000個
print("サンプルサイズ=%d, 平均=%f, 標本分散=%f, 標本分散*サンプルサイズ=%f"
% (k,meanArray.mean(), meanArray.var(),meanArray.var()*k))
以上の結果から、サンプルサイズを変えても標本平均の平均が母平均にほぼ一致すること、また標本平均の分散は、母分散/サンプルサイズであることが確かめられました。
ある推定量の標本分布の平均が、推定しようとしている母数の値と一致するとき、その推定量は不偏性があるといいます。また、 不偏性がある推定量のことを不偏推定量と言います。 推定量の不偏性は「標本分布が母数の本当の値を中心として分布しているか」に対応した概念です。
「推定量の標本分布の標準偏差(へんさ)」を標準誤差(ごさ)と言います。 これは「標本分布の横への広がり」を表す値です。標準誤差が小さいということは、 標本分布の広がりが小さいことを表し、これは『ばらつきが少ない』ことを意味します。
例: $N(50,10^2)$の正規母集団から n=10 の標本を抽出したとき、 標本平均の標本分布は$N(50,10^2/n)$に従うので、標準誤差は$\sqrt{n}$、つまり(n=10なので) およそ3.2となります。 ということは、標本平均から母集団の平均を推定しようとすると、おおよそ68%の確率で46.8から53.2の間の値が得られますが、運が悪いと5以上の誤差が生じることがあることがわかります。
この例から、(1)母集団の分散(母分散)が大きければ大きいほど標本の平均値は母平均から離れた値を取りやすい(標準誤差が大きくなりやすい)、(2)標本のサイズが小さいほど標準誤差が大きい、ということがわかります。 言い換えれば、標準誤差を小さくするにはサンプルサイズnを大きくするのが有効な方法であると言えます。
一般に、どのような標本統計量でも、サンプルサイズが大きいほど標準誤差は小さくなります
中心極限定理: 母集団が正規分布でなくとも、標本平均の標本分布は「サンプルサイズnが大きい時」はほぼ正規分布になる
標本分布を求めるの項目ではn=10の標本を10,000個作り、その標本平均を求めた。 この課題ではn=100として10,000個の標本平均を求めてみよう。そしてn=10の場合のヒストグラムとn=100の場合のヒストグラムを比較し、標準誤差が小さくなっていることを確かめてみよ。
%matplotlib inline
import matplotlib.pyplot as plt
import numpy.random as random
# 計算・解答用
「平均」以外の母数のうち、分散と中央値の標本分布について述べます。
母分散の推定量としての標本分散と不偏分散の違いを実験で示しましょう。
import numpy as np
import numpy.random as random
SampleVars = []
UnbiasedVars = []
for i in range(10000):
Sample = random.normal(50, 10.0,10) # N(50, 10^2)に従う乱数10個
SampleVars = SampleVars+[Sample.var()] # 標本分散
UnbiasedVars = UnbiasedVars + [Sample.var(ddof=1.0)]
SampleVars = np.array(SampleVars)
UnbiasedVars = np.array(UnbiasedVars)
print(SampleVars.mean())
print(UnbiasedVars.mean())
平均的に母分散の推定値として不偏分散のほうが近い値であることがわかります。次にその標準偏差を求めてみましょう。
print(SampleVars.std())
print(UnbiasedVars.std())
不偏分散のほうが推定値のばらつきが大きいことがわかります。
plt.subplot(2,2,1)
plt.hist(SampleVars)
plt.title(u'Sample Variance')
plt.subplot(2,2,2)
plt.hist(UnbiasedVars)
plt.title(u'Unbiased Variance')
plt.show()
plt.subplot(2,2,1)
plt.hist(SampleVars,bins=50)
plt.title(u'Sample Variance')
plt.subplot(2,2,2)
plt.hist(UnbiasedVars,bins=50)
plt.title(u'Unbiased Variance')
plt.show()
実験: 標本中央値の標本分布を求め、標本平均と比較する(どちらが母平均の推定量としてよいか?)
import numpy.random as random
SampleAverage = []
SampleMedian = []
for i in range(10000):
Sample = random.normal(50, 10.0, 10) # N(50, 10^2)に従う乱数10個
SampleAverage = SampleAverage+[np.mean(Sample)] # 標本平均
SampleMedian = SampleMedian + [np.median(Sample)] # 中央値
SampleAverage = np.array(SampleAverage)
SampleMedian = np.array(SampleMedian)
print(SampleAverage.mean())
print(SampleMedian.mean())
標本中央値も標本平均も母平均とほぼ一致---母平均の推定量としての候補として優劣が付けられない
print(np.std(SampleAverage))
print(np.std(SampleMedian))
標本平均のほうが標準誤差が小さい---標本平均は母平均の不偏推定量の中で</font>標準誤差が最小
plt.subplot(2,2,1)
plt.hist(SampleAverage,bins=50)
plt.title(u'Sample Mean')
plt.subplot(2,2,2)
plt.hist(SampleMedian,bins=50)
plt.title(u'Sample Median')
plt.show()
目的 | 関数名とモジュール | 使い方 |
---|---|---|
一様分布に従う乱数(復習) | random.random(個数) | random.random(10) # [0,1)の乱数を10個 |
正規分布に従う乱数 | random.normal(平均,標準偏差,個数) | random.normal(50,10,10) # N(50, 10^2)に従う乱数10個 |
指定の範囲の数の生成 | np.linspace(最小値,最大値,個数) | np.linspace(35, 65, 10000) # [35,65]の範囲の数を10000個 |
正規分布の密度関数 | st.norm.pdf(xデータ,平均,標準偏差) | plt.plot(x, st.norm.pdf(x,50,10.0**0.5)) # $N(50,\sqrt{10})$の正規分布の関数描画 |
$N(50,10^2)$の正規母集団から$n = 20$の標本抽出を5,000回繰り返すことにより、標本平均の経験的な標本分布を求めよ。またこの経験的な標本分布から得られるグラフ(ヒストグラム)に、理論的に導かれる標本分布を重ねあわせて表示せよ。
この結果下のようなグラフが得られればよい:
# 計算・解答用
# 計算・解答用
# 計算・解答用