母集団と標本の関係
まず重要な用語をまとめておきましょう:
- 母集団: 対象のデータ全体
- 標本: 母集団から一部のデータを取り出した(抽出した)もの
- 標本抽出: 母集団から標本をとりだすこと
- 母数: 母集団の性質をあらわす統計的指数（比率、平均、分散、相関係数など） </UL> 例えば、日本の中学生全体のテスト得点データを母集団とすると、名古屋市昭和区の(一部の)中学生のテスト得点データが「標本」になります。そして、平均点や分散などが母数とみなせます。
推測統計
手元にある標本から母集団の性質(母数)を推測することを推測統計といいます。 推測統計は推定と検定に分けられます。
「推定」は、母数の値に対し具体的な値を求める点推定、つまり一つの値で推定の結果を表すものと、区間推定、つまりある幅を持った区間で結果を表すものに分けられます。
また「検定」は母集団について述べた異なる立場の二つの主張(仮説)のうちどれを採用するかを決定するものです。例えば、日本の中学生のテスト得点が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 は、平均値を求めるメソッド(関数)でした。
print(Height.mean())
推定の基礎
まず用語を確認しましょう:
- 標本統計量：標本データから計算される平均や分散などの値
- 推定値：標本データを用いて計算した母数の推定量
</A>
</CENTER>
ここで、標本統計量からどのように母数を推定するかを表にしておきましょう：
| 母数 | 推定量 | 関係する関数 | オプション |
|---|---|---|---|
| 母平均 | 標本平均 | np.mean | |
| 母分散 | 不偏分散 | np.var | ddof=1 |
| 母標準偏差 | 不偏分散の正の平方根 | np.std | ddof=1 |
| 母相関係数 | 標本相関係数 | np.corrcoef | 対応する要素 |
練習問題 前項で例に出した「身長データ」を用いて、母集団の分散と標準偏差を点推定した値を求めよ。 またそれに用いたコードをあわせて答えよ。
# Exercise1 計算・解答用
# 母分散と母集団の標準偏差の点推定
Height.var(ddof=1), Height.std(ddof=1)
# 比較 (標本分散, 標準偏差)
# こちらのほうが小さい値になっている
Height.var(), Height.std()
標本抽出に伴う誤差
標本から推定値を得る手順よりも大事なこと:
- 実際の母数の値にどの位近い推定値が得られるか？
- 推定の結果はどのくらい信用できるか？
- 標本誤差： 母数の値と推定値とのズレ--- データを沢山集めれば、生じうる誤差は小さくできる
- 標本分布: どのような推定値が得られる可能性があるか、についての手がかり
推定値の信頼性を調べる方法
推定値を計算する元になる、標本に含まれる「個々のデータの値」の決定方法について見ていきましょう。
- 標本抽出の方法---
単純無作為抽出と呼ばれる方法が前提
- 単純無作為抽出: 母集団の中のどのデータにも平等に選ばれる可能性を持っているような標本抽出の方法
- 無作為標本: 単純無作為抽出によって得られた標本
- 単純無作為抽出によって得られるデータの性質としての確率変数
- 確率変数: 実際に結果が得られるまでどのような値が得られるかが決まっていない変数
例： サイコロを振った時の出目
同じ手続でデータを得ても結果に再現性がない⇒ 単純無作為抽出によって得られるデータは確率変数といえる
- 確率変数: 実際に結果が得られるまでどのような値が得られるかが決まっていない変数
- 確率変数がどのような値をとるのかを示す確率分布
- 確率変数の実現値：抽出の結果、確率変数がとる値
- 確率分布: 確率変数がどのような値をどのような確率でとるかということを表した分布
- 確率分布を用いた母集団の表現としての母集団分布
- 母集団分布：ある変数の母集団における分布
無作為抽出の場合、標本データの確率分布は母集団分布と同じになります
- 母集団分布：ある変数の母集団における分布
- 代表的な母集団分布である正規分布
- 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 |
正規分布
正規分布（normal distribution, ガウス分布 Gaussian distributionともいう)は統計学で最もよく用いられる確率分布なので、ここで復習しておこう
正規分布の表現方法： N (μ,σ2) --- μは平均値、σ2は分散 (σは標準偏差)
</A>
</CENTER>
|
|
| |||||||||||||||||||||
| 正規分布 N (0,1)のグラフ | 正規分布 N (1,1) のグラフを重ね書き | </TR> </TABLE>
</A></TD>
</A> </TD>
</TR>|
|
|
|
| μ-σからμ+σ の範囲 | μ-2σからμ+2σ の範囲 | μ-3σからμ+3σの範囲 |
</A></TD>
</A> </TD>
</A> </TD></TR>
| 目的 | 関数名とモジュール | 使い方 |
|---|---|---|
| 一様分布に従う乱数(復習) | 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})$の正規分布の関数描画 |
演習問題
演習問題1
$N(50,10^2)$の正規母集団から$n = 20$の標本抽出を5,000回繰り返すことにより、標本平均の経験的な標本分布を求めよ。またこの経験的な標本分布から得られるグラフ(ヒストグラム)に、理論的に導かれる標本分布を重ねあわせて表示せよ。
この結果下のようなグラフが得られればよい:
# 計算・解答用
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as random
import scipy.stats as st
data = np.array([np.mean([random.normal(50,10,20)]) for i in range(5000)])
plt.hist(data,color='y',density=True,bins=20, ec='k')
x = np.linspace(40, 60, 1000) # x = np.arange(35,65,0.01)でもよい
plt.plot(x, st.norm.pdf(x,50,(100/20)**0.5), ls='--',label=r'$N(50,\frac{10^2}{20})$')
plt.xlabel('mean of samples')
plt.ylabel('Density')
plt.title('Histogram of means of samples')
plt.legend()
plt.plot()
演習問題2
理論的な標本分布について、サンプルサイズを$n=1, 4, 9, 16, 25$と変化させた時に、標本分布の形状がどのように代わるかを調べてみよ。ただし、母集団分布は標準正規分布$N(0, 1^2)$とする。
これによって次のような図が得られれば良い。
![](http://lang.sist.chukyo-u.ac.jp/Classes/DataScience/Figs/RS-3-B.png")
# 計算・解答用
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as st
x = np.linspace(-3,3,1000)
for n in range(1,6):
plt.plot(x,st.norm.pdf(x,0,(1/n)),label=r"$n=%i^2$"%(n))
plt.xlabel('x')
plt.ylabel('norm(x,0,np.sqrt(1/n))')
plt.legend()
plt.show()
?plt.plot
演習問題3
最小値0、最大値100の一様分布から$n=20$の標本抽出を5,000回繰り返すことにより、標本平均の経験的な標本分布を求めよ。中心極限定理からこれは正規分布に近づくはずである。そこでこのヒストグラムに、標本平均データの平均と標準偏差による正規分布のグラフを重ねせて表示せよ。
この結果、次のようなグラフが得られればよい:
# 計算・解答用
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as random
import scipy.stats as st
data = np.array([np.mean(random.random(20)*100) for i in range(5000)])
plt.hist(data,color='y',density=True,bins=20,ec='k')
x = np.linspace(25, 75, 1000) # x = np.arange(35,65,0.01)でもよい
plt.plot(x, st.norm.pdf(x,50,(10000/(12*20))**0.5), ls='--',label=r'$N(50,\frac{100^2}{12*20})$')
plt.ylabel('Density')
plt.xlabel('Mean of samples')
plt.title('Histogram of means of samples')
plt.legend()
plt.plot()
区間推定¶
この項は次の書籍を参考にしている: 馬場孝之(2014)『統計学キャンパスゼミ』マセマ
- 予備知識:
- 信頼係数とは、母数の推定において、ある区間のなかに母数が入っている確率が1に十分近い値$1-\alpha$になるようにしたときの値。またこの$\alpha$は(本来は『検定」の用語である)有意水準に対応する
- 逆に$\alpha$を有意水準とすると、$1-\alpha$は信頼係数。$\alpha$として1%や5%を用いることが慣習的
- 区間推定とは、推定する母数$\theta$に対し、$P(\theta_1 \leq \theta \leq \theta_2) = 1-\alpha$を満たす信頼区間$\theta_1$と$\theta_2$の値を求めるもの
以下は、正規分布$N(\mu, \sigma^2)$に従う母集団から無作為抽出した標本$X_1, X_2, \ldots, X_n$に対する推定についてのみ述べる。$\bar{X}=\frac{1}{n}\sum_{i=1}^n X_i$とする。
- 中心極限定理（復習）
- 互いに独立な$n$個の確率変数$X_1, X_2, \ldots,X_n$が平均$\mu$, 分散$\sigma^2$の同一の確率分布に従うとき、 $\bar{X}=\frac{1}{n} (X_1+X_2+\cdots+X_n)$(つまり$\bar{X}$はその平均)は$n\rightarrow \infty$のとき 正規分布$N(\mu, \sigma^2/n)$に従う。
- そして$\bar{X}$を標準化した確率変数$Z=\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}$は、$n \rightarrow \infty$のとき標準正規分布$N(0,1)$に従う
- 正規分布$N(\mu, \sigma^2)$：復習
- 平均値$\mu=0$、分散$\sigma^2=1$の正規分布を標準正規分布と呼ぶ。これは、正規分布に従う確率変数$X$を標準化: $Z=\frac{X-\mu}{\sigma}$ した変数が従う分布である。
- 互いに独立な$n$個の確率変数$X_1, X_2, \ldots,X_n$がすべて同一の正規分布$N(\mu,\sigma^2)$に従うとき、 $\bar{X}=\frac{1}{n} (X_1+X_2+\cdots+X_n)$(つまり$\bar{X}$はその平均)は正規分布$N(\mu, \sigma^2/n)$に従う。 (注：中心極限定理と違い、$n \rightarrow \infty$が不要であることに注意）
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# Gaussian : N(0,1)
def Gaussian (x):
return (np.exp(-x**2/2)) / np.sqrt(2*np.pi)
# [-6.0, 6.0] の区間で等差数列を作成
x = np.linspace(-6.0, 6.0, 10000)
# f(x)の結果を得る: [Gaussian(-6.0),...,Gaussian(6.0)]
y = [Gaussian(i) for i in x]
# グラフに表示
plt.plot(x,y)
plt.axis([-6,6,-0.02,0.45]) # 描画範囲を指定する x:[-6.0,6.0], y:[-0.02, 0.45]
plt.title("standard normal distribution: N(0,1)")
plt.grid()
- ここで取り上げる区間推定は、母集団が正規分布$N(\mu,\sigma^2)$に従うもので、そこから無作為に抽出した標本$X_1, X_2, \ldots, X_n$の母数を区間推定するもの。
- 推定する母数が平均値か分散か、また$\sigma^2$が既知か未知か、で条件が異なる
$\sigma^2$が既知のときの母平均$\mu$の推定¶
- 母集団が正規分布$N(\mu,\sigma^2)$に従い($\mu$は未知、$\sigma^2$は既知)、そこから無作為に抽出した標本$X_1, X_2, \ldots, X_n$から$\mu$を推定する
- $\bar{X}=\frac{1}{n}\sum_{i=1}^n X_i$とするーーー正規分布の性質から$\bar{X}$は$N(\mu, \frac{\sigma^2}{n})$に従う
- 新たな確率変数$Z=\frac{\bar{X}-\mu}{\sqrt{\sigma^2/n}}$を定義すると$Z$は$N(0,1)$に従う($Z \sim N(0,1)$とも書ける)
- これを用いて、標準正規分布のpercent point関数(Pythonではscipy.stats.norm.ppf関数)ppf(x)を用いると、 $$P(-{\rm ppf}(\frac{\alpha}{2}) \leq Z \leq {\rm ppf}(\frac{\alpha}{2})) = 1-\alpha$$
であるから信頼区間$(1-\alpha)$を満たすには $$-{\rm ppf}(\frac{\alpha}{2}) \leq \frac{\bar{X}-\mu}{\sqrt{\sigma^2/n}} \leq {\rm ppf}(\frac{\alpha}{2})$$ これから、$\bar{X} - {\rm ppf}(\frac{\alpha}{2})\sqrt{\frac{\sigma^2}{n}} \leq \mu \leq \bar{X} +{\rm ppf}(\frac{\alpha}{2})\sqrt{\frac{\sigma^2}{n}}$
- 参考: ppf(0.025)=-1.96, ppf(1-0.025)=1.96
from scipy import stats as st
import matplotlib.pyplot as plt
import numpy as np
lower=st.norm.ppf(0.025)
upper=st.norm.ppf(1-0.025)
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(lower,c='b',ls='--')
plt.axvline(upper,c='b',ls='--')
x = np.linspace(lower, upper, 10000)
plt.fill_between(x, st.norm.pdf(x,0,1),color="green")
plt.legend()
plt.title("95% confidence interval")
# 標準正規分布のppf関数の0.025(2.5%)と0.975(97.5%)点
import scipy.stats as st
print("ppf(0.025)={:5.3f}, ppf(0.975)={:5.3f}".format(st.norm.ppf(0.025),st.norm.ppf(1-0.025)))
練習問題2-1
$N(\mu, 16)$に従う母集団から9個の標本$X_1,X_2,\cdots,X_9$を無作為抽出したところ、その標本平均は$\bar{X}=10.0$であった。これから、母平均$\mu$の95%信頼区間を求める。
from scipy import stats as st
import matplotlib.pyplot as plt
import numpy as np
# 有意水準
alpha=0.05
# 標本平均のサイズ
n=9
# 標本平均
xBar=10.0
# 母分散
var=16
# ppf(alpha/2) * (var/n)**0.5
print(st.norm.ppf(0.025)*(var/n)**0.5)
print(st.norm.ppf(0.975)*(var/n)**0.5)
# 以上から
print("{:6.3f}<=μ<={:6.3f}".format(xBar+st.norm.ppf(0.025)*(var/n)**0.5, xBar+st.norm.ppf(0.975)*(var/n)**0.5))
問題: 上と同じ状況において、母平均μの99%信頼区間を求めよ。
# 計算
print("{:6.3f}<=μ<={:6.3f}".format(xBar+st.norm.ppf(0.005)*(var/n)**0.5, xBar+st.norm.ppf(0.995)*(var/n)**0.5))
練習問題2-2
$N(\mu, 25)$に従う母集団から16個の標本$X_1,X_2,\cdots,X_{16}$を無作為抽出したところ、その標本平均は$\bar{X}=63.0$であった。これから、母平均$\mu$の95%信頼区間を求めよ。
# 計算
# 有意水準
alpha=0.05
# 標本平均のサイズ
n=16
# 標本平均
xBar=63.0
# 母分散
var=25
# 以上から
print("{:6.3f}<=μ<={:6.3f}".format(xBar+st.norm.ppf(0.025)*(var/n)**0.5, xBar+st.norm.ppf(0.975)*(var/n)**0.5))
確率分布の復習¶
$\chi^2$分布
互いに独立な$n$個の確率変数$X_1, X_2, \ldots,X_n$が標準正規分布$N(0,1)$に従うとき、その2乗和$Z=X_1^2+\cdots+X_n^2$は、自由度$n$の$\chi^2$分布に従う
$t$分布
2つの独立な確率変数$X$と$Y$があり、$X$は標準正規分布$N(0,1)$に、$Y$は自由度$n$の$\chi^2$分布に従うとき、$Z=\frac{X}{\sqrt{Y/n}}$とすると、$Z$は自由度$n$の$t$分布に従う
$\sigma^2$が未知のときの母平均$\mu$の推定¶
- 母集団が正規分布$N(\mu,\sigma^2)$に従い($\mu$も$\sigma^2$も未知)、そこから無作為に抽出した標本$X_1, X_2, \ldots, X_n$から$\mu$を推定する
- $\bar{X}=\frac{1}{n}\sum_{i=1}^n X_i$とするーーー正規分布の性質から$\bar{X}$は$N(\mu, \frac{\sigma^2}{n})$に従うが、$\sigma^2$は未知！
- 標本の不偏分散$S^2$を母分散$\sigma^2$の代わりに用いる。確率変数$Z=\frac{\bar{X}-\mu}{\sqrt{S^2/n}}$を定義すると $Z$は 自由度$n-1$のt分布に従う
その根拠:
- $Z=\frac{ \frac{\bar{X}-\mu}{\sqrt{\sigma^2/n}}}{\sqrt{S^2/\sigma^2}}$と書き換える
- ここで、$X=\frac{\bar{X}-\mu}{\sqrt{\sigma^2/n}}$とおくと$Z=\frac{X}{\sqrt{S^2/\sigma^2}}$と書け、$X$は$N(0,1)$に従う
- $Y = (n-1)\cdot\frac{S^2}{\sigma^2}$ とおくと(つまり$\frac{S^2}{\sigma^2}=\frac{Y}{n-1}$）、 $Z=\frac{X}{\sqrt{Y/(n-1)}}$の形になる
- $S^2$は不偏分散であるから $$Y=\frac{n-1}{\sigma^2}S^2 = \frac{n-1}{\sigma^2}\cdot \frac{1}{n-1}\sum_{i=1}^n (X_i-\bar{X})^2 = \sum_{i=1}^n \left( \frac{X_i-\bar{X}}{\sigma}\right)^2$$
- $\frac{X_i - \bar{X}}{\sigma}$は$X_i$の標準化と呼ばれるもの、これは$N(0,1)$に従う
- したがって$Y$はこの２乗和であるから、$Y$は($\frac{X_i-\bar{X}}{\sigma}$は線形従属なので自由度が$n$より1減って) 自由度$n-1$の$\chi^2$分布に従う
- 分散が既知の場合と同様、信頼区間$1-\alpha$とすると、 $t$分布のpercent point関数(Pythonではscipy.stats.t.ppf(x,自由度)、以下ではt.ppfと表す)を用いて $$\bar{X} - {\rm t.ppf}(\frac{\alpha}{2},n-1)\sqrt{\frac{S^2}{n}} \leq \mu \leq \bar{X} +{\rm t.ppf}(\frac{\alpha}{2},n-1)\sqrt{\frac{S^2}{n}}$$
# 自由度9のときの2.5%点と97.5%点
import scipy.stats as st
print("t.ppf(0.025,3)={:5.3f}, t.ppf(0.975,3)={:5.3f}".format(st.t.ppf(0.025,3),st.t.ppf(1-0.025,3)))
print("t.ppf(0.025,9)={:5.3f}, t.ppf(0.975,9)={:5.3f}".format(st.t.ppf(0.025,9),st.t.ppf(1-0.025,9)))
print("t.ppf(0.025,100)={:5.3f}, t.ppf(0.975,100)={:5.3f}".format(st.t.ppf(0.025,100),st.t.ppf(1-0.025,100)))
# 作図
%matplotlib inline
from scipy import stats as st
import numpy as np
import matplotlib.pyplot as plt
fig,ax = plt.subplots(1,1)
x = np.linspace(-6, 6, 10000)
linestyles = [':', '--', '-.','-']
deg_of_freedom = [1,3,9,100]
for df, ls in zip(deg_of_freedom, linestyles):
ax.plot(x, st.t.pdf(x, df), linestyle=ls, label=r'$df=%i$' % df)
plt.xlabel('$t$')
plt.ylabel('$f(t)$')
plt.title('t-Distribution')
plt.axis([-6,6,-0.02,0.6])
plt.legend()
plt.grid()
練習問題3-1
$N(\mu, \sigma^2)$に従う母集団から10個の標本を無作為抽出したところ、その値は9.8, 10.3, 11.9, 8.7, 9.5, 10.2, 10.9,8.2, 9.4, 10.5であった。これから、母平均$\mu$の95%信頼区間を求める。
# 自由度7のt分布のppf関数の0.025と0.975に対する値
import scipy.stats as st
import numpy as np
# 標本
Xs = np.array([9.8, 10.3, 11.9, 8.7, 9.5, 10.2, 10.9,8.2, 9.4, 10.5])
N=Xs.size
# 平均
xBar=Xs.mean()
N,xBar
# 自由度
ddof=N-1
# 不偏分散
var=Xs.var(ddof=1)
var
# t.ppf(alpha/2,ddof) * (var/n)**0.5
print(st.t.ppf(0.025,ddof)*(var/N)**0.5)
print(st.t.ppf(0.975,ddof)*(var/N)**0.5)
# 以上から
print("{:6.3f}<=μ<={:6.3f}".format(xBar+st.t.ppf(0.025,ddof)*(var/N)**0.5, xBar+st.t.ppf(0.975,ddof)*(var/N)**0.5))
%matplotlib inline
from scipy import stats as st
import numpy as np
import matplotlib.pyplot as plt
lower=st.t.ppf(0.025,ddof)
upper=st.t.ppf(0.975,ddof)
df=ddof
x = np.linspace(-6, 6, 10000)
plt.plot(x, st.t.pdf(x, df))
plt.xlabel('$t$')
plt.ylabel('$f(t)$')
plt.title('t-Distribution')
plt.axis([-6,6,-0.02,0.6])
plt.axvline(lower,c='yellow',ls='--')
plt.axvline(upper,c='orange',ls='--')
plt.axhline(0,c='k',ls='-')
x = np.linspace(lower, upper, 10000)
plt.fill_between(x, st.t.pdf(x,df),color="green")
plt.grid()
問題: 上と同じ状況において、母平均μの99%信頼区間を求めよ。
# 標本
Xs = np.array([9.8, 10.3, 11.9, 8.7, 9.5, 10.2, 10.9,8.2, 9.4, 10.5])
N=Xs.size
# 平均
xBar=Xs.mean()
# 自由度
ddof=N-1
# 不偏分散
var=Xs.var(ddof=1)
print("{:6.3f}<=μ<={:6.3f}".format(xBar+st.t.ppf(0.005,ddof)*(var/N)**0.5, xBar+st.t.ppf(0.995,ddof)*(var/N)**0.5))
練習問題3-2
ある一年間の新生児から８名の体重を無作為抽出した結果が以下である:(単位はkg)
2.82, 3.01, 3.24, 3.16, 3.34, 3.57, 2.98, 3.28
全国の新生児の体重を母集団とし、それが正規分布に従うと仮定したとき、母平均の95%信頼区間を推定せよ。
Xs=np.array([2.82, 3.01, 3.24, 3.16, 3.34, 3.57, 2.98, 3.28])
N=Xs.size
# 平均
xBar=Xs.mean()
# 自由度
ddof=N-1
# 不偏分散
var=Xs.var(ddof=1)
print("Xbar = {:6.3f}".format(xBar))
print("{:6.3f}<=μ<={:6.3f}".format(xBar+st.t.ppf(0.025,ddof)*(var/N)**0.5, xBar+st.t.ppf(0.975,ddof)*(var/N)**0.5))
母分散$\sigma^2$の推定($\mu$も未知)¶
- 母集団が正規分布$N(\mu,\sigma^2)$に従い($\mu$も$\sigma^2$も未知)、そこから無作為に抽出した標本$X_1, X_2, \ldots, X_n$から$\sigma^2$を推定する
- $V=\sum_{i=1}^n (\frac{X_i-\bar{X}}{\sigma})^2$とすると(前述から)$V$は自由度$n-1$の$\chi^2$分布に従う。このことから$\sigma^2$の区間推定を行う
- $\chi^2$分布のpercent point関数(Pythonではscipy.stats.chi2.ppf(x,自由度), 以下ではchi2.ppfと表す)を用いて $$ P(chi2.{\rm ppf}(1-\frac{\alpha}{2},n-1) \leq V \leq chi2.{\rm ppf}(\frac{\alpha}{2},n-1)) = 1-\alpha$$より、 $$chi2.{\rm ppf}(1-\frac{\alpha}{2},n-1) \leq \sum_{i=1}^n (\frac{X_i-\bar{X}}{\sigma})^2 \leq chi2.{\rm ppf}(\frac{\alpha}{2},n-1)$$
- ここで標本からの不偏分散$S^2 = \frac{1}{n-1}\sum_{i=0}^n (X_i - \bar{X}^2)$を用いると $$V=\sum_{i=1}^n (\frac{X_i-\bar{X}}{\sigma})^2=\frac{n-1}{\sigma^2} \cdot \frac{1}{n-1}\sum_{i=0}^n (X_i-\bar{X})^2 = \frac{n-1}{\sigma^2} S^2$$
- 以上から、$\frac{(n-1)S^2}{chi2.{\rm ppf}(\frac{\alpha}{2},n-1)} \leq \sigma^2 \leq \frac{(n-1)S^2}{chi2.{\rm ppf}(1-\frac{\alpha}{2},n-1)}$
練習問題4-1
$𝑁(\mu,\sigma^2)$に従う母集団から10個の標本を無作為抽出したところ、その値は9.8, 10.3, 11.9, 8.7, 9.5, 10.2, 10.9,8.2, 9.4, 10.5であった。これから、母分散$\sigma^2$の95%信頼区間を求める。
# 自由度7のt分布のppf関数の0.025と0.975に対する値
import scipy.stats as st
import numpy as np
# 標本
Xs = np.array([9.8, 10.3, 11.9, 8.7, 9.5, 10.2, 10.9,8.2, 9.4, 10.5])
N=Xs.size
# 平均
xBar=Xs.mean()
N,xBar
# 自由度
ddof=N-1
# 不偏分散
var=Xs.var(ddof=1)
var
# chi2.ppf(alpha/2,ddof)
print(st.chi2.ppf(0.025,ddof))
print(st.chi2.ppf(0.975,ddof))
# 以上から
print(" σ**2={:6.3f}".format(var))
print("{:6.3f}<=σ**2<={:6.3f}".format(var*ddof/st.chi2.ppf(0.975,ddof), var*ddof/st.chi2.ppf(0.025,ddof)))
%matplotlib inline
from scipy import stats as st
import numpy as np
import matplotlib.pyplot as plt
df=ddof
upper=st.chi2.ppf(0.025,ddof)
lower=st.chi2.ppf(0.975,ddof)
x = np.linspace(0, 20, 10000)
plt.plot(x, st.chi2.pdf(x, df))
plt.xlabel('$\chi^2$')
plt.ylabel(r'$f(\chi^2)$')
plt.title(r'$\chi^2\ \mathrm{Distribution}$')
plt.axis([-0.2,20,-0.02,0.2])
plt.axvline(lower,c='yellow',ls='--')
plt.axvline(upper,c='orange',ls='--')
plt.axhline(0,c='k',ls='-')
x = np.linspace(lower, upper, 10000)
plt.fill_between(x, st.chi2.pdf(x,df),color="green")
plt.grid()
問題: 上と同じ状況において、母分散の99%信頼区間を求めよ。
print(" σ**2={:6.3f}".format(var))
print("{:6.3f}<=σ**2<={:6.3f}".format(var*ddof/st.chi2.ppf(0.995,ddof), var*ddof/st.chi2.ppf(0.005,ddof)))
練習問題4-2
ある一年間の新生児から８名の体重を無作為抽出した結果が以下である:(単位はkg)
2.82, 3.01, 3.24, 3.16, 3.34, 3.57, 2.98, 3.28
全国の新生児の体重を母集団とし、それが正規分布に従うと仮定したとき、母分散の95%信頼区間を推定せよ。
Xs=np.array([2.82, 3.01, 3.24, 3.16, 3.34, 3.57, 2.98, 3.28])
N=Xs.size
# 平均
xBar=Xs.mean()
# 自由度
ddof=N-1
# 不偏分散
var=Xs.var(ddof=1)
# 以上から
print(" xBar={:6.3f}".format(xBar))
print("{:6.3f}<=μ<={:6.3f}".format(xBar+st.t.ppf(0.025,ddof)*(var/N)**0.5, xBar+st.t.ppf(0.975,ddof)*(var/N)**0.5))
print(" σ**2={:6.3f}".format(var))
print("{:6.3f}<=σ**2<={:6.3f}".format(var*ddof/st.chi2.ppf(0.975,ddof), var*ddof/st.chi2.ppf(0.025,ddof)))