Pythonで統計学を学ぶ(3)

この内容は山田、杉澤、村井(2008)「R」によるやさしい統計学を参考にしています。

この講義では、「母集団と標本」をとりあげます。 これは、大きな集団から一部を取り出した少数のデータの情報を用いて、もとの集団の性質を推測する推測統計の基本的な理論の学習がねらいです。 特に 標本分布の理解を目的とします。

学習項目です

母集団と標本の関係

まず重要な用語をまとめておきましょう:

  • 母集団: 対象のデータ全体               
  • 標本: 母集団から一部のデータを取り出した(抽出した)もの                
  • 標本抽出: 母集団から標本をとりだすこと
  • 母数: 母集団の性質をあらわす統計的指数(比率、平均、分散、相関係数など) </UL> 例えば、日本の中学生全体のテスト得点データを母集団とすると、名古屋市昭和区の(一部の)中学生のテスト得点データが「標本」になります。そして、平均点や分散などが母数とみなせます。

推測統計

手元にある標本から母集団の性質(母数)を推測することを推測統計といいます。 推測統計は推定検定に分けられます。

「推定」は、母数の値に対し具体的な値を求める点推定、つまり一つの値で推定の結果を表すものと、区間推定、つまりある幅を持った区間で結果を表すものに分けられます。

また「検定」は母集団について述べた異なる立場の二つの主張(仮説)のうちどれを採用するかを決定するものです。例えば、日本の中学生のテスト得点が5年前と現在とで変化したという仮説と、変化していないという仮説がある場合、どちらを採択するかをきめるのが検定です。

点推定

母集団から抽出した標本を用いて、母数の点推定を行う手順についてみていきます。 ここで、標本のデータの個数のことをサンプルサイズ、 もしくは標本の大きさと言います。

例えば、10人の17歳男性の身長データがあるとします。これは日本人の17歳男性全体という『母集団』から抽出された『標本』とみなせます。そしてこのとき標本のデータの個数が10ですから、n = 10と書きます。ここでnはサンプルサイズを表す記号です。

In [2]:
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の関数でしたね。

In [3]:
print(Height.mean())
169.188888889

推定の基礎

まず用語を確認しましょう:

  • 標本統計量:標本データから計算される平均や分散などの値
  • 推定値:標本データを用いて計算した母数の推定量
母集団、標本、標本統計量、母数の関係を図で表すと以下のようになります: 母集団から抽出された標本データを分析し、標本統計量を得ます、それを使って母数を推定するのです

</A> </CENTER> ここで、標本統計量からどのように母数を推定するかを表にしておきましょう:

母数 推定量 関係する関数 オプション
母平均 標本平均 np.mean
母分散 不偏分散 np.var ddof=1
母標準偏差 不偏分散の正の平方根 np.std ddof=1
母相関係数 標本相関係数 np.corrcoef 対応する要素

課題3-1

前項で例に出した「身長データ」を用いて、母集団の分散と標準偏差を点推定した値を求めよ。 またそれに用いたコードをあわせて答えよ。

In [ ]:
# kadai 3-1 計算・解答用

標本抽出に伴う誤差

標本から推定値を得る手順よりも大事なこと:

  • 実際の母数の値にどの位近い推定値が得られるか?
  • 推定の結果はどのくらい信用できるか?

  • 標本誤差: 母数の値と推定値とのズレ--- データを沢山集めれば、生じうる誤差は小さくできる
  • 標本分布: どのような推定値が得られる可能性があるか、についての手がかり

推定値の信頼性を調べる方法

推定値を計算する元になる、標本に含まれる「個々のデータの値」の決定方法について見ていきましょう。

  1. 標本抽出の方法--- 単純無作為抽出と呼ばれる方法が前提
    • 単純無作為抽出: 母集団の中のどのデータにも平等に選ばれる可能性を持っているような標本抽出の方法
    • 無作為標本: 単純無作為抽出によって得られた標本
  2. 単純無作為抽出によって得られるデータの性質としての確率変数
    • 確率変数: 実際に結果が得られるまでどのような値が得られるかが決まっていない変数
      例: サイコロを振った時の出目
      同じ手続でデータを得ても結果に再現性がない⇒ 単純無作為抽出によって得られるデータは確率変数といえる
  3. 確率変数がどのような値をとるのかを示す確率分布
    • 確率変数の実現値:抽出の結果、確率変数がとる値
    • 確率分布: 確率変数がどのような値をどのような確率でとるかということを表した分布
    • 確率分布を用いた母集団の表現としての母集団分布
      • 母集団分布:ある変数の母集団における分布
        無作為抽出の場合、標本データの確率分布は母集団分布と同じになります
    • 代表的な母集団分布である正規分布
    • Pythonを使って正規分布の母集団から標本を抽出する方法

