7 母集団と標本

佐久間 智広
sakuma@port.kobe-u.ac.jp

神戸大学大学院経営学研究科

2024/05/27


0.1 準備

Code 
rm(list=ls()); gc();  gc(); 
1if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, magrittr,
2               estimatr ,modelsummary)
1
複数のパッケージを一度に読み込めるpacmanパッケージが入っていない場合は，インストールする
2
パッケージの読み込み


pacman

複数のパッケージを読み込めるパッケージ。そのパッケージの中のp_loadコマンドを使うと，（）内に指示したパッケージを読み込める。さらにそのパッケージがそもそもインストールされていない場合にはインストールした後に読み込んでくれる

tidyverse

Rのコードを楽にするたくさんのパッケージをパッケージにしたもの

magrittr

tidyverseに含まれる%>%などのコードをさらに拡張するもの。これがあると%$%とか%<>%が使える

1 概要

データを分析する際に，知りたい人やもののすべてのデータを使うことは稀です。

代わりに，知りたい集団の一部のデータを使って全体を予測します。

これがなぜ達成可能なのかについてその背後の論理を説明します。

例

メディアでは，選挙の日に出口調査を行い，開票完了前に「当選確実」を出します。出口調査はすべての選挙会場（総務省によると全国47,000箇所以上あるそうです）に調査員を配置して，投票したすべての有権者（5,000万人以上）に投票内容を聞くのは現実的ではありません。全国の限られた数の投票所の限られた数の有権者からの意見でなぜ「当選確実」を出せるのでしょうか？

2 母集団と標本

  • 社会調査や学術研究では，知りたい対象と手元のデータが必ずしも一致していません。そもそも，知りたい対象全体のデータを取るような調査はほとんど行われません。
    • 唯一の例外は国勢調査
  1. 例えば，世論調査は日本在住者の政治に対する意見などを知るためにデータを集めますが，全国民を対象としているわけではありません。
    • 全国民の意見を聞いていたら時間もお金もかかります。そもそも全員と連絡をつけることは非現実的です。
    • そこで，例えば1,000人に聞いた意見から日本全体の世論を推測します。
  2. 例えば，新しい薬の効果を検証したい場合，知りたい対象は，過去・現在・未来に渡る特定の病気の患者への効果です。調査時点で過去の患者や未来の患者に対する効果を知ることはできません。
    • そこで，治験対象者のデータから，人類一般に対する効果を推測します。
  • ここで，知りたい対象を母集団，検証に使う母集団の一部を標本と言います。また，標本に含まれる個体の数を標本サイズ（サンプルサイズ）と言います。

特に経済学などでは，上記1つ目の例のように日本人全体，といった特定の母集団を想定しないことが多いとされます(星野ほか 2023; 高橋 2022)。その代わりに，サンプルを生み出す特定の確率分布を母集団と呼びます。

データ解析の目的は，母集団がどんな形になっているのか，標本を発生させた母集団はそもそもどんな形をしているのかを知ることにあります。

  • 日本全体の内閣支持率は何％なのか？
  • ラーメンが好きな日本人の割合は？

2.1 サンプリング（無作為抽出）

  • 確率論を基礎として，十分な数の適切に選ばれた標本を適切な方法で分析すると，非常に高い精度で母集団の傾向を推測できることが知られています。
    • 逆に，うまくサンプルを取らなかったら当然ズレます

例

東京大学の学生の1,000人のTOEICの点数の平均を使って，日本国民の英語力を推定する。

  • この場合，母集団は東京大学の学生？

適切にデータを選ぶ，という際に重要なのが，標本が，母集団から偏りなく抽出されているということです。

  • 母集団に含まれる各個体が一様に等しい確率で選ばれているということ
  • 無作為抽出もしくはランダムサンプリングという
  • 無作為抽出によって選ばれた標本を無作為標本という

2.2 確率モデル

データ解析の目的は，無作為標本（母集団からすると限られた情報）を使って，母集団全体の傾向を推測することにあります。

その際，母集団をモデル（様々な仮定を置いた数式）で表します。

例

  • 野球の打席→ヒットor凡退
  • ヒットを1，凡退を0，ヒットを打つ確率をpとして，n回打席に立ってk回ヒットが出る確率を数式で表す（モデルで記述する）と

\[ P(X=k) = {}_n C_k p^k(1-p)^{n-k} \]

標本の数値と確率モデルを利用して，母集団の分布の形を把握したい

確率モデル（確率分布）にはいくつかの便利な性質があり，それを利用することで未知の母集団の分布の形が推測できる

3 確率分布

母集団を数学的に表現する確率分布には，様々な種類があります。

  • ベルヌーイ分布
  • 二項分布
  • 多項分布
  • ポアソン分布
  • 正規分布
  • t分布
  • \(\chi^2\)分布
  • F分布
  • …

後述する理由で，一般的な分析はほとんど正規分布とそれを元にできる分布（t, \(\chi^2\), Fなど）を使います。

3.1 二項分布

\[ P(X=k) = {}_n C_k p^k(1-p)^{n-k} \]

打率3割5分の人と1割5分の人がそれぞれ10回打席に立ってヒットが何本出るか？

Code 
# パラメータを指定
phi <- 0.35

# 試行回数を指定
M <- 10

x_vals <- 0:M

prob_df <- tidyr::tibble(
  x = x_vals, # 確率変数
  probability = dbinom(x = x_vals, size = M, prob = phi) # 確率
)

ggplot(data = prob_df, mapping = aes(x = x, y = probability)) + # データ
  geom_bar(stat = "identity", position = "dodge", fill = "#00A968") + # 棒グラフ
  scale_x_continuous(breaks = x_vals, labels = x_vals) + # x軸目盛
  labs(title = "Binomial Distribution", 
       subtitle = paste0("phi=", phi, ", M=", M), # (文字列表記用)
       #subtitle = parse(text = paste0("list(phi==", phi, ", M==", M, ")")), # (数式表記用)
       x = "x", y = "probability") # ラベル

Code 
# パラメータを指定
phi <- 0.15

# 試行回数を指定
M <- 10

x_vals <- 0:M

prob_df <- tibble(
  x = x_vals, # 確率変数
  probability = dbinom(x = x_vals, 
                       size = M, 
                       prob = phi) # 確率
)

ggplot(data = prob_df, 
       mapping = aes(x = x, 
                     y = probability)) + # データ
  geom_bar(stat = "identity", 
           position = "dodge", 
           fill = "#00A968") + # 棒グラフ
  scale_x_continuous(breaks = x_vals, 
                     labels = x_vals) + # x軸目盛
  labs(title = "Binomial Distribution", 
       subtitle = paste0("phi=", phi, ", M=", M), # (文字列表記用)
       #subtitle = parse(text = paste0("list(phi==", phi, ", M==", M, ")")), # (数式表記用)
       x = "x", y = "probability") # ラベル

3.2 正規分布

身長などの身体測定値や測定誤差など様々な場面に現れる分布。あとで使う望ましい性質がある。平均 \(\mu\)，分散\(\sigma^2\) として，\(\mathcal{N}(\mu,\sigma^2)\) と標記

\[ f(x)= \frac{1}{\sqrt {2 \pi \sigma^2 }} \exp \left(- \frac{(x- \mu)^2}{2\sigma^2} \right) \]

特に平均0、分散1の正規分布を標準正規分布と呼ぶ

  • 左右対称な釣鐘型
  • 頂点が平均（\(\mu\)）
  • その点から標準偏差（\(\sigma\)）の大きさごとに確率が決まっている
    • 1シグマ範囲 \(\mu - \sigma \le X \le \mu + \sigma\) 確率68.3%
    • 2シグマ範囲 \(\mu - 2\sigma \le X \le \mu + 2\sigma\) 確率95.4%
    • 3シグマ範囲 \(\mu - 3\sigma \le X \le \mu + 3\sigma\) 確率99.7%


例

例えば，小学校1年生の児童の身長が，平均\(\mu=118\), 分散\(\sigma^2=6^2\)の正規分布に従っているとすると, \(118 \pm 6 cm\)の範囲に68.3%， \(118 \pm 12 cm\)の範囲に95.4%の児童が入ることが予測できる。



例2

あるテストの結果が平均点が55点，標準偏差15の正規分布に従っているとすると，\(55 \pm 15\)つまり偏差値40-60の範囲内に68.3%が入る。55+15 = 70点の人は偏差値60で，上位15.85％あたり( %6%)。85点の人は，偏差値70で，上位2.3%あたり。

3.3 確率分布の形なんて…

  • 正規分布だとわかっていれば，どの範囲にどのぐらいの数字が現れるのかがわかる
  • でも，母集団の分布がわからへんから困っているわけであって，\(\mu\)も\(\sigma\)も何分布かもわかってない

サンプリングをうまくやって，確率分布の望ましい性質を使って解決する！

4 標本平均と不偏分散

4.1 (復習+α)

  • 平均はデータの特徴を表す代表的な値の一つでした。変数 \(X\) の平均 \(\bar X\)は，

\[ \bar X = (X_1 + X_2 + ... X_n )/n = \frac{1}{n} \sum_{i=1}^{n} X_i \]

  • 分散はデータのばらつきを表します。データは平均より高かったり低かったりするので，各データの平均からのずれ \((X_i - \bar X)\)を2乗して

\[ \begin{split} s_X^2 &= \frac{1}{n-1}\{(X_1 - \bar X)^2+(X_2 - \bar X)^2+(X_3 - \bar X)^2+...(X_n - \bar X)^2 \} \\ &= \frac{1}{n-1}\Sigma (X_i - \bar X)^2 \end{split} \]

前回と違って，\(n-1\)で割っています。前回のものを標本分散，今回のものを不偏分散と言います。標本を使って母集団の分散を推定する際には，不偏分散の方が良い（より正確になる）ためです。

なんで\(n-1\)で割るのかは，この授業の範囲を超えるので省略します。例えばこちらが参考になります。

4.2 標本平均と不偏分散の性質

  • ある母集団から得た標本があります。母集団全体の平均（母平均）を \(\mu\) とします。この \(\mu\) は未知です。

  • ここで，標本観測と母平均との関係を以下のように想定します。

\[ X_i = \mu+u_i, u_i  \sim N(0,\sigma ^2) \]

　ただし\(u_1,u_2,u_3...u_n\)は互いに独立(前の数字に後からでてくる数字が影響を受けない)

このとき，以下の性質が知られています。

\[E(\bar X) = \mu\]
  • 標本平均の期待値は，母平均となる（不偏性）

標本平均の不偏性は，どんな分布でも成り立つことがわかっています。

要するに…うまくサンプルをとることができれば，標本平均を使って母平均を推測できるということです。

  • 分散についても同じで，不偏分散の期待値は，母集団の分散（母分散）になります。

要するに…うまくサンプルをとることができれば，不偏分散を使って母分散を推測できるということです。

4.3 大数の法則

大数の法則

母集団から無作為に抽出された標本の平均値は，標本の大きさが大きいほど，母集団平均に近い値を取る法則

\(n\)を大きくすると，標本平均\(\bar X = (X_1 + X_2 + ... X_n )/n = \frac{1}{n} \sum_{i=1}^{n} X_i\) は真の平均に近づく

\[ \bar X \to \mu \]

4.3.1 大数の法則のシミュレーション

サイコロには１から6までの目があり，それぞれが出る確率は1/6です。したがって，母集団の平均は，

\[ \frac{1}{6}\sum^6_{i=1}i=\frac{21}{6}=3.5 \]

サイコロを10回振った時

Code 
1s10 <- sample(1:6, 10, replace = TRUE)
2mean(s10)
1
sample()は，決められた数の標本を無作為に選び出す関数です。ここでは，1から6までの数を無作為に10個選び出しています。replaceは復元抽出をするかしないかを指定するオプションです。「１から6までの数を無作為に10個復元抽出で選び取る」という命令をしていることになり，これをs10と名づけて保存しています。
2
このs10の平均を計算しています。 
[1] 2.8

同じように，1,000回サイコロを振ることを以下のように再現してみます。

Code 
s1000 <- sample(1:6, 1000, replace = TRUE)
mean(s1000)
[1] 3.514

（無作為抽出しているので，結果は個々人で違うと思いますが）おそらく10回の時よりも1,000回の時の方が期待値である3.5に近いはずです。

これをさらに以下のように検討してみます。

Code 
#サイコロを10回振るのを1000回繰り返す
1S <- 1000
2rec10 <- numeric(S)
for(i in 1:S){
  rec10[i] <- sample(1:6, 10, replace = TRUE)|>
    mean()
3}
4summary(rec10)
1
Sという名前で数字の1000を保存しています。
2
numericというコマンドで空（0）の行をSの数だけ作って，rec10と名付けています。
3
1から6までの数字を無作為に10個（サイコロを10回振る）時の平均を計算して，rec10のi行目に入れる，という作業をS回(1000行分)繰り返しています。
4
そうして作られたrec10の最小値，最大値，平均，中央値等を出しています。 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.800   3.100   3.500   3.517   3.900   5.100

これをサイコロを10000回振るのを1000回繰り返すのに変えてみます。

Code 
#サイコロを10000回振るのを1000回繰り返す
rec10000 <- numeric(S)
for(i in 1:S){
  rec10000[i] <- sample(1:6, 10000, replace = TRUE)|>
    mean()
}
summary(rec10000)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.441   3.489   3.500   3.499   3.510   3.546
Code 
hist(rec10000)

3.5にかなり収束しました。

サンプルサイズ100のものも加えてヒストグラムを書きました。一回一回のサンプルサイズが大きいほど，母平均に近い値を取っていることが見えます。

5 標本平均の分布

5.1 正規母集団と標本平均の分布

母集団が正規分布であれば

\(X_1, X_2,X3,\cdot\cdot ,X_n\)は正規分布\(N(\mu,\sigma^2)\)からの無作為標本であるとする。このとき標本平均\(\bar X =(X_1+ X_2+X_3+\cdot\cdot +X_n)/n\)の分布は正規分布\(N(\mu,\sigma^2/n)\)

つまり，母集団が正規分布の時は，標本平均も正規分布に従うということが知られています。

例

20代男性の身長の母集団が正規分布\(N(\mu,5^2)\) として，無作為に選ばれた100人の身長の測定値が\(X_1, X_2,X_3,\cdot\cdot ,X_{100}\)とする。標本平均\(\bar X\)で母集団の平均\(\mu\)を推定したい時の測定誤差\(|\bar X - \mu|\)は，上の定理から正規分布\(N(\bar X,0.5^2)\) に従う。これがわかると，2σ区間が \(\mu - 1 \le \bar X \le \mu + 1\)。正規分布で2σ区間に入る確率は95.4%だったので，95.4%の確率で測定誤差が1cm以下となることがわかります。

5.2 中心極限定理

重要なことはこの「標本平均が正規分布になる」という性質は母集団が正規分布じゃなくてもデータが十分に大きい時には成り立つ，ということです。これを中心極限定理と言います。


\(X_1, X_2,X3,\cdot\cdot ,X_n\)は母平均が\(\mu\)母分散が\(\sigma^2\)からの大きさ\(n\)の無作為標本とする。

\(n\)が大きくなるに従って，標本平均\(\bar X=(X_1+ X_2+X_3+\cdot\cdot +X_n)/n\) の分布は，\(N(\mu,\sigma^2/n)\)に近づく（収束する）


中心極限定理を使うと，母集団が正規分布でなくとも，データが大きくなる時は標本平均が正規分布に従うとして計算できる

5.2.1 シミュレーション

サイコロには１から6までの目があり，それぞれが出る確率は1/6です。したがって，母集団の平均は，

\[ \frac{1}{6}\Sigma^6_{i=1}i=\frac{21}{6}=3.5 \]

これも，標本が大きくなると，平均が正規分布になる？

サイコロを1回振った時，平均＝出た目。これを1000回繰り返す

Code 
S <- 1000
rec1 <- numeric(S)
for(i in 1:S){
  rec1[i] <- sample(1:6, 1, replace = TRUE)|>
    mean()
}
sim1 <- tibble(avg = rec1)
ggplot(sim1, aes(x = avg)) +
  geom_histogram(binwidth = 1,
                 color = "black")

同様にサイコロを2回振って平均を取る

Code 
S <- 1000
rec2 <- numeric(S)
for(i in 1:S){
  rec2[i] <- sample(1:6, 2, replace = TRUE)|>
    mean()
}
sim2 <- tibble(avg = rec2)
ggplot(sim2, aes(x = avg)) +
  geom_histogram(binwidth = 0.5,
                 color = "black")

10回

Code 
S <- 1000
rec10 <- numeric(S)
for(i in 1:S){
  rec10[i] <- sample(1:6, 10, replace = TRUE)|>
    mean()
}
sim10 <- tibble(avg = rec10)
ggplot(sim10, aes(x = avg)) +
  geom_histogram(binwidth = 0.5,
                 color = "black")

100回

Code 
S <- 1000
rec100 <- numeric(S)
for(i in 1:S){
  rec100[i] <- sample(1:6, 100, replace = TRUE)|>
    mean()
}
sim100 <- tibble(avg = rec100)
ggplot(sim100, aes(x = avg)) +
  geom_histogram(binwidth = 0.1,
                 color = "black")

5.2.2 シミュレーション2：イカサマコイン

（高校までの確率では「歪みのないコイン」が使われますが）あるコインは，表が出る確率が35％になるような歪み方をしています。

コインを1回投げる，という作業を1000回繰り返してみます。

Code 
# パラメータを指定
phi <- 0.35

# データ数を指定
N <- 1

S <- 1000
rec1 <- numeric(S)
for(i in 1:S){
  rec1[i] <- rbinom(n = N, size = 1, prob = phi)|>
    mean()
}
sim1 <- tibble(avg = rec1)
ggplot(sim1, aes(x = avg)) +
  geom_histogram(binwidth = 1,
                 color = "black")

コインを2回投げて平均をとります。プログラムのsizeの部分を2に変えています。

Code 
rec2 <- numeric(S) 
for(i in 1:S){
  rec2[i] <- rbinom(n = N, size = 2, prob = phi)|>
    mean()
} 
sim2 <- tibble(avg = rec2)
ggplot(sim2, aes(x = avg)) +
  geom_histogram(binwidth = 1,
                 color = "black")

同様に10回投げた平均

Code 
rec10 <- numeric(S)
for(i in 1:S){
  rec10[i] <- rbinom(n = N, size = 10, prob = phi)|>
    mean()
}
sim10 <- tibble(avg = rec10)
ggplot(sim10, aes(x = avg)) +
  geom_histogram(binwidth = 1,
                 color = "black")

100回

Code 
rec100 <- numeric(S)
for(i in 1:S){
  rec100[i] <- rbinom(n = N, size = 100, prob = phi)|>
    mean()
}
sim100 <- tibble(avg = rec100)
ggplot(sim100, aes(x = avg)) +
  geom_histogram(binwidth = 1,
                 color = "black")

10000回

Code 
rec1000 <- numeric(S) 
for(i in 1:S){
  rec1000[i] <- rbinom(n = N, size = 10000, prob = phi)|>
    mean()
} 
sim1000 <- tibble(avg = rec1000)
ggplot(sim1000, aes(x = avg)) +
  geom_histogram(binwidth = 10,
                 color = "black")

6 まとめ

ここまでをまとめると

  1. 母集団から無作為に標本を集める（無作為抽出）
  2. 母集団がどんな分布でも，標本平均は不偏（不偏性）で，なおかつ正規分布に従う（中心極限定理）
  3. 正規分布に従うとわかると，標本平均がどれぐらいの信頼度で母集団を推定できているかを把握できる

無作為抽出という手順と，確率の性質を組み合わせて利用することで，未知の母集団を推測できる！

  • 次回は，この戦略を使って，母集団の形を推測する手続きについてやります。

7 参考

7.1 大数の法則

  • 一般に，もしくは直感的に，アンケートでもなんでも「データをたくさん集めれば集めるほど正確な結果が得られる」と考えられています。
  • この「サンプルサイズ(n)が大きくなるほど精度が高まる（ばらつきが小さく）」ことは，大数の法則として知られています。
大数の法則

母集団から無作為に抽出された標本の平均値は，標本の大きさが大きいほど，母集団平均に近い値を取る法則。

7.1.1 標本分散の性質

標本分散や標準偏差についても大数の法則が成立します

Code 
S <- 1000
srec10 <- numeric(S)
for(i in 1:S){
  srec10[i] <- sd(rnorm(10,50,10))}
srec100 <- numeric(S)
for(i in 1:S){
  srec100[i] <- sd(rnorm(100,50,10))}
srec10000 <- numeric(S)
for(i in 1:S){
  srec10000[i] <- sd(rnorm(10000,50,10))}

hist(srec10000, 
     col="#FF00007F", 
     xlim=c(1,18),ann=F, main="", xlab="" 
     ,breaks=seq(1,20,1))
par(new=T)
hist(srec10, 
     breaks=seq(1,20,1),
     col="#0000FF7F",ann=F, add=T)
hist(srec100, 
     breaks=seq(1,20,1),
     col="#32CD32",ann=F, add=T)
legend("topright", 
       legend=c("rec10000", "rec10", "rec100"),
       fill=c("#FF00007F", "#0000FF7F","#32CD32"))

7.2 普遍性と一致性

あるサンプルデータを用いて推定を行う時，その推定の良さを考える上で重要な性質が2つあります。

不偏性(Unbiasness)
  • 推定値の期待値が真の値と同じになる
  • サンプルサイズに依存しない
  • 不偏性が担保されないと，サンプルをいくら増やしても真値（母集団のパラメータ）が推定できない
一致性(Consistency) 1
  • 推定値がある値に確率収束する
  • サンプルサイズの大きさが重要
  • サンプルサイズｎが大きくなるに従って、推定量が、母集団の真の値に近づいていく

2024社会調査法（立命館大学）