10. 推測統計学

• このセクションで使う R パッケージ一覧

library(tidyverse)
library(patchwork)

• 統計学は記述統計学 (descriptive statistics) と推測統計学 (inferential statistics) の二つに分類できる 

• 記述統計学では、平均や分散（データの散らばり）などを使って、データの特徴を整理したり記述する 

• 推測統計学では、母集団 (population) から標本 (sample) を無作為に抽出し、その標本によって得られた標本平均や不偏分散などの統計量を使って、母集団の母数（母平均や母分散）を検定し推定する 

• 母集団の全てを調査すること・・・全数調査（センサスまたは悉皆 [しっかい] 調査）
  

• しかし、現実的に全数調査を行うことが出来ないことが多い → 母集団の一部を切り取る（標本：サンプル）

• 推測統計学の基本的な考えは、母集団からランダムに抽出した標本を増やし、無限に試行を繰り返せば、全体の一部である標本から、巨大で未知の母集団を推測できる、ということ（= 統計的推論）　　
  

• 推測統計学では、分析対象が確率分布すると考える 

• ここでは推測統計学の基礎を演習する
  

1. 推定と検定

• 母集団から無作為抽出された標本を使って、抽出元全体である母集団の特徴を推測する統計学を推測統計学 (inferential statistics) という 
• 推測統計学は「推定」と「検定」に分類できる 

• 推定 (estimation)・・・具体的な値を用いて「母集団の平均値は○○くらいだろう」と予想すること 

1. 点推定・・・標本から判断して、具体的な母数の値を予想する  ex: 「早稲田大学全体で女子学生の平均身長は 160 cm くらいだろう」
   
2. 区間推定・・・標本から判断して（ある程度の幅を持った区間で）具体的な母数の値を予想する  ex: 「早稲田大学全体で女子学生の平均身長は 160 cm から 165 cm くらいだろう」

• 検定 (test)・・・母集団について述べた異なる 2 つの「仮説」のうちどちらかを採用すること
  

2. 母集団と標本

2.1 「標本数」と「標本サイズ」の違い

• 母集団 (population): 明確に定義された分析対象
• 標本 (sample): 母集団から抜き出される一部
  
• 標本は母集団のミニチュア版

• 推測統計学においては、標本サイズ (sample size) と標本数 (the number of samples) が混同されがちなので、十分注意が必要
  

• 標本サイズ (sample size)・・・・・・・n で表す抽出した観測数のこと
  

• 標本数 (the number of samples) ・・・ n 個の観測数で抽出した標本セット数

例：　　

• 矢内さんは人口328万人の郡山市に行って、東日本大震災に関して二つの標本をとった 
• 標本 A では 1,000 人を対象に世論調査を実施し、標本 B では 2,000 人に対して世論調査を行った 

→ 標本数は 2（標本 A と標本 B）であり、標本サイズはそれぞれ 2,000 と 1,000 である　　　　

母数と統計量:
推測統統計学では、母数 (parameter) と統計量 (statistic) とを明確に区別しているので注意が必要

母数

• 母平均: \((μ)\)（ミュー）— population mean
  
• 母比率: \((π)\)（パイ）— population ratio
  
• 母分散: \((σ^2)\) — population variance
  
• 母標準偏差: \((σ)\)（シグマ）— population standard deviation

標本統計量

• 標本平均: \(\bar{x}\) — sample mean
  
• 標本比率: \((p)\) — sample ratio
  
• 標本分散: \((s^2)\) — sample variance
  
• 標本不偏分散: \((u_x^2)\) — unbiased variance
  
• 標本標準偏差: \((s)\) — sample variance
  
• 標本不偏標準偏差: \((u_x)\) — unbiased standard deviation

2.2無作為抽出法 (Random Sampling)

(1) 単純ランダム・サンプリング法 (simple random sampling):

• 偏りのない標本を選ぶ代表的な方法
  
• 単純無作為抽出された標本は、母集団の偏りのない縮図と見なせる
  → その標本を使って、母集団を推定できる

(2) 多段抽出法:

• 何段階かに分けて標本の抽出を行う方法
  →例: 日本の有権者から1000人を標本として抽出する
  1 段階: 市町村を無作為に選ぶ
  2 段階: 選ばれた市町村から無作為に標本を選ぶ
  
• 母集団分布が広範囲に分布している場合に多段抽出法を使う
  →多段抽出法の一例・・・層化 2 段抽出
  
• 各市町村を大都市、中小都市、郡部と「層化」する
  
• 各層の大きさに応じたウェイ卜（例えば 3 : 2 : 1 ) をつけて市町村を無作為に選ぶ

【母平均 \((μ)\) と標本平均 \(\bar{x}\) の関係】

• 母集団に関する母数は未知のことが多いが、ここでは R を使って人工的に母集団を作成してみる
  

• 作成する母集団は、母平均 10、母標準偏差 5（母分散 = 25） 

• dnorm()関数を使うと、母集団を人工的に作成でき、その分布密度を表示できる 

# 母平均 10、母標準偏差 5 と指定
# -10 から 30 の範囲を指定
curve(dnorm(x, 10, 5), from = -10, to = 30) 

• この母集団からランダムに標本を抽出する 
  

• 抽出する標本サイズは 20 (N = 20)で、その標本に \((x1)\) と名前を付ける 

• ここで人工的に作成した母集団からランダムサンプリングする時は rnorm() 関数を使う 

# 平均10、標準偏差 5 (sd = 5) の母集団から 20 個の数字を無作為抽出して x1 と名付け、データとヒストグラムを表示する

x1 <- rnorm(20, mean = 10, sd =5) 
x1

 [1]  7.191427 17.707697  9.749040 12.720787 15.094785 15.396032 10.072061
 [8]  6.798442  8.820112  6.858025 10.179139  3.921366  4.517479 13.924562
[15] 10.191817  4.302003 18.414976  9.708356  7.567451  6.455296

• 表示したい桁数は round() 関数を使って指定できる  例えば、整数を表示したければ、小数点以下の表示桁数 (digits) を 0 と指定すればよい 

round(x1, digits = 0)　# 小数点以下表示したい桁数を指定

 [1]  7 18 10 13 15 15 10  7  9  7 10  4  5 14 10  4 18 10  8  6

hist(x1)　　　　　　　# 取り出したサンプルをヒストグラムで表す

標本 \((x)\) の標本平均 \(\bar{x}\) は次の式で求めることができる: \[\bar{x} = \frac{\sum_{i=1}^n x_i}{n}\]

・R では mean( ) を使えば標本平均を求めることができる

mean(x1)              # 取り出したサンプル x1 の平均値

[1] 9.979543

2.3 「推定量」と「推定値」の違い

\[\bar{x} = \frac{\sum_{i=1}^n x_i}{n}\]

• 推定量 (estimator) :
  母数の推定のために利用され、標本の関数として表現される・・・ここでは上記 \(\bar{x}\) の式
  

• 推定値 (estimate) :
  実際に標本の数値を代入して計算した結果の数値・・・ここでは上記 \(\bar{x}\) の 値