[参考: サイコロの出目の確率分布] 理論的には、サイコロの出目の確率分布はつぎのようになります:

サイコロの出る目 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

(2) 60,000回サイコロを振った場合:

サイコロの目 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
このように、非常に多くの実現値が得られれば、度数分布は確率分布とほぼ同じになります

In [30]:
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()

正規母集団から単純無作為抽出を行う

正規母集団から単純無作為抽出を行う場合、特にサンプルサイズが小さいと、 確率変数が正規分布に従っているかどうかは標本データのヒストグラムを見てもわからない(左図)。 しかし大きなサイズになれば、ほぼ正規分布に従うことわかる(右図).

In [103]:
import numpy.random as random
help(random.normal)
Help on built-in function normal:

normal(...)
    normal(loc=0.0, scale=1.0, size=None)
    
    Draw random samples from a normal (Gaussian) distribution.
    
    The probability density function of the normal distribution, first
    derived by De Moivre and 200 years later by both Gauss and Laplace
    independently [2]_, is often called the bell curve because of
    its characteristic shape (see the example below).
    
    The normal distributions occurs often in nature.  For example, it
    describes the commonly occurring distribution of samples influenced
    by a large number of tiny, random disturbances, each with its own
    unique distribution [2]_.
    
    Parameters
    ----------
    loc : float
        Mean ("centre") of the distribution.
    scale : float
        Standard deviation (spread or "width") of the distribution.
    size : int or tuple of ints, optional
        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
        ``m * n * k`` samples are drawn.  Default is None, in which case a
        single value is returned.
    
    See Also
    --------
    scipy.stats.distributions.norm : probability density function,
        distribution or cumulative density function, etc.
    
    Notes
    -----
    The probability density for the Gaussian distribution is
    
    .. math:: p(x) = \frac{1}{\sqrt{ 2 \pi \sigma^2 }}
                     e^{ - \frac{ (x - \mu)^2 } {2 \sigma^2} },
    
    where :math:`\mu` is the mean and :math:`\sigma` the standard
    deviation. The square of the standard deviation, :math:`\sigma^2`,
    is called the variance.
    
    The function has its peak at the mean, and its "spread" increases with
    the standard deviation (the function reaches 0.607 times its maximum at
    :math:`x + \sigma` and :math:`x - \sigma` [2]_).  This implies that
    `numpy.random.normal` is more likely to return samples lying close to
    the mean, rather than those far away.
    
    References
    ----------
    .. [1] Wikipedia, "Normal distribution",
           http://en.wikipedia.org/wiki/Normal_distribution
    .. [2] P. R. Peebles Jr., "Central Limit Theorem" in "Probability,
           Random Variables and Random Signal Principles", 4th ed., 2001,
           pp. 51, 51, 125.
    
    Examples
    --------
    Draw samples from the distribution:
    
    >>> mu, sigma = 0, 0.1 # mean and standard deviation
    >>> s = np.random.normal(mu, sigma, 1000)
    
    Verify the mean and the variance:
    
    >>> abs(mu - np.mean(s)) < 0.01
    True
    
    >>> abs(sigma - np.std(s, ddof=1)) < 0.01
    True
    
    Display the histogram of the samples, along with
    the probability density function:
    
    >>> import matplotlib.pyplot as plt
    >>> count, bins, ignored = plt.hist(s, 30, normed=True)
    >>> plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *
    ...                np.exp( - (bins - mu)**2 / (2 * sigma**2) ),
    ...          linewidth=2, color='r')
    >>> plt.show()

In [106]:
samplesSmall = random.normal(50, 10.0, 10)  # N(50,10^2)に従う乱数を10個生成
samplesSmall
Out[106]:
array([ 33.09788212,  58.53310324,  46.89549989,  73.32767525,
        59.72873003,  57.4554278 ,  54.04003881,  35.12864043,
        40.43082595,  56.46944824])
In [107]:
plt.hist(samplesSmall)
plt.show()
In [108]:
samplesLarge = random.normal(50, 10.0, 10000)  # N(50,10^2)に従う乱数を10000個生成
plt.hist(samplesLarge)
plt.show()

標本分布

標本分布とは、標本統計量(標本平均、標本分散など)に関する確率分布のこと

標本分布は標本における個々のデータの実現値を表した度数分布ではなく、標本統計量の確率分布
標本分布は母集団分布と標本統計量の種類、そしてサンプルサイズが決まると理論的(数学的)に導かれる。 実際のデータから作成されるわけではない

標本分布から何がわかるのか

標本分布を調べることで、推定値の性質が分かる可能性がある


