教育データ分析入門2（第4回）

• 4 統計的推定2
  • 4.1 区間推定（分散が既知で標本が複数）
  • 4.2 標準化
  • 4.3 区間推定（分散が未知の場合）
  • 4.4 中心極限定理
  • 4.5 区間推定（母集団分布が未知の場合）
  • 練習問題1
  • 練習問題2
  • 練習問題3

import math
import numpy as np
from scipy.stats import t
import matplotlib.pyplot as plt
%matplotlib inline

4 統計的推定2

4.1 区間推定（分散が既知で標本が複数）

• 母集団から複数の標本を抽出して、「X（統計量）は○○%の確率でY1からY2の範囲に含まれる」と主張したい。

  • 例: 「東京都新宿区にある高校に通う生徒のあるテストの平均点は95%の確率で82点から87点の間に含まれる」

問2: 正規分布に従うある母集団($\mu$が未知、$\sigma^2=25$)から、無作為にx人の点数を取り出した。この時、信頼度95%で母集団の平均を推定せよ。

• $X$の分布と$\bar{X}$の分布を以下に示す。
• $\bar{X}$の分布は$X$の分布と比べるとより平均に寄っている。
• $\bar{X}$は平均$\mu$、分散$\frac{\sigma^2}{n}$の正規分布に従うことが知られている。

実際に$\bar{X}$と$X$の分布を確認してみよう

# 平均が100で標準偏差が10の正規分布から100000回のサンプリング
# これを母集団とみなす。
x = np.random.normal(100, 10, 100000)

# 母集団から100000回１００個のデータを抽出して平均を算出
M = []

for i in range(100000):
    L = np.random.choice(x,100)
    M.append(np.average(L))

• $X$と$\bar{X}$の分布

plt.hist(x,bins=100,alpha=0.5)
plt.hist(M,bins=100,alpha=0.5)
plt.show()

# Mの分散
np.var(M)

0.9972307856494186

# Mの分散の理論値
# 母集団の分散が100、nが100
100/100

1.0

• $\bar{X}$は平均$\mu$、分散$\frac{\sigma^2}{n}$の正規分布に従うので、９５%の確率で以下の不等式が成り立つ。

$$ \mu - 1.96 \times \frac{\sigma}{\sqrt{n}} \leq \bar{X} \leq \mu + 1.96 \times \frac{\sigma}{\sqrt{n}} $$

• これを$\mu$について整理すると

$$ \bar{X} - 1.96 \times \frac{\sigma}{\sqrt{n}} \leq \mu \leq \bar{X} + 1.96 \times \frac{\sigma}{\sqrt{n}} $$

• ひとつのデータに基づいて行う区間推定では、以下の式が95%の確率で成り立つので、

$$ X -1.96\sigma \leq \mu \leq X +1.96\sigma$$

• 標本平均を用いた区間推定では、ひとつのデータを用いた区間推定に比べると、$\mu$の推定区間が$\frac{1}{\sqrt{n}}$だけ小さくなることが分かる。
• また、$n$を大きくすることで推定区間は小さくなる。

サンプルサイズを増やして推定区間が狭まるか確認してみよう

for i in range(100,10001,100):
    s = np.random.choice(x,i)
    x_mean = np.average(s)
    lower = x_mean - 1.96*(10/np.sqrt(i))
    upper = x_mean + 1.96*(10/np.sqrt(i))
    print(lower,upper)

98.69265877002417 102.61265877002415
99.81147704047244 102.58333562272371
99.41120890564176 101.67442196086509
99.806268724807 101.766268724807
99.2598207889833 101.01289808334313
99.63623239920999 101.23656569782834
99.66526093211515 101.14688166631133
99.57216358287246 100.95809287399808
99.1854165073879 100.49208317405457
99.194793087419 100.43440593020502
99.86691783157336 101.04884230231819
98.85180055226458 99.98340707987624
99.53521522036014 100.6224276049616
99.42602484912207 100.47368891741878
99.86815411967612 100.880293767485
99.68601654838551 100.6660165483855
99.81372098924761 100.76446063939004
99.60095487689475 100.52490773764518
99.46919860323257 100.36850828000516
99.58478160817417 100.4613202553541
99.68854964216128 100.54396377188637
99.93662190845714 100.77236871657112
99.36832961078888 100.1857061138926
99.65080916793632 100.45097581724549
99.67301571040704 100.45701571040703
99.16806211581269 99.93683736555437
99.46439903359442 100.21880338533552
99.50586670163673 100.24667706873481
99.40897869128327 100.13690441693731
99.54428623624466 100.25997704471808
99.92865839821759 100.63271118216254
99.35927578834057 100.05224043390339
99.31738590229413 99.99977031364047
99.58399192246206 100.2562663762207
99.51623930962211 100.17884024532927
99.47958608510535 100.13291941843869
99.7083252096235 100.35276923986049
99.90908996968422 100.54499794051681
99.76127689601225 100.3889792589282
99.74140097489304 100.36120739628605
99.79525439169686 100.40745553830018
99.72994151787816 100.33481064972955
99.703525265826 100.3013196615298
99.5300615012539 100.12102373662631
99.61325319341864 100.19761229153858
99.86083812127859 100.43881058940579
99.82191074442892 100.39370151110066
99.5622462958037 100.12804955960954
99.62869254658557 100.18869254658557
99.96730313212467 100.52167484857492
99.84409541602844 100.39300520896742
99.74676970695518 100.2903758992559
99.68080106589382 100.21925447657269
99.56639827059888 100.09984270347167
99.72980486984873 100.25837756201993
99.80115712630676 100.3249891604551
99.99775969657092 100.51697638054041
99.51146665938177 100.02618787619187
99.76961945807298 100.2799599891178
99.75892625639665 100.26499608030107
99.80351831158252 100.30542288091947
99.73256545208334 100.23040594992408
99.57848157515936 100.07235515322476
99.79869618722005 100.28869618722005
99.67353951411714 100.15975567370684
99.62901845582489 100.11153710046388
99.71449021046985 100.19339443234657
99.69789679691075 100.17326662198198
99.71161464108795 100.18352718518418
99.68475120801995 100.15328082287903
99.4307773270126 99.8959957370246
99.63602014507661 100.09799657545183
99.8615599737334 100.32036127074223
99.81937751300435 100.27506825688013
99.6771462371899 100.12978884823457
99.79317588351373 100.24283072190002
99.84136615137083 100.28809161109263
99.73988575818532 100.18373835557
99.58605667774195 100.02709113490495
99.55853085279439 99.99680017638435
99.67415051379255 100.10970606934812
99.62774620450837 100.06063778672176
99.80610196457413 100.23637786373372
99.71733489233398 100.14504195719653
99.86005772404417 100.28524142136874
99.65604038421225 100.07874485516973
99.84601769883119 100.26628581248761
99.86502076847667 100.28289417253366
99.68530957551377 100.10082874447627
99.71709462056666 100.13029890149534
99.7494180125868 100.16034566858178
99.92859830345017 100.33728655500204
99.79098821863514 100.19747328296947
99.77913737385451 100.18345450238796
99.50740341821226 99.90958693222963
99.9217195028384 100.32180282749299
99.86121477767753 100.25923047440992
99.724841964212 100.12082176167647
99.79062817143485 100.18460299501645
99.73597710042638 100.12797710042638

4.2 標準化

• 標準化とは平均を0、分散（標準偏差）を1にすることでした。
• このプロセスを詳細に見ていきます。
• 以下のデータ(L)は平均が６３、標準偏差が23.8です。

L = [25,30,45,55,65,65,70,80,95,100]

L_ave = np.average(L)
L_std = np.std(L)
print(L_ave,L_std)

63.0 23.790754506740637

• このデータの平均を0にするためにここのデータから平均を引きます。
• L2の平均値は0になります。

L2 = []

for i in L:
    L2.append(i - L_ave)
    
np.average(L2)

0.0

• 標準偏差を標準偏差で割ることにより標準偏差は1になります。