• 標本サイズ 20 の標本を無作為抽出して x2 と名付け、標本平均を求めてみる
  

• 母集団の平均は 10 (mean = 10)、標準偏差は 5 (sd = 5) と指定

x2 <- rnorm(20, mean = 10, sd =5) 
x2

 [1]  2.2102410  5.6472328 18.7305019 16.6984417  4.6051185  8.3373866
 [7] 10.1534012 10.3106054  8.3268239 12.1513965  4.8177404  6.1054748
[13] 13.6114709 10.4737223 -0.6118851 12.5638541  6.2480742 10.8739064
[19]  1.1506101 20.2324380

round(x2, digits = 0)      # 小数点以下表示したい桁数を指定

 [1]  2  6 19 17  5  8 10 10  8 12  5  6 14 10 -1 13  6 11  1 20

hist(x2)                 　# 抽出したサンプルをヒストグラムで表す

mean(x2)                  # 取り出したサンプル x2 の平均値

[1] 9.131828

• R は母集団から無作為に数値を抽出しているので、サンプリングにおいて「偏り」はないが「誤差」はある 
• ここでは母平均を 10 に設定したにもかかわらず、一回目のサンプル x1 と二回目のサンプル x2 の標本平均がぴったり 10 にはならない 
• しかし、サンプリングを何度も繰り返すと（つまり x1, x2, ..... xn)、母平均より大きかったり小さかったりする標本平均が得られ、それらの値を平均すると母平均に一致する 
• → 標本から得られた統計量を使って母集団のパラメタを統計的に推定できる根拠 
• 標本平均以外にも min( ), median( ), max( )などのコードを使うと最小値、中央値、最大値などの記述統計が簡単に得られる 

min(x2)　　# min = 最小値

[1] -0.6118851

median(x2) # median = 中央値

[1] 9.245394

max(x2)    # max = 最大値

[1] 20.23244

• 5数は次のコマンドで一度に表示できる
  

quantile(x2, c(0, .25, .5, .75, 1))

        0%        25%        50%        75%       100% 
-0.6118851  5.4398597  9.2453939 12.2545109 20.2324380 

3. 中心極限定理

3.1 中心極限定理の定義

• 中心極限定理 (CLT: Central Limit Theorem) とは、母集団が正規分布していなくても「標本サイズ（サンプルサイズ）が大きくなれば、標本平均は正規分布で近似できる」という定理
  
• 標本サイズを無限大にすると標本平均を標準化したものが標準正規分布に収束する」という定理
  
• 統計的推定や検定を行うときに標準正規分布を使えるのはこの定理のおかげである
  
• この定理のおかげで、世論調査などで統計的推定をするために必要な標本サイズを算出することが可能 

正規分布 (Normal Distribution)

• 「平均値の付近に集積するようなデータの分布を表した連続的な変数に関する確率分布」
  

• 様々な分布を正規分布に見立てると、その内容の範囲を確率的に予想することが可能になる 

• 例えば、95%の信頼区間が \(([\bar{x} - 1.96\mathrm{SE}, \bar{x} + 1.96\mathrm{SE}])\) （SEは標準誤差で、標本平均の標準偏差のこと）の間に存在することを推定するために「正規分布では母平均 \((\mu)\) を中心に±1.96標準偏差の中に全体の 95% が収まる（\((\mu \pm 1.96 \sigma)\)）」という事実を使えるのは中心極限定理の裏付けがあるから 

• ここでは R を使ってシミュレーションを実行し、中心極限定理の信憑性を「実感」することが目的 

3.2 中心極限定理のシミュレーション

標本サイズが大きくなると、標本平均は近似的に正規分布に従うのか？

バックに入れた10枚のカードを使ったシミュレーション

• 0 から 9 までの数字が書かれた 10 枚のカードを作成しバックに入れる 
• このバックを母集団だとする（人工的に母集団を作成） 
• 母集団から一枚カードを無作為に選ぶ 
• それぞれのカードが選ばれる確率は等しく 1 / 10 
• カードを無作為に選んだ時の数の分布は離散一様 (discrete uniform) 
• 母集団の10 枚のカードの母平均は 4.5 なので、標本平均も 4.5 になるのか？
  
• 標本平均から、通常未知である（ここでは既知である）母平均を推計する試み 
• カードは一枚ずつ引き、引いた番号を記録し、カードをバッグに戻す 
• カードをシャッフルして、再びカードを一枚引きその番号を記録し、カードをバックに戻す
• 無作為に選んだ 2 枚のカードの平均点を記録する
  
• 以上の作業を n 回繰り返す
  

• 1 回目のカードの選び方は 10 通り （0から9まで10枚カードがあるから）
  
• 2 回目のカードの選び方は 10 通り （0から9まで10枚カードがあるから）
  →10枚のカードを 2 枚引いた時の 2 枚のカードの組み合わせは 100通り
  
• カードに書かれている数は 0 から 9 までの 10 個の整数
  →可能な合計値は 0 から 18 までの 19 通り（例：選ばれた数字が 2 枚とも 9 なら合計値は 18）
  →可能な平均値は 0 から 9 までの 19 通り（例：選ばれた数字が 2 枚とも 9 なら平均値は 9）
  
• 平均値が得られる確率は次のように表すことができる

出典：https://yukiyanai.github.io/jp/classes/stat2/contents/R/rng-n-clt.html

• 例えば、ここで想定したような 0 と 8 の 2 枚のカードを引いた時の標本平均 (4) が得られる確率は 0.8 (= 8%) 程度だとわかる
  
• 母平均である 4.5 が得られる確率は 0.1 (= 10%) である
  
• R を使って、実際に 1 回引いてみる

bag <- 0:9                    # 人工的に 0 から 9 までのカードを作成  
exp_1 <- sample(bag,          # カードが入っている bag を使う  
              size = 2,       # 引くカードは 2 枚と指定（標本サイズ N = 2）
              replace = TRUE) # 取り出したカードは元に戻すと指定  
mean(exp_1)                   # 平均値を計算して表示  

[1] 8.5

• 実験をやる度に異なる平均値が得られる（何回か試して下さい）

• この実験を 1,000 回繰り返してみよう
  

• 「それぞれの回で得られた平均値の分布」はどんな形か確かめる
  

• R で同じ作業を繰り返し行う簡単な方法は、for ループを使う 　　

for ループの使い方：

• 例えば、0 からスタートして「10を足す」という作業を 5 回繰り返すには、次のようにする

• まず、スタート時の数である 0 を作成し、保存する

A <- 0 # スタート時の数である A を 0 と指定し保存  

• 次に、結果を保存するために、result という名前の入れ物を用意する

result <- rep(NA,              # NA は「中身が空」であることを表す
              length.out = 5) # 保存場所として何カ所必要かを指定  

• result の中身を確認してみる

result

[1] NA NA NA NA NA

• ここでは入れ物を作っただけなので、中には何も入っていない (= NA)
  
• forループを使って、数を 10 だけ加える作業を 5 回繰り返してみる
  
• for の直後の丸カッコ ( ) の中に繰り返し回数を指定し、ループさせる内容を中括弧 \(\{ \}\) で囲む

for(i in 1:5){   # i が 1 から 5 まで繰り返すという意味  
  A <- A + 10    # A に 10 を足す
  result[i] <- A # i 番目の作業（＝足し算）の結果を result[i] に入れる  
}
result           # 入れ物の中身を確認

[1] 10 20 30 40 50

→入れ物 (result) の中に想定どおりの数値 （0 からスタートして 10ずつ足した数） が入っていることがわかる

for ループの使った中心極限定理のシミュレーション：
sim1: 標本サイズ = 2, 繰り返し回数 = 10

• このfor ループを利用して、皆さんが授業の実習で実施したように、カードを 2 回引く作業を 10 回繰り返してみる

N <- 2                                # 標本サイズ（カードを 2 枚選ぶ）
trials <- 10                          # 実験の繰り返し回数 = 10
sim1 <- rep(NA, length.out = trials)  # 結果の保存容器
for (i in 1:trials) {
  experiment <- sample(bag, 
                       size = N, 
                       replace = TRUE) # 復元抽出
  sim1[i] <- mean(experiment)          # i 回目の平均値を保存
}

• 結果をヒストグラムにしてみる

文字化けしないためのコマンド（マックユーザのみ）

theme_set(theme_bw(base_size = 14,base_family = "HiraKakuPro-W3"))

df_sim1 <- tibble(avg = sim1)
h_sim1 <- ggplot(df_sim1, aes(x = avg)) +
  geom_histogram(binwidth = 1, 
                 boundary = 0.5,
                 color = "black") +
  labs(x = "2枚のカードの平均値", 
       y = "度数") +
  ggtitle("SIM1: N = 2, 繰返回数 = 10") +
  scale_x_continuous(breaks = 0:9) +
  geom_vline(xintercept = mean(df_sim1$avg),  # 平均値に縦線を引く  
             col = "lightgreen")

plot(h_sim1)

• 正規分布には見えない

df_sim1

# A tibble: 10 x 1
     avg
   <dbl>
 1   7.5
 2   1.5
 3   4  
 4   5.5
 5   9  
 6   7  
 7   1.5
 8   1.5
 9   7  
10   3  

• この実験で得られた平均値を知りたければ次のコマンドを実行すればよい

mean(df_sim1$avg)

[1] 4.75

sim2: 標本サイズ = 5, 繰り返し回数 = 100

• カードを 5 枚引く作業を 100 回繰り返してみる

N <- 5                                # 標本サイズ（カードを 5 枚選ぶ）
trials <- 100                          # 実験の繰り返し回数 = 100
sim2 <- rep(NA, length.out = trials)   # 結果の保存容器
for (i in 1:trials) {
  experiment <- sample(bag, 
                       size = N, 
                       replace = TRUE)  # 復元抽出
  sim2[i] <- mean(experiment)         # i 回目の平均値を保存
}

• 結果をヒストグラムにしてみる

文字化けしないためのコマンド（マックユーザのみ）

theme_set(theme_bw(base_size = 14,base_family = "HiraKakuPro-W3"))

df_sim2 <- tibble(avg = sim2)
h_sim2 <- ggplot(df_sim2, aes(x = avg)) +
  geom_histogram(binwidth = 1, 
                 boundary = 0.5,
                 color = "black") +
  labs(x = "5 枚のカードの平均値", 
       y = "度数") +
  ggtitle("SIM2: N = 5, 繰返回数 = 100") +
  scale_x_continuous(breaks = 0:9) +
  geom_vline(xintercept = mean(df_sim2$avg),  # 平均値に縦線を引く  
             col = "lightgreen")

plot(h_sim2)

• 正規分布には見えない

df_sim2

# A tibble: 100 x 1
     avg
   <dbl>
 1   5.8
 2   5.2
 3   3.8
 4   6  
 5   2.6
 6   3  
 7   5.2
 8   4.2
 9   6.2
10   3.2
# … with 90 more rows

• この実験で得られた平均値を知りたければ次のコマンドを実行すればよい

mean(df_sim2$avg)

[1] 4.602

sim3: 標本サイズ = 100, 繰り返し回数 = 1000

• カードを 100 枚選ぶ作業を 1000 回繰り返してみる

N <- 100                             # 標本サイズ（カードを 100 枚選ぶ）
trials <- 1000                       # 実験の繰り返し回数
sim3 <- rep(NA, length.out = trials) # 結果の保存容器
for (i in 1:trials) {
  experiment <- sample(bag, size = N, replace = TRUE)  # 復元抽出
  sim3[i] <- mean(experiment)        # i 回目の平均値を保存
}

• 結果をヒストグラムにしてみる

df_sim3 <- tibble(avg = sim3)
h_sim3 <- ggplot(df_sim3, aes(x = avg)) +
  geom_histogram(binwidth = 0.125, 
                 color = "black") +
  labs(x = "100枚のカードの平均値", y = "度数")+
  ggtitle("SIM3: N = 100, 繰返回数 = 1000") +
  scale_x_continuous(breaks = 0:9) +
  geom_vline(xintercept = mean(df_sim3$avg),  # 平均値に縦線を引く  
             col = "lightgreen")

plot(h_sim3)

df_sim3

# A tibble: 1,000 x 1
     avg
   <dbl>
 1  4.72
 2  4.94
 3  4   
 4  4.83
 5  4.85
 6  4.59
 7  4.87
 8  4.23
 9  4.99
10  4.26
# … with 990 more rows

• この実験で得られた平均値を知りたければ次のコマンドを実行すればよい

mean(df_sim3$avg)

[1] 4.51027

• N = 2 と N = かなり正規分布に近くなったことがわかる

library(patchwork)

h_sim1 + h_sim2 + h_sim3

まとめ ・元になる母集の分布は一様分布でも、サンプルサイズ N を増やすと「平均値の分布」は正規分布に近づく

→サンプルサイズ N が十分大きい（大まかな目安はN >100）とき、正規分布を使って統計的推定や検定を行うことが許される

4. 標本分布

4.1 統計量は分布する？

• 標本分布 (sampling distribution) : - ある標本サイズで標本を抽出したとき、そこで得られる統計量の分布
  
• ランダムサンプリングにおいて、標本を選ぶ方法はたくさんある 
• 例えば、10 人から 3 人を選ぶ方法 \((_{10}C_3)\) は 120 通りある 
• R では choose(　) コードを使って次のように計算できる 

choose(10, 3)  # 10 人から 3 人を選ぶ組み合わせ

[1] 120

• 母集団が 1 億人の有権者から 10 人の標本を選ぶ方法を計算してみると、\((2.7)\) x \((10^{73})\) 通り以上あることがわかる 

choose(100000000, 10)

[1] 2.755731e+73

• 特定の母集団について特定のサイズの標本を選ぶ方法は何通りもある 
• それぞれの標本について統計量の値を求めることは可能 
• しかし、統計量の値が全て一致するとは限らない 
• 同じ大きさの標本を抽出する作業を繰り返せば、統計量は分布する 
• 一定の標本サイズで全ての組み合わせで標本を抽出したとき、そこで得られる統計量の分布を標本分布と呼ぶ 
• 標本分布・・・統計量と母数を確率でつなぐ統計的推定の肝 

4.2 「母平均」と「標本平均」の関係

• 母集団を人工的に作る
  
• 3人の議員（Aさん、Bさん、Cさん）それぞれの当選回数 

・当選回数 3, 6, 12回の母集団を作成し、母平均と母標準偏差を計算する 

win <- c(3, 6, 12)  # 当選回数 3, 6, 12回の母集団を作成
mean(win)　　　　　 # 当選回数の母平均

[1] 7

• 分散を計算する

var1 <-  ((7-3)^2 +(7-6)^2 + (7-12)^2)/3
var1

[1] 14

• 標準偏差は分散の平方根なので、標準偏差は

sqrt(var1)

[1] 3.741657

• 平均当選回数は 7 回、当選回数の標準偏差は 3.7
  

• 従って、母平均 \((μ = 7)\)、母標準偏差 \((σ = 3.7)\)
  

• この母集団からサイズ = 2 の標本を全て選ぶ 

• A, B, C の 3 人から 2 人が選ばれるのは次の 3 通り：「A と B」、「A と C」、「B と C」 
  

• 例えば「選ばれる人物」が「A と B」の当選回数に関する標本平均と標本標準偏差を R で計算してみると、

ab <- c(3, 6)  # A と B が選ばれたサンプルを ab と名付ける
mean(ab)　　　 # サンプル ab の標本平均

[1] 4.5

• 分散を計算する

var2 <-  ((4.5-3)^2 + (4.5-6)^2)/2
var2

[1] 2.25

• 標準偏差は分散の平方根なので、標準偏差は

sqrt(var2)

[1] 1.5

• 上の表の「標本平均」を見ると、3 つの標本の標本平均 (4.5, 7.5, 9) が全て異なっている 
• しかも、事前に設定した母集団の母平均 \((μ = 7)\)とも異なる 
• 標本から母数について推定できること・・・標本が母集団の偏りのない縮図であること
  
• 「偏りのない」＝ひとつひとつの統計量の値が母数と異なっても、得られた値の平均値が母数に一致する
• 3 つの標本平均の平均値 \((μ)\) を求めてみると

\[{μ} = \frac{4.5+7.5+9}{3} = 7\]

• 抽出可能なすべての組み合わせで標本を抽出し、各標本から得られた平均値の平均値を求めると、母平均に一致する 
• → 不編性 (unbiasedness)　　
• → 不編推定量 (unbiased estimator): 不偏性のある推定量のこと
• → 「標本平均 \(\bar{x}\) は母平均 \((μ)\) の不偏推定量であるといえる」
  → 母平均 \((μ)\) は標本平均 \(\bar{x}\) を使って推定する 

4.3 「母標準偏差」と「標本標準偏差」の関係

• しかし、3 つの標本標準偏差の平均値を求めてみると、母標準偏差 (3.7) に一致しない 

\[\frac{1.5+4.5+3}{3} = 3 < {3.7}.\]

• 標本標準偏差の平均は、常に母標準偏差より小さい 
• 標本標準偏差（統計量）は母標準偏差（母数）を過小評価するという体系的な bias（偏り）をもつ 
• → 「標本標準偏差 \((s)\)は母標準偏差 \((σ)\) の不偏推定量とはいえない」 
• → 標本不偏標準偏差 \((u_x)\)を使って母標準偏差 \((σ)\) を推定する 
  → 母標準偏差 \((σ)\) は標本不偏標準偏差 \((u_x)\) を計算できる
  
• sd( ) 関数 を使って推定する 

sd(win) # win の不偏標準偏差を求める

[1] 4.582576

母分散 \((\sigma^2)\) を求める式・・・（母標準偏差 \((σ)\) は \((σ^2)\) の平方根）

\[\sigma^2 = \frac{\sum_{i=1}^N (x_i - \mu)^2}{N}\]

不偏分散 \((u_x^2)\) を求める式・・・（不偏標準偏差 \((u_x)\) は \((u_x^2)\) の平方根）

\[u_x^2 = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{N-1}\]

• R を使って計算する時に使うコード:

母分散 \((\sigma^2)\) ・・・var(x) - (length(x) - 1) / length(x)

母標準偏差 \((\sigma)\) ・・・sqrt(var(x) - (length(x) - 1) / length(x))

不偏分散 \((u_x^2)\) ・・・var(x)

不偏標準偏差 \((u_x)\) ・・・sd(x)

4.5 大きな母集団を想定したシミュレーション

• 人工的に日本人成人男性の身長に関する正規分布する母集団（n = 10万人、母平均 = 170cm、母標準偏差 = 6）を作り、height と名付ける 

height <- rnorm(100000, 170, 6)
hist(height)

• 10 万人から 10 人選ぶ方法は、\((2.7)\) x \((10^{43})\) 通り以上あることがわかる 

choose(100000, 10)

[1] 2.754492e+43

• \((2.7)\) x \((10^{43})\) 通りの組み合わせで標本平均を求めるのは不可能 
• 代替策 → 10 万人の母集団から無作為に 10 人を選び sample1 と名付け、標本平均を求め、ヒストグラムで示す 

sample1 <- sample(height, 10, replace = FALSE) # 標本サイズ 10 で無作為抽出する
hist(sample1)　　　　　　　　　　　　　　　　　# ヒストグラムで表す

mean(sample1)　　　　　　　　　　　　　　　　  # sample1 の平均値を表示する

[1] 164.0849

• 母集団から 10 人を無作為に選ぶという作業を 500 回行い、その分布を調べる 
• 標本サイズを 10 と 90 の二種類 (n = 10, 90)、標本数は 2 つ作成し、アニメーションで表示する 

アニメーションを使ったシミュレーション (n = 10, 90)

• 人工的に日本人成人男性の身長に関する母集団（n = 10万人、母平均 = 170cm、母標準偏差 = 6）を作り、height と名付ける 

height <- rnorm(100000, 170, 6)  
## Parameters: n = 標本サイズ
##              T = 標本数
height_experiment <- function(n = 10, T = 500, seed = 2016-01-24) {
    means <- rep(NA, T)
    s_mat <- matrix(NA, ncol = n, nrow = T)
    set.seed(seed)
    height <- rnorm(100000, 170, 6)
    for (i in 1:T) {
        s <- sample(height, n, replace = TRUE)
        s_mat[i,] <- s
        means[i] <- mean(s)
        par(family = "HiraKakuPro-W3", cex.lab = 1.4, cex.axis = 1.2, cex.main = 1.6, mex = 1.4)
        hist(means[1:i], axes = FALSE, freq = FALSE, col = "gray",
             xlim = c(min(height), max(height)),
             xlab = "sample mean", ylab = "probability density",
             main = paste0("sample size = ", n, ", sample id = ", i))
        axis(1, min(height):max(height))
        abline(v = 170, lwd = 2, col = "red")
    }
    axis(2)
    return(s_mat)
}

• シミュレーションの結果を表示する際には、次のコマンドを入力する
  

source("<path-of-the-directory>/height_sim.R")
sim_1 <- height_experiment(n = 10) # サイズ 10 の標本を 500 個抽出
sim_2 <- height_experiment(n = 90)　# サイズ 90 の標本を 500 個抽出

• 母集団からサイズ 10 の標本を 500 個抽出した時の標本平均 \(\bar{x}\) の分布

• 母集団からサイズ 90 の標本を 500 個抽出した時の標本平均 \(\bar{x}\) の分布

1. どちらの標本においても、ヒストグラムの中心が母平均 \(({μ} = 170)\) に一致  → どんなに母集団が大きくても、標本平均は母平均の普遍推定量である 
2. 母集団のばらつきより、標本平均のばらつきが小さい  その理由 → 標本平均では極端な値を取りにくいから 

• 標本平均は母集団の分布の両端付近の値をとりにくく、母集団よりも狭い範囲に分布する
• 分布のばらつき・・・標準偏差で測る 
• 標本平均のばらつき・・・次の式で表す
  

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

• 式の分母に \(n\) が含まれているので、標本サイズ \(n\) が大きくなるほど、標本平均はますます極端な値を取りにくくなる 
• 母集団の身長の母標準偏差 \(({\sigma = 6})\)、抽出した 500 個の標本サイズは 10 だから、

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

\[= \frac{6}{\sqrt{10}}\] \[= 1.9\]

• 母集団の身長の母標準偏差 \(({\sigma = 6})\)、抽出した 500 個の標本サイズを 90 に増やすと、

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

\[= \frac{6}{\sqrt{90}}\] \[= 0.63\]

• 標本サイズは 10 のときの標準偏差 1.90 のおよそ3 分の 1 
• 一致性 (consistency): 標本サイズを大きくするほど、ひとつの標本から計算される標本平均が母平均近くの値をとる確率が大きくなる
  
• → 一致推定量 (consistent estimator)

5 母平均の推定と信頼区間

5.1 \(z\) 値と \(t\) 値

• 点推定 (point estimation): 例）「母平均は 45 点」
  

• 区間推定 (interval estimation): 例）「母平均は 43 点以上、47 点以下」
  

• 信頼区間 (confidence interval: CI)・・・推定に伴う不確実性を表すために利用される 

• 標本平均の標準偏差 \(({SD (\bar{x}})\)）がわかる
  → 標本平均の不確実性を捉えることができる
  → 推定の不確実性を明らかにできる
  

• 標本平均の標準偏差 \(({SD (\bar{x}})\)）を求める式は、

\[SD (\bar{x})= \frac{\sigma}{\sqrt{n}}\]

• 母標準偏差 \({\sigma}\) と 標本サイズ \(n\) の値がわかれば、標本平均の標準偏差 \(({SD (\bar{x}})\)）がわかる 

• 標本サイズ \(n\) はわかるが、母標準偏差 \({\sigma}\) は通常わからない  → 母標準偏差 \({\sigma}\) の代わりに不偏標準偏差 \(u_x\) を使って、標本平均の標準偏差 \(({SD (\bar{x}})\)）を推定 
  → このページの「標本分布」の項目を参照

\[SE = \frac{u_x}{\sqrt{n}}\]

• この \(SE\) が標準誤差 (standard Error: Std. Errとも書く) 

• 母標準偏差 \({\sigma}\) の値がわからない時には、標本平均の標準偏差 \({SD (\bar{x}})\)の推定値として標準誤差 \(SE\) を使う
  

• \(t\) 値は、標本平均のバラツキを表す標準誤差 \(SE\) を使って、次の式で求めることができる 

\[t = \frac{\bar{x} - μ}{SE} = \frac{\bar{x} - μ}{u_x / \sqrt{n}} \]
・ \(\bar{x}\) : 標本平均
・ \(μ\) : 母平均
・ \(n\) : 標本サイズ
・ \(u_x\) : 不偏標準偏差（= 標本標準偏差）
・ \(SE\) : 標準誤差 (standard Error)

参考

• 母集団の標準偏差 \({\sigma}\) がわかっている場合 
• \(t\) 分布ではなく標準正規分布 \(z\) 分布)を利用できる 
• \(z\) 値は、次の式で求めることができる 

\[Z = \frac{\bar{x} - μ_0}{σ / \sqrt{n}} \]

• \(z\) 値は、\({\sigma}\) を使って\(\bar{x}\) を変形（＝標準化）した値

• しかし現実のデータ分析で母集団の標準偏差 \({\sigma}\) が既知であることはほとんどない 　　

• 実際には \(t\) 値を使って分析する 

• 自由度 (df: degree of freedom) の大きさごとに 4 種類の \(t\) 分布図を描いてみる

x <- seq(-4, 4, length=100)
hx <- dnorm(x)

df <- c(1, 3, 8, 30)
colors <- c("blue", "red", "green", "gold", "black")
labels <- c("df=1", "df=3", "df=8", "df=30", "正規分布")

plot(x, hx, type="l", lty=2, xlab="t value",
  ylab="Density", main="t Distribution and degree of freedom")

for (i in 1:4){
  lines(x, dt(x,df[i]), lwd=2, col=colors[i])
}

legend("topright", inset=.05, title="Degree of Freedom",
  labels, lwd=2, lty=c(1, 1, 1, 1, 2), col=colors)

• \(t\) 分布の特徴は次のとおり 

1. 0 を中心とした左右対称な分布 
2. 自由度 (degree of freedom: df) の数値によって分布の形が変わる 
3. 自由度が小さい（サンプルサイズが小さいと）\(t\) 値が大きく変動する  → \(t\) 分布表に自由度が小さい時の\(t\) 値を詳しく載せている理由 
   
4. 自由度が大きくなるにつれて（サンプルサイズが大きくなるにつれて）分布が中央に集まる 
   
5. \(t\) 値を求める式の分子は \(({\bar{x} - μ})\) なので、\(({\bar{x} = μ})\) の時（つまり、標本平均が母平均と同じ時） \(t\) 値は 0 になる 
6. 標本平均の推定では、自由度が \(({n - 1})\) の \(t\) 分布表を使う 

• 例えば、標本サイズ n = 100 なら、自由度 99 の \(t\) 分布を使う 
• 「標本平均 \(\bar{x}\) を変形した \(t\) が自由度 99 の \(t\) 分布に従う」ことを利用して統計的推定を行う 
• \(t\) 分布図の黒い実線（自由度 = 99）が示していること
  → 標本を何度も取り直した時、自由度 99 の分布に従う \(t\) の分布の状態 → 確率密度曲線 (probability density curve)

自由度99の \(t\) 分布の確率密度曲線

確率密度曲線と横軸（確率密度 = 0）に囲まれる部分の面積が確率を表す

• 自由度 \(({n - 1 = 99})\) のとき、\(t\) の95％は \(({-1.98})\) 以上、\(({1.98})\) 以下の範囲に収まる（この範囲を\(({[-1.98, 1.98]})\) と表す）
  

• \(t\) 分布の形は自由度によって変わる 
  → 自由度が変われば、\(({-1.98})\) や \(({1.98})\) 以外の値を使う 

• 一般的に、\(t\) が自由度 \(({n - 1})\) の自由度 \(t\) 分布に従うとき
  → \(t\) の 95% は区間を\(({[-t_{n-1,0.025}, t_{n-1,0.025}]})\) に収まる
  → 曲線左右両端の空白部分の面積はそれぞれ \(({0.025 (0.025*2=0.05: 有意水準（α）= 0.05) })\)
  → \(t\) の 90% は区間を\(({[-t_{n-1,0.05}, t_{n-1,0.05}]})\) に収まる
  → 曲線左右両端の空白部分の面積はそれぞれ \(({0.05 (0.05*2=0.10: 有意水準（α）= 0.10) })\)

• 有意水準 (\(α\): level of significance)
  どれくらい低い確率であれば「現実には起こりえない」と言えるのかという基準  

• \(α = 0.05 (5 パーセント)\)を使うのが一般的 

• qt( ) 関数を使うと、有意水準を指定して棄却域を求めることができる
  

• 自由度が 99 で、有意水準 5%、両側検定の時の棄却域は、

qt(0.025, 99)  # t分布で下側確率が0.025、自由度99のt値

[1] -1.984217

qt(0.975, 99)  # t分布で下側確率が0.975、自由度99のt値

[1] 1.984217

5.2 Excercise

Q1: 自由度が 99 で、有意水準 10%、両側検定の時の棄却域を R を使って求めなさい 

Q2: 自由度が 99 で、有意水準 1%、両側検定の時の棄却域を R を使って求めなさい 

6. 信頼区間

6.1 信頼区間の計算方法

• 日本の有権者の民主党に対する感情温度を調査 
• 単純無作為抽出によって n = 100人の標本を入手 
• 不偏標準偏差 \(u_x\) = 21.36
  
• 民主党に対する感情温度: 平均感情温度 \(\bar{x}\) = 55.13 度 
• 有権者全体における母平均感情温度は 50 度だといえるのか？

標本平均 \(\bar{x}\)、 標準誤差 \((SE = \frac{u_x}{\sqrt{n}})\) の時、母平均 \((μ)\) の 95％信頼区間は次の式で求められる 

\(({[\bar{x}-t_{n-1,0.025} - SE, \bar{x} + t_{n-1,0.025} - SE]})\)

• 平均感情温度 \(\bar{x}\) = 55.13度
  
• \(({t_{99,0.025}})\) = 1.98

\[SE = \frac{u_x}{\sqrt{n}}= \frac{21.36}{\sqrt{100}} = 2.136\]

• 以上を代入すると、母平均 \((μ)\) の 95％信頼区間は

\(({[\bar{x}-t_{n-1,0.025} - SE < μ < \bar{x} + t_{n-1,0.025} - SE]})\)

55.13 - 1.98 - 2.136 ≦ \((μ)\) ≦ 55.13 + 1.98 - 2.136 :::

50.9 ≦ \((μ)\) ≦ 59.4

結論 ・有権者全体の民主党に対する感情温度の 95% 信頼区間は 50.9度から 59.4度
→　「50度」は感情温度の 95% 信頼区間に含まれない  
→　母集団での感情温度の平均は 50度だとはいえない

信頼区間のまとめ

1. 複数の標本を抽出すると、標本平均 \(\bar{x}\) を変形した \(t\) が自由度 \(({n - 1})\) の \(t\) 分布に従う
   

\[t = \frac{\bar{x} - μ}{SE} = \frac{\bar{x} - μ}{u_x / \sqrt{n}} \]

2. \(t\) 分布の特性により、すべての標本の 95%について次の式が成り立つ 

\[-t_{n-1,0.025} ≦ \frac{\bar{x} - μ}{SE} ≦ t_{n-1,0.025}\]

3. \(μ\) を中心に変形すると、

\(({[\bar{x}-t_{n-1,0.025} - SE < μ < \bar{x} + t_{n-1,0.025} - SE]})\)

　　

4. 得られた標本数全体の 95% については \(({[\bar{x}-t_{n-1,0.025} - SE]})\) と \(({[\bar{x}+t_{n-1,0.025} - SE]})\) の区間に、真の母平均 \((μ)\) の値が含まれるということ 

5. こうして求められた区間を「95％信頼区間」と呼ぶ 

6.2 信頼区間の解釈

• 「95％信頼区間」の意味：
  

• 「母数の値がこの区間に含まれている確率が95%」は間違い

• それぞれの標本から得られた95%信頼区間に母数が含まれている確率は 0% か 100%のどちらか

  → 上から1 番目のサンプルの95%信頼区間に母数\((μ = 45)\) が含まれる確率・・・100%

  → 上から11番目のサンプルの95%信頼区間に母数\((μ = 45)\)が含まれる確率・・・ 0%

• 標本サイズ100で、標本数20個をとり、20個の「95％信頼区間」が得られたとする
  

• そのうちの95%（つまり19個）は母数をとらえるが、5%（つまり1個）は母数を捉え損ねる

• R を使って感情温度の問題を確かめてみる
  

• 日本の有権者の民主党に対する感情温度を調査 

• 単純無作為抽出によって n = 100人の標本を入手 

• 不偏標準偏差 \(u_x\) = 21.36
  

• 民主党に対する感情温度: 平均感情温度 \(\bar{x}\) = 55.13 度

• 有権者全体における母平均感情温度は 50 度だといえるのか？
  ・母集団における平均感情温度を 55 度、標準偏差を 20 と想定し、100人の標本をランダムに抽出し、記述統計を表示 

set.seed(2017)
emotion <- rnorm(100, mean=55, sd=20)
summary(emotion)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.14   39.56   55.20   55.13   69.00  117.75 

hist(emotion)

sd(emotion)

[1] 21.36289

• R を使って母平均 = 50 という条件で t 検定する 

t.test(emotion, mu=50)

    One Sample t-test

data:  emotion
t = 2.4016, df = 99, p-value = 0.01819
alternative hypothesis: true mean is not equal to 50
95 percent confidence interval:
 50.89171 59.36943
sample estimates:
mean of x 
 55.13057 

• 下から4 〜 5 行目の結果「95 percent confidence interval: 50.89171 − 59.36943」に注目
  →これが「感情温度の 95% 信頼区間」

結論 ・「有権者全体の民主党に対する感情温度の 95% 信頼区間は 50.9度から 59.4度」
→「50度」は感情温度の 95% 信頼区間に含まれない
→ 母集団での感情温度の平均は 50度だとはいえない

• p-値を使った検定の解説は「11. t 検定（1グループ）」を参照

6.3 Exercise

Q1: R を使って、母平均 50、母標準偏差 10 の母集団から sample size 5 の sample を 3 つ (number of samples = 3) 無作為抽出し、3つの sample それぞれのヒストグラム、平均、標準偏差を示しなさい

Q2: Rを使って、母平均 50、母標準偏差 10 の母集団からsample size500 の sample を 3 つ (number of samples = 3`) 無作為抽出し、3つの sample それぞれのヒストグラム、平均、標準偏差を示しなさい

Q3: 上の二種類のシミュレーションからどのようなことが言えるか 

7. 確率変数・確率分布・確率密度関数

7.1 基本事項の確認

• コイン投げ試行のように、ある確率である事象が起こるか起こらないかが定まっている試行をベルヌーイ試行という　　
• コインの表がでる確率を小文字の \((p)\) で表すと、ベルヌーイ試行は次の式で表せる
  

\[P(X = Head) = p\] \[P(X = Tail) = 1 - p\]

• \((P (X))\) は確率を表し 0 から 1 の範囲の実数をとる 
• コインの表裏を表す変数\((X)\) を確率変数と呼ぶ 
• コインの表が出るか裏が出るか、コインを投げ終わらないと確定しない  → このような変数を「確率変数」という
  
• 確率変数のとる各値（この場合は “Head”" と “Tail”）に対応して、それらが起こる確率が与えられる時、その対応を確率分布と呼ぶ 
• 確率変数 \((X)\) とは「確率分布に従って特定の値になる変数」
• 確率変数 \((X)\) は、確率 \((p)\) で表（“Head”）になり、確率 \((1 - p)\) で裏（“Tail”）になる
  
• ベルヌーイ試行を繰り返す 
• コインの表(“Head”)が出る回数の分布を二項分布という 
• R は二項分布の確率密度関数 dbinom( ) を備えている
  
• 10回コインを投げるとき、表が 5 回でる二項分布の確率密度（＝確率）は、

dbinom(5, 10, 0.5)  # dbinom(x, size, probability)

[1] 0.2460938

• 10回コインを投げたとき、表が出る回数が 0 から 10 それぞれの確率を求めると、

dbinom(0:10, 10, 0.5)

 [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
 [6] 0.2460937500 0.2050781250 0.1171875000 0.0439453125 0.0097656250
[11] 0.0009765625

• [1] のすぐ右隣に示されている数字がコインが 0 回出る確率、その右隣の数字が 1 回出る確率と続き、最後の[11] のすぐ右隣に示されている数字が 10 回出る確率 
• 予想通り、5 が出る確率が最も大きく (0.246)、4 と 6 がほぼ同じ程度 (0.205)、3 と 7 がほぼ同じ程度 (0.117) 
• 5 から左右に離れるにつれて確率が小さくなっていることがわかる 
• コインの表が出る確率 0.5 でベルヌーイ試行を 10 回繰り返すと、表が 5 回出る確率が最も高いことがわかる
  
• 10回コインを投げたとき、表が出る「回数」ではなく、表が出る回数の「確率」をプロットすると、

# コインが出る回数を 0 から 10 まで指定、試行は 10 回、表が出る確率は 0.5と指定
coin10 <- dbinom(0:10, 10, 0.5) 
plot(0:10, coin10, type = "h", lwd = 5) # "h"は図に垂直な線を指定、lwd では線の太さを指定

• 縦棒の高さが確率密度（＝確率） 
• （連続量の場合には確率密度＝確率ではなく、区間の面積＝確率）
• ここに示されているデータの出現確率を確率分布という 
• その確率分布を計算する式を確率密度関数 (probability density function) といい、次の式で表す
  

\[f(X) = _NC_XP^X{(1 - P)^{N-X}}\]

• この式に試行回数 \((N = 10)\)、成功回数 \((X = 5)\) を代入してみる 
  
• 10回サイコロを投げて 5 回表が出る組み合わせ \((_NC_X)\) は choose(　) コードを使って計算できる 

choose(10, 5)  # 10回中 5 回表が出る組み合わせ

[1] 252

• \((P^X{(1 - P)^{N-X}})\)の値は、

(0.5)^5*(1-0.5)^5

[1] 0.0009765625

• 従って、\((f(X))\)の値は、

(252)*(0.0009765625)

[1] 0.2460938

• この値は上記ヒストグラムの縦棒の高さ（確率密度＝確率） 

• この確率密度（＝確率）は R では次のように求めることができる 

dbinom(5, 10, 0.5)  # dbinom(x, size, probability)

[1] 0.2460938

7.2 コイントスのシミュレーション

• 前節のシミュレーションでは、コインを 10 回投げて表が出た回数（0 から 10 までの整数）を数えた （＝離散値） 

• ここでは連続量を使ったシミュレーションを行う　　

• 連続量のデータの確率は、ある特定の値ではなく、その値を含むある「区間」を取り、その面積で表す 

• 一例として、最小値が 140、最大値が 185 で一様分布している連続変数（早稲田大学の学生の身長）から 10 人の身長を無作為抽出してみる 

• ここで使うコードはrunif(n, a, b)で \((a)\) が最小値、\((b)\) が最大値、そして \((n)\) が抽出数である
  

runif(10, 140, 185)

 [1] 162.8133 170.4660 181.2077 177.7272 161.2677 171.0087 176.0036 156.1114
 [9] 147.5133 150.7378

• runifの r は random（無作為）、unif は uniform（一様分布）を表している 
• これは無作為抽出なので、実行する度に、異なる数値が抽出される 

runif(10, 140, 185)

 [1] 141.0996 142.8775 171.5061 140.8381 167.1321 170.8180 177.8217 150.8456
 [9] 151.2468 175.7331

• 学術論文を仕上げる際、無作為抽出の結果として同一の結果を得たいのために、無作為抽出する際に set.seed( ) を使って seed を設定できる 

set.seed(2016-01-09)
runif(10, 140, 185)

 [1] 179.9503 179.6794 149.5198 162.7430 148.9673 175.3110 171.7291 168.4825
 [9] 154.0057 140.8939

• 同じ値を seed 設定すると、実行する度に、同一の数値が抽出される
  

set.seed(2016-01-09)
runif(10, 140, 185)

 [1] 179.9503 179.6794 149.5198 162.7430 148.9673 175.3110 171.7291 168.4825
 [9] 154.0057 140.8939

• 無作為抽出で得られた10個の数値をヒストグラムにすると

set.seed(2016-01-09)
hist(runif(10, 140, 185))

• 無作為抽出する数値の数を100個に増やして、それらの数値をヒストグラムにすると

set.seed(2016-01-09)
hist(runif(100, 140, 185))

• 無作為抽出する数値の数を1000個に増やして、それらの数値をヒストグラムにすると

set.seed(2016-01-09)
hist(runif(1000, 140, 185))

• 母集団は最小値が 140、最大値が 185 で一様分布している連続変数だから、抽出する標本サイズを増やす程、母集団の形と似てくることがわかる 
• 一様分布だけではなく、R では 次のコマンドを使うことで様々な分布から無作為抽出することが可能 
• よく使われる分布は次のとおりである
  

rnorm()	: 正規分布
rt()	: \(t\) 分布
rchisq()	: \((\chi^2)\)（カイ二乗）分布

• それぞれの分布から無作為抽出する際、それぞれ特有の引数 (arguments) があるため注意が必要 

• 例えば、正規分布では rnorm() のカッコの中に平均値 (mean) と標準偏差 (sd) を、\(t\) 分布と \((\chi^2)\)分布ではカッコの中に自由度 (df) を指定しなければならない 

• rnorm() を使って無作為抽出してみる 

• 正規分布する母集団 N \((\mu)\), \((\sigma^2)\) から n 個の数値を無作為抽出する時：
  → rnorm(n, mean = mu, sd = sigma) に必要な情報を入力 

• 母平均 (mu) 160、母分散 16 (sigma = 4) で正規分布する母集団から n = 100 個 の数値を無作為抽出する場合、次のように入力すればよい
  

rnorm(100, 160, 4)

  [1] 159.3848 159.4990 164.7559 155.6251 153.1917 158.3015 163.9680 153.7831
  [9] 158.1598 160.8625 161.2840 161.1273 164.9580 164.3071 157.1850 152.9570
 [17] 157.9721 165.1655 161.9859 157.9363 159.3190 153.6560 160.8890 160.0436
 [25] 152.7490 156.1923 166.4202 159.1252 161.6708 158.8871 158.8819 163.8937
 [33] 154.8815 163.1115 159.2016 152.5330 157.0690 152.9178 158.2415 162.3454
 [41] 160.8212 160.1693 162.7420 157.0018 162.6093 156.4365 163.4589 157.0134
 [49] 161.0999 163.0928 154.6236 167.2826 165.1174 152.7759 156.5711 156.8450
 [57] 160.3522 162.3860 159.9502 162.2307 159.0250 163.0924 164.7416 161.5864
 [65] 157.7696 167.6563 159.5414 163.9407 164.6144 156.0849 156.9475 164.7797
 [73] 156.1790 157.2057 157.2734 159.0295 157.1015 157.6920 159.8148 160.7410
 [81] 164.8714 157.5850 160.9043 163.3799 157.6321 156.9417 154.8651 154.0923
 [89] 161.5911 163.2804 165.0042 159.4423 163.6247 161.3810 159.2093 160.0146
 [97] 153.1990 164.3284 150.2262 158.1045

# ヒストグラムを描くと
hist(rnorm(100, 160, 4))

• 母集団が正規分布しているため、標本サイズが大きくなるほど、無作為抽出される標本も正規分布の形に似てくる 

# 5数は
quantile(rnorm(100, 160, 4), c(0, .25, .5, .75, 1))

      0%      25%      50%      75%     100% 
150.0572 157.2044 159.8816 163.0106 168.9840 

# boxplotを描くと
boxplot(rnorm(100, 160, 4))

8. 期待値

8.1 期待値 = 平均値

• 期待値 \((E(X_i))\) とは「ある事柄が起こる確率 \((P(X))\)」と「その結果の数値 \((X_i)\)」を全て掛け合わせ、合計した値のこと　　
• 期待値は「平均値」とも呼ばれる
  
• 数式で表すと、

\[E(X) = \frac{\sum_{i=1}^N p(X_i)X_i}{N}\]

• 事例）
  サイコロを振ったときの期待値は「それぞれの目がでる確率 (1/6)」と「その結果の数値 (1, 2, 3, 4, 5, 6)」を掛け合わせ、合計する
  

\[E(X) = \frac{\sum_{i=1}^N p(X_i)X_i}{N}\]

\[= \frac{(1/6)*1 + (1/6)*2 + (1/6)*3 + (1/6)*4 + (1/6)*5 + (1/6)*6}{N}\]

\[= 3.5\]

• R を使って、サイコロを 10 回振り、その期待値を求めてみると、

x <- 1:6
y <- sample(x, 10, replace = TRUE)
mean(y)

[1] 2.7

• 理論的にはサイコロを振る際の期待値は 3.5 だが、実際にはぴったり 3.5 が得られるわけではない 
• もう 10 回振ってみる 

y <- sample(x, 10, replace = TRUE)
mean(y)

[1] 4.3

• さらにサイコロを振る回数を 10000 回に増やして、期待値を求めてみると 

y <- sample(x, 10000, replace = TRUE)
mean(y)

[1] 3.5204

• サイコロを振る回数を限りなく増やすと、期待値が 3.5 に限りなく近づく（＝大数の法則） 

8.2 期待値のシミュレーション

• 確率分布・・・データの出現状況を確率に置き換えること
  ex) ①サイコロを振る、②コインを投げる、etc
  

• ある確率である事象が起こる時、その事柄を「事象」(event)という

• R には一定の確率で分布する母集団から無作為に標本を抽出する機能が備えられている  【二項分布のシミュレーション】

• R を使ってコインを 10 回投げるシミュレーションをしてみる
  

• まず、コインの表を “H”、裏を “T” と定義して架空のコイン (coin) を作る
  

• コインを 10 回投げて、表と裏が出た回数（整数）をそれぞれ記録する 

• このような整数を離散値と呼ぶ 

coin <- c("H", "T") 
flip <- sample(coin, 10, replace = TRUE) # 10回コインを投げる
table(flip)　# 結果をテーブルに表す

flip
H T 
2 8 

• コインの表 (“H”) がでる確率は理論的には 50 % のはずだが、実際に10 回コインを投げると、ぴったり 5 回表が出るとは限らない 

• コインを 10 回投げる試みを 20 回やってみて、その結果をヒストグラムで表してみる
  

x <- numeric(20) # ベクトル x を作り、そこに 20 個の 0 を入れる
for(i in 1:20) { # コインを 10 回投げる度、変数iを 1 から 20 まで増やす
  flip <- sample(coin, 10, rep = TRUE) # コインを 10 回投げるコマンド
  x[i] <- sum(flip == "H") # 右側の結果"H"を左の x に代入するコマンド
}
table(x) # 結果 x を表にするコマンド

x
4 5 6 7 8 
2 5 9 3 1 

hist(x, breaks = 0:10)　# ヒストグラムの区間を 0 から 10 に限定 

• シミュレーションをやる度に得られる結果は異なるが、10 回のコイン投げを 20 回試みると、必ずしも 5 回が最頻値であるとは限らず、いびつな分布になることが多い 
• 10 回のコイン投げを 10000 回やってみる 

x <- numeric(10000) # ベクトル x を作り、そこに 10000 個の 0 を入れる
for(i in 1:10000) { # コインを 10 回投げる度、変数iを 1 から 10000 まで増やす
  flip <- sample(coin, 10, rep = TRUE) # コインを 10 回投げるコマンド
  x[i] <- sum(flip == "H") # 右側の結果"H"を左の x に代入するコマンド
}
table(x) # 結果 x を表にするコマンド

x
   0    1    2    3    4    5    6    7    8    9   10 
  13  110  427 1213 2021 2429 2072 1139  454  114    8 

hist(x, breaks = 0:10)　# ヒストグラムの区間を 0 から 10 に限定 
abline(v = mean(x), col = "red") # 平均値に赤い縦線を引く

mean(x)

[1] 4.9975

• 結果をヒストグラムで表してみると、平均値の 5 に確実に近づいていることが明らか
  

参考文献

宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」

土井翔平（北海道大学公共政策大学院）「Rで計量政治学入門」

矢内勇生（高知工科大学）統計学2

浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年

浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年

Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.

Kieran Healy, DATA VISUALIZATION, Princeton, 2019

Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017