標本分布を調べるときの観点:

  1. 標本分布が母数の本当の値を中心として分布しているか⇒平均
  2. 標本分布が横に大きく広がっていないか⇒標準偏差
  (1)が満たされないと推定値と母数との誤差が大きくなる
  (2) は推定値にどの程度誤差が生じるか調べるもの

標本分布を「経験的」に求める

Pythonを用いて「経験的」に標本分布を求める方法---「経験的」とは「現実に得られたデータに基づく」という意味、 理論的な標本分布を「近似的」に再現


母集団からサンプルサイズnの標本を何度も繰り返し抽出し、そのたびに標本統計量の実現値を計算して記録する


これを実行するときの問題

  • 母集団分布がどのような分布になるかが不明  
  • 母集団分布がわからないと、標本から得られる推定値の信頼性も不明
    ⇒ 母集団が正規分布であると仮定して「もし母集団分布がこのような正規分布ならば、このくらいあてになる推定値が得られる」ということを検討

正規母集団の母平均の推定

仮定:母集団分布は$N(50,10^2)$、 サンプルサイズ n =10

In [109]:
samples = random.normal(50, 10.0, 10)  # N(50,10^2)に従う乱数を10個生成

samples # データの確認
Out[109]:
array([ 48.19415007,  63.26636104,  57.49745353,  39.96272656,
        51.64073697,  60.36259109,  58.71608459,  63.91610309,
        65.66073452,  45.49907224])
In [110]:
np.mean(samples)   # 平均
Out[110]:
55.471601370394922

この結果、母平均の推定値は 55.5

[疑問: 55.4716... なのに答は 55.5?] 注意: 55.4716... と表示されたとしても、答として、これをこのまま答えないこと。

有効数字の桁数を考慮して、55.5 と答える(有効数字は2桁、計算途中なので3桁とする)。

母平均は50なので、推定値との差は 5.5 --- この差が標本誤差(誤差)

一回きりの標本抽出で推定値がどれくらいの誤差をもつかは、いろいろ

標本分布を求める

先の操作を10,000回繰り返してみよう:

In [3]:
import numpy as np
import numpy.random as random
meanArray = np.array([ np.mean(random.normal(50, 10.0, 10)) for _ in range(10000)])

こうして得られた標本平均のヒストグラムと母集団の正規分布とを、次のコードにより重ねあわせてみると、かなり近い関係にあることが分かります:

In [4]:
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以下の範囲で得られた回数を調べてみましょう:

In [113]:
# 45以上55以下のものを1それ以外を0としたリストを作り、そこで0でない(つまり1)の個数を数える
np.count_nonzero([1 if 45 <= x <= 55 else 0 for x in meanArray]) 
Out[113]:
8868

全体の89%は本当の母数±5に収まっていることがわかります。理論的にはこのデータは$N(50,10^2)$に従うので、次章の知識を先取りしscipy.statsモジュールの関数cdfを用いた理論的な値と比較してみましょう:

In [60]:
from scipy.stats import norm
help(norm.cdf)
Help on method cdf in module scipy.stats._distn_infrastructure:

cdf(self, x, *args, **kwds) method of scipy.stats._continuous_distns.norm_gen instance
    Cumulative distribution function of the given RV.
    
    Parameters
    ----------
    x : array_like
        quantiles
    arg1, arg2, arg3,... : array_like
        The shape parameter(s) for the distribution (see docstring of the
        instance object for more information)
    loc : array_like, optional
        location parameter (default=0)
    scale : array_like, optional
        scale parameter (default=1)
    
    Returns
    -------
    cdf : ndarray
        Cumulative distribution function evaluated at `x`

In [63]:
norm.cdf(55, loc=50, scale=10.0**0.5) - norm.cdf(45,loc=50,scale=10.0**0.5)
Out[63]:
0.88615370199334187

このように、 上位二桁まで一致していることが確かめられました。

平均$\mu$、分散$\sigma^2$の正規分布 $N(\mu,\sigma^2)$ に従う母集団からサンプルサイズ$n$の標本を抽出したとき、標本平均の標本分布は、理論的に$N(\mu,\sigma^2/n)$に従うことは、次のようにして実験的に確かめることができます。

In [114]:
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))
サンプルサイズ=50, 平均=50.001157, 標本分散=2.018868, 標本分散*サンプルサイズ=100.943413
サンプルサイズ=60, 平均=49.989521, 標本分散=1.637274, 標本分散*サンプルサイズ=98.236410
サンプルサイズ=70, 平均=50.011332, 標本分散=1.432251, 標本分散*サンプルサイズ=100.257560
サンプルサイズ=80, 平均=50.002751, 標本分散=1.255702, 標本分散*サンプルサイズ=100.456120
サンプルサイズ=90, 平均=50.011754, 標本分散=1.103067, 標本分散*サンプルサイズ=99.276044
サンプルサイズ=100, 平均=49.991631, 標本分散=1.002772, 標本分散*サンプルサイズ=100.277168
サンプルサイズ=150, 平均=49.996742, 標本分散=0.664137, 標本分散*サンプルサイズ=99.620536
サンプルサイズ=200, 平均=50.005825, 標本分散=0.500871, 標本分散*サンプルサイズ=100.174264
サンプルサイズ=250, 平均=50.000568, 標本分散=0.405689, 標本分散*サンプルサイズ=101.422301
サンプルサイズ=300, 平均=49.997763, 標本分散=0.329371, 標本分散*サンプルサイズ=98.811296

以上の結果から、サンプルサイズを変えても標本平均の平均が母平均にほぼ一致すること、また標本平均の分散は、母分散/サンプルサイズであることが確かめられました。

不偏性

ある推定量の標本分布の平均が、推定しようとしている母数の値と一致するとき、その推定量は不偏性があるといいます。また、 不偏性がある推定量のことを不偏推定量と言います。 推定量の不偏性は「標本分布が母数の本当の値を中心として分布しているか」に対応した概念です。

標準誤差

「推定量の標本分布の標準偏差(へんさ)」を標準誤差(ごさ)と言います。 これは「標本分布の横への広がり」を表す値です。標準誤差が小さいということは、 標本分布の広がりが小さいことを表し、これは『ばらつきが少ない』ことを意味します。

例: $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が大きい時」はほぼ正規分布になる

課題3-2

標本分布を求めるの項目ではn=10の標本を10,000個作り、その標本平均を求めた。 この課題ではn=100として10,000個の標本平均を求めてみよう。そしてn=10の場合のヒストグラムとn=100の場合のヒストグラムを比較し、標準誤差が小さくなっていることを確かめてみよ。

In [115]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy.random as random
# 計算・解答用
std of sample10=3.137151, sample100=0.990142

標本平均以外の標本分布

「平均」以外の母数のうち、分散と中央値の標本分布について述べます。

標本分散と不偏分散の標本分布</H4>

母分散の推定量としての標本分散と不偏分散の違いを実験で示しましょう。

In [117]:
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())
90.2677729196
100.297525466

平均的に母分散の推定値として不偏分散のほうが近い値であることがわかります。次にその標準偏差を求めてみましょう。

In [118]:
print(SampleVars.std())
print(UnbiasedVars.std())
42.4670221473
47.1855801637

不偏分散のほうが推定値のばらつきが大きいことがわかります。

In [119]:
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()
In [120]:
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()
  • 不偏分散は不偏性をもつという点では望ましい性質をもっているが、 標本分散に比べて過大な推定値が得られる可能性が高い
  • 不偏分散の平方根は母標準偏差の不偏推定量ではない

中央値の標本分布</H4>

実験: 標本中央値の標本分布を求め、標本平均と比較する(どちらが母平均の推定量としてよいか?)

In [121]:
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())
50.0321506814
50.0576939024

標本中央値も標本平均も母平均とほぼ一致---母平均の推定量としての候補として優劣が付けられない

In [122]:
print(np.std(SampleAverage))
print(np.std(SampleMedian))
3.15787180632
3.6923757235

標本平均のほうが標準誤差が小さい---標本平均は母平均の不偏推定量の中で</font>標準誤差が最小

In [123]:
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()

関数のまとめ

注: numpyをnp, np.randomを random、 matplotlib.pyplotをplt、pandasをpd、scipy.statsをstと略記する

目的 関数名とモジュール 使い方
一様分布に従う乱数(復習) 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})$の正規分布の関数描画


演習問題3

演習問題3-1

$N(50,10^2)$の正規母集団から$n = 20$の標本抽出を5,000回繰り返すことにより、標本平均の経験的な標本分布を求めよ。またこの経験的な標本分布から得られるグラフ(ヒストグラム)に、理論的に導かれる標本分布を重ねあわせて表示せよ。

この結果下のようなグラフが得られればよい:

In [ ]:
# 計算・解答用

演習問題3-2

理論的な標本分布について、サンプルサイズを$n=1, 4, 9, 16, 25$と変化させた時に、標本分布の形状がどのように代わるかを調べてみよ。ただし、母集団分布は標準正規分布$N(0, 1^2)$とする。


これによって次のような図が得られれば良い。

In [ ]:
# 計算・解答用

演習問題3-3

最小値0、最大値100の一様分布から$n=20$の標本抽出を5,000回繰り返すことにより、標本平均の経験的な標本分布を求めよ。中心極限定理からこれは正規分布に近づくはずである。そこでこのヒストグラムに、標本平均データの平均と標準偏差による正規分布のグラフを重ねせて表示せよ。
この結果、次のようなグラフが得られればよい:

In [ ]:
# 計算・解答用