標準偏差の定義は以下でした。 $$ S = \sqrt{\frac{\sum{(x - \bar{x})^2}}{n}} $$

上でデータ全体の平均を0にしたので、上式の$\bar{x}$は0なので、このデータの標準偏差は以下で求められます。

$$ S = \sqrt{\frac{\sum{x^2}}{n}} $$

さて、ここで、すべてのデータを$a$倍したときの標準偏差を考えてみましょう。

$$ \sqrt{\frac{\sum{(ax)^2}}{n}} $$

$\sum$の係数を展開します。

$$\sqrt{\frac{\sum{a^2x^2}}{n}}$$

$a^2$を$\sum$の前に出します。

$$ \sqrt{\frac{a^2\sum{x^2}}{n}} $$

$a^2$をルートの外に出します。

$$ a\sqrt{\frac{\sum{x^2}}{n}} $$

• すべてのデータを$a$倍すると、標準偏差も$a$倍になることがわかりました。
• 標準偏差を$\frac{1}{5}$にするためには個々のデータに$\frac{1}{5}$をかけます。
• 標準偏差を1にするには、標準偏差に$\frac{1}{S}$をかけます。

$$ S \times \frac{1}{S} = 1 $$

• 確認してみましょう。
• 各データを10倍すると標準偏差は10倍になりますが、平均は0のままです。

L3 = []

for i in L2:
    L3.append(i*10)
    
np.std(L3)

237.90754506740637

np.average(L3)

0.0

• 標本平均$\bar{X}$の分布を標準化すると

$$ Z = \frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}} $$

• となり、この$Z$は平均0、標準偏差1の標準正規分布に従うので、95%の確率で以下の式が成り立つ。

$$ -1.96 \leq \frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}} \leq 1.96 $$

• $\mu$について整理すると

$$ \bar{X} - 1.96 \times \frac{\sigma}{\sqrt{n}} \leq \mu \leq \bar{X} + 1.96 \times \frac{\sigma}{\sqrt{n}} $$

• $Z$の値は、母分散が既知の場合は求められるが、未知の場合は母分散を不偏分散で代用した以下の$T$という値を用いる。

$$ T = \frac{\bar{X} - \mu}{\frac{U}{\sqrt{n}}} $$

• この$T$は自由度($n - 1$)に応じたt分布に従うことが知られている。

• 以下は標準正規分布と自由度に応じたt分布である。

• t分布の特徴

  • 左右対称
  • 正規分布よりピークが低い。
  • 正規分布より両端が浮いている。
  • nが大きくなると正規分布に近づく。

• 以下のグラフは$T$の自由度($n - 1$)を1, 5, 10と増やした場合と正規分布との違いを示している。

x = np.arange(-5,5,0.1)
y = []
for i in x:
    y.append((1/np.sqrt(2*math.pi))*math.exp(-(i**2)/2))
    
n_deg = [1,5,10]
for i in n_deg:
    plt.plot(x, t(i).pdf(x),label="n={}".format(i))
    
plt.plot(x,y,label="Normal Distribution")
plt.legend(loc="upper right")

<matplotlib.legend.Legend at 0x12c92dac0>

4.3 区間推定（分散が未知の場合）

問3: ある正規分布に従う母集団($\mu$、$\sigma^2$ともに未知)から、無作為にx人の点数を取り出した。この時、信頼度95%で母集団の平均を推定せよ。

• t分布のパーセント点（信頼区間）は自由度によって異なります。
• 以下は$ n = 3$の例です。

• 自由度3のt分布では95%の確率で以下の式が成り立つ。

$$ -3.18 \leq \frac{\bar{X} - \mu}{\frac{U}{\sqrt{n}}} \leq 3.18 $$

• これを$\mu$について整理すると

$$ \bar{X} - 3.18 \times \frac{U}{\sqrt{n}} \leq \mu \leq \bar{X} + 3.18 \times \frac{U}{\sqrt{n}} $$

• となり、この式が95%の確率で成り立つ。

# 平均が100で標準偏差が10の正規分布から100000回のサンプリング
# これを母集団とみなす。
x = np.random.normal(100, 10, 100000)

自由度は以下のように求めます。

# 自由度3の95%点
t.ppf(0.975,3)

3.182446305284263

サンプルサイズを増やしていって、推定される区間の変化を見てみましょう。

for i in range(100,10001,100):
    s = np.random.choice(x,i)
    x_mean = np.average(s)
    x_std = np.std(s,ddof=1)
    t_p = t.ppf(0.975, i - 1)
    lower = x_mean - t_p*(x_std/np.sqrt(i))
    upper = x_mean + t_p*(x_std/np.sqrt(i))
    print(lower,upper)

98.14910023280373 102.23982169151702
98.28159318047524 101.23078563281125
99.7548039238907 102.18138279539349
99.19798068603406 101.18986690619211
98.95932453037616 100.62122393169906
98.92915710301428 100.50199507188417
99.05856895229262 100.50620986689667
99.35304637148984 100.7493603781924
99.31155668974078 100.64577916515124
99.20938891725346 100.48681277260751
99.3801041112124 100.55046115899808
99.3206892695063 100.4540434507782
99.72303424845215 100.79963718579317
99.41306663513036 100.45936478409847
99.53081093098749 100.52480675502952
99.84506518835778 100.83462971455208
99.3126984137373 100.24654533856786
99.45571162547667 100.3849717896815
99.5360707705797 100.41808158627802
99.38441843963528 100.23395404931543
99.32837371253947 100.1780366822342
99.68776083546094 100.51892127962641
99.22756518686562 100.04505095993373
99.59954571790946 100.40432027734481
99.29156149620947 100.06286083538977
99.67532026362868 100.44555069320752
99.55071762190202 100.28540773744903
99.56333443830378 100.31503285508238
99.88087242656731 100.59766587960435
99.74219254084973 100.47395401273603
99.60327453034962 100.30912010588374
99.44325073488514 100.13008030204188
99.60301039815015 100.2862636861821
99.50532253627128 100.17658504464816
100.00151090676307 100.66675656810837
99.82808890492933 100.48541695363303
99.5587792066063 100.20334401924278
99.56423458841955 100.18611979602308
99.72382197683307 100.34524385848322
99.49659426504974 100.11548607635865
99.70765338982119 100.3169174549422
99.81621829118212 100.4241802047327
99.71339167482014 100.3077098411419
99.50211634557108 100.10103923350634
99.63035453440568 100.21530019304825
99.8339065462094 100.4149783851769
99.6474703820871 100.23077818525539
99.85021887078031 100.41589580912677
99.70534645763996 100.26130337193288
99.92465825083886 100.46758592064427
99.84896897675529 100.3876627202317
99.67504294131221 100.21501617098964
99.76601813521444 100.30190884085226
99.82347363559371 100.36058117111146
99.72763518333072 100.25468985104307
99.96299638325999 100.47864273468176
99.62789197384623 100.15005488321006
99.63126807808658 100.15256232085005
99.66629925154787 100.17352132285792
99.72628126850229 100.23589434005784
99.67003845852894 100.17378765159002
99.55850001489539 100.04634563691282
99.68433266603087 100.17975003112083
99.7256790955022 100.21981591326251
99.68957920687139 100.17909662964321
99.77424617644388 100.26234412727081
99.85006219865011 100.33258628854458
99.9616237389137 100.43512269719498
99.57040335474795 100.05116054797179
99.86353837603963 100.34004348460277
99.74497820813662 100.21419137381758
99.916196912511 100.37815269071189
99.70207162283872 100.16079375493977
99.75823424612318 100.22123398600615
99.94254837991093 100.39395188515238
99.84196792945515 100.29209898167586
99.69540545102413 100.1437371408883
99.62954005833016 100.06902811056749
99.83152486479004 100.27531996803576
99.84686202338212 100.28032724806589
99.75467754271877 100.18364943037305
99.65969967861712 100.09585349995378
99.91069776268026 100.33886117803999
99.64863834411938 100.08492671286353
99.8682645714727 100.28992388645052
99.88078004326444 100.30693173006667
99.9075213039425 100.32785640479216
99.80943021218062 100.22952040951301
99.98947884177008 100.40807590813839
99.65839360533789 100.06778072932144
99.907417088514 100.3219818407542
99.77828948146005 100.18286252340049
99.66461496590215 100.072336244563
99.62850231640091 100.03213917319557
99.7064000203679 100.10659384156175
99.87361343288437 100.27669550003087
100.00154601304686 100.39984203857347
99.80169162915864 100.19550712886014
99.97307782899748 100.36325256174467
99.86286054769322 100.25349032380022

4.4 中心極限定理

• 分布が分からない母集団から無作為に抽出した標本の$n$が十分に大きいとき、以下の式が成り立つという定理。

$ \sum_{n=1}^i X_i $ ~ $N(n\mu, n\sigma)$

$ \bar{X} $ ~ $ N(\mu, \frac{\sigma}{\sqrt{n}}) $

• 以下ではサイコロを10000回振った。

D1= {1:0,2:0,3:0,4:0,5:0,6:0}

for i in range(10000):
    d= np.random.randint(1,7)
    D1[d] += 1

left = [1,2,3,4,5,6]
height = ([D1[1],D1[2],D1[3],D1[4],D1[5],D1[6]])
plt.bar(left,height,width=1.0,edgecolor="k")

<BarContainer object of 6 artists>

• サイコロのそれぞれの目が出る確率は等確率であるが、サイコロを1000回降ってその平均値を求めるという行為を1000回繰り返し、平均値を算出するとその平均値の分布は正規分布と近似する。

D = []

for i in range(1000):
    tmp = []
    for i in range(1000):
        tmp.append(np.random.randint(1,7))
    D.append(np.average(tmp))

plt.hist(D,bins=25)
plt.show()

4.5 区間推定（母集団分布が未知の場合）

問4: 分布が未知である母集団から、無作為にx人の点数を取り出した。この時、信頼度95%で母集団の平均を推定せよ。

• $n$が十分大きいとき、$\bar{X}$は正規分布に従うので、95%の確率で以下の式が成り立つ。

$$ \bar{X} - 1.96 \times \frac{\sigma}{\sqrt{n}} \leq \mu \leq \bar{X} + 1.96 \times \frac{\sigma}{\sqrt{n}} $$

• また、$n$が十分大きいとき、母分散$\sigma$を不偏分散$U^2$を用いた標準偏差$U$で代用できるので、以下を95%信頼区間とする。

$$ \bar{X} - 1.96 \times \frac{U}{\sqrt{n}} \leq \mu \leq \bar{X} + 1.96 \times \frac{U}{\sqrt{n}} $$

# 平均が100で標準偏差が10の正規分布から100000回のサンプリング
# これを母集団とみなす。
x = np.random.normal(100, 10, 100000)

• $n$を30とした場合

L = np.random.choice(x,30)

L_mean = np.mean(L)
L_std = np.std(L,ddof=1)

lower = L_mean - 1.96 * L_std/np.sqrt(30)
upper = L_mean + 1.96 * L_std/np.sqrt(30)

print(lower,upper)

98.45679598660443 105.78097200696389

練習問題1

正規分布に従う母集団（$\mu$が未知、$\sigma^2=200$)から30人の生徒の点数を取り出したとき、その平均点は115点であった。信頼度95%で母集団の平均を推定せよ。

練習問題2

ある正規分布に従う母集団($\mu$、$\sigma^2$ともに未知)から、無作為に10人の点数を取り出した(X)。この時、信頼度95%で母集団の平均を推定せよ。

X = [95,80,65,50,90,35,40,70,70,55]

練習問題3

分布が未知である母集団から、無作為に20人の点数を取り出した(X)。この時、信頼度95%で母集団の平均を推定せよ。

X = [95,80,65,50,90,35,40,70,70,55,50,750,90,50,90,25,20,75,70,65]