統計学2

乱数生成と中心極限定理

今回の目標

• 乱数を作ろう！
• 中心極限定理を理解しよう！

準備

必要なパッケージを読み込む。

library(tidyverse)

次に、日本語が正しく表示されるようにする。

## 図のなかで日本語を使えるようにする
## フォントの設定はお好みで
## （Unix/Linux ではIPAexフォントのインストールが必要かも）
if (.Platform$OS.type == "windows") { 
  # Window
  if (require(fontregisterer)) {
    my_font <- "Yu Gothic"
  } else {
    my_font <- "Japan1"
  }
} else if (capabilities("aqua")) {
  # macOS
  my_font <- "HiraginoSans-W3"
} else {
  # Unix/Linux
  my_font <- "IPAexGothic"
}

theme_set(theme_gray(base_size   = 9,
                     base_family = my_font))

乱数の生成

Rを使って乱数を生成してみよう。

コイン投げ

自分で決めた特定の対象から、ランダムにどれかを選びたいときは、sample() 関数を利用する。 たとえば、この関数を使ってコイン投げを実行したいときは、

coin <- c("表", "裏")  # コインを定義する
sample(coin, size = 1)

## [1] "裏"

とする。size で何回選ぶか（何回コインを投げるか）を指定している。

これを何度か実行してみよう。

sample(coin, size = 1)

## [1] "裏"

sample(coin, size = 1)

## [1] "裏"

sample(coin, size = 1)

## [1] "裏"

sample(coin, size = 1)

## [1] "裏"

sample(coin, size = 1)

## [1] "表"

選ぶ回数を変えてみよう。

sample(coin, size = 2)

## [1] "裏" "表"

もう一度やってみよう。

sample(coin)

## [1] "表" "裏"

さらに、もう一度やってみよう。

sample(coin, size = 2)

## [1] "裏" "表"

これを何度やっても、1回目が表なら2回目は裏、1回目が裏なら2回目が表になる。つまり、2回目のコイン投げはランダムではない。これは、sample() が決められた対象から1つずつ順番に選ぶという作業をしているためである。私たちが定義した coin の中身は「表」と「裏」の2つしかない。この2つから順番に選ぶ作業をすると、1つ目に表（裏）が出れば、2回目に残されているのは裏（表）だけなので、2回目がランダムではなくなってしまう。

ためしに選ぶ回数を3回にしてみよう。

sample(coin, size = 3)

## Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'

エラーが出た。選ぶ対象が2つしかないのに3つは選べないのでエラーになる（エラーをよく読むと、そのように書いてある）。

したがって、sample() を使ってコイン投げを2回以上行うには、少し工夫が必要になる。sample() は私たちが特に指示をしないと非復元抽出（sampling without replacement）を行う。非復元抽出というのは、1度選んだものは選択肢から外すという選び方である。何度もコイン投げを繰り返すには、復元抽出 (sampling with replacement) を行えばよい。復元抽出では、1度選んだものも選択肢の中に戻す（選択肢として復元する）という選び方である。 sapmle() で復元抽出を実行するために、replace = TRUE という指示を加える。

sample(coin, size = 5, replace = TRUE)

## [1] "表" "裏" "表" "裏" "裏"

これでコイン投げを複数回実行できる。

試しにコインを10回投げて、その結果を coin10 という名前で保存してみよう。（%>% print()は結果をっ表示するためにつけているだけなので、結果を表示する必要がないならなくてもよい。）

coin10 <- sample(coin, size = 10, replace = TRUE) %>% 
  print()

##  [1] "裏" "表" "裏" "裏" "表" "裏" "表" "表" "裏" "裏"

10回のうち、表は何回出ただろうか？Rを使って数えてみよう。数えるために、以下の方法をとる。

1. 特定のコイン投げ（1回目, 2回目, \(\dots\), 10回目）が表だったかどうか調べる
2. 表の回数を数える

Rである特定の値（数または文字列）になっているか調べたいときは、== （二重等号）を使う（= [等号1つ] は <- と同じで右側の内容を左側に保存してしまうので注意すること）。「二重等号の左側は右側と同じかどうか」を調べ、同じときは TRUE（真）、異なるときは FALSE（偽）という答えが返ってくる。簡単な例で確認してみよう。

a <- 2    # aに2を入れる
a == 2

## [1] TRUE

a == 3

## [1] FALSE

a == "裏" # 文字列と比べるときは文字列を引用符で囲む

## [1] FALSE

これを使って、上で行った10回のコイン投げが表だったかどうか確かめよう。

coin10 == "表"

##  [1] FALSE  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE

このように、各回について、表 (TRUE) か裏 (FALSE) かを教えてくれる。

私たちは表が何回出たかを知りたいので、TRUE の回数を数えればいい。ここでは、10回しか投げていないので、自力で数えることもできるが、数が増えると数えるのは面倒である。そこで、Rを使って数える。RではTRUE とFALSE を数として扱うと、TRUE は1、FALSE = 0 とみなされる。したがって、TRUEの数を数えたいなら、上の結果を合計すればよい。合計は sum() で求められるので、次のようにする。

sum(coin10 == "表")

## [1] 4

表は10回中4回だったことがわかる。

サイコロ投げ

同じ関数を使って、サイコロ (die) 投げを実行してみよう。 まずは、サイコロを定義する。

die <- 1:6  # 1から6までの整数
# 以下のような書き方も可能
# die <- c(1, 2, 3, 4, 5, 6)
# die <- seq(from = 1, to = 6, by = 1)

これを復元抽出すれば、サイコロ投げを何度も行える。100回投げてみよう。

die100 <- sample(die, size = 100, replace = TRUE)

3は何回出ただろうか？

sum(die100 == 3)

## [1] 16

16回出たことがわかる。

ここまでは、「正しいコイン (a fair coin)」や「正しいサイコロ (a fair die)」を想定してきたが、特定の目が出やすいサイコロ（やコイン）を使うこともできる。「1が出る確率だけ他の目の4倍」というサイコロを1,000回投げてみよう。そのために、prob で各目が出る比率を指定する（確率を指定してもよい）。

unfair1000 <- sample(die, 
                     size = 1000, 
                     replace = TRUE,
                     prob = c(4, 1, 1, 1, 1, 1))

1が何回出たか確かめてみよう。

sum(unfair1000 == 1)

## [1] 439

1,000回のうち、439回1の目が出たようだ。

このように、それぞれの選択肢が選べる確率を自由に設定して実験することができる。

また、選ぶ対象も自由に設定できる。たとえば、

sample(c('高知工科大', '高知県立大', '高知大'), 
       size = 1, 
       replace = TRUE, 
       prob = c(3, 2, 1))

## [1] "高知工科大"

のようなこともできる。

実習課題 (1)：

• 1から6までの目がある「正しい」サイコロを2個振るという作業を1,000回繰り返し、出た目の合計が9になる回数を数えてみよう。
• 「正しくない」コイン（表が出る確率が0.5ではないコイン）を500回投げ、表が出た回数を数えよう。表が出る確率は自由に設定してよい（ただし、0, 0.5, 1 を除く）。

確率分布からの乱数生成

Rでは、代表的な確率分布から乱数を生成することが可能である。基本的には、 r（randomの頭文字）と、分布名の最初の数文字を組み合わせた関数を使う。

一様分布

一様分布 (uniform distribution) からの乱数生成には、runif() を使う。この関数では、以下の3つの引数（ひきすう）を指定する。

1. n: 生成する乱数の個数（必ず指定する）
2. min: 最小値（指定しないと 0 に設定される）
3. max: 最大値（指定しないと 1 に設定される）

時計のように1から12（12は0とも考えらえる）の数字が書いてあるような円盤の上でランダムにルーレットを回すことを考える。そうすると、0から12の間のどの位置にルーレットが止まる確率も等しいと考えられる。このような状況を、最小値が0で最大値が12の（連続な）一様分布（連続一様分布; continuous uniform distribution）で表すことができる。このルーレットを100回使ってみよう。

a1 <- runif(n = 100, min = 0, max = 12)

結果をヒストグラムにしてみよう。

df1 <- tibble(a1)  # 結果をデータフレーム (tibble) に入れる
h1 <- ggplot(df1, aes(x = a1)) +
  geom_histogram(binwidth = 1, 
                 boundary = 0,
                 fill = "royalblue", 
                 color = "black") +
  labs(x = "時計盤上のルーレットからの乱数", 
       y = "度数") +
  scale_x_continuous(breaks = 0:12)  # x軸の目盛を調整
plot(h1)

乱数の個数を増やしてみよう。

a2 <- runif(n = 10000, min = 0, max = 12)
df2 <- tibble(a2)
h2 <- ggplot(df2, aes(x = a2)) +
  geom_histogram(binwidth = 1, 
                 boundary = 0,
                 fill = "royalblue", 
                 color = "black") +
  labs(x = "時計盤上のルーレットからの乱数", 
       y = "度数") +
  scale_x_continuous(breaks = 0:12)  # x軸の目盛を調整
plot(h2)

連続ではない一様分布（離散一様分布; discrete uniform distribution）は、上で使った sample() で実行する。「正しいサイコロ」投げは、「1から6までの整数の（非連続な）一様分布」である。

a3 <- sample(1:6, size = 100000, replace = TRUE)
df3 <- tibble(a3)
p3 <- ggplot(df3, aes(x = a3)) +
  geom_bar(width = 0.5) +
  labs(x = "出た目", y = "回数") +
  scale_x_continuous(breaks = 1:6)
plot(p3)

二項分布

コイン投げの結果は、表か裏の2パタン（通常、「成功」と「失敗」と呼ばれる。コイン投げでは、表と裏のどちらを成功と呼んでもよいが、ここでは表を成功、裏を失敗としておく）しかない。 また、同じコインを何度か投げることを繰り返すとき、1回1回のコイン投げで表がでる確率は一定であると考えられる。

このように、結果が成功と失敗の2種類しかなく、成功確率が \(\theta\) で一定（したがって、失敗確率は “\(1-\theta\)” で一定）であるような試行を \(N\) 回繰り返したとき、その成功回数の分布を、「試行回数 \(N\) で成功確率 \(\theta\) の二項分布 (binomial distribution)」と呼ぶ。

理論的には、二項分布の平均値（期待値）は \(N\theta\)、分散は \(N \theta (1 - \theta)\) になる。したがって、標準偏差は \(\sqrt{N \theta (1 - \theta)}\) になる。また、最小値は0（1ではないので注意）、最大値は \(N\) である。

たとえば、「試行回数5で成功確率0.3の二項分布」の平均値は\(5 \cdot 0.3 = 1.5\)、分散は \(5 \cdot 0.3 (1 - 0.3) = 1.05\)、標準偏差は \(\sqrt{1.05} \approx 1.02\) である（「\(\approx\)」は「ほぼ等しい」という意味）。最小値は0、最大値は5である。 この分布は、次の図のような形をしている。

sample() を使わずに、Rでコイン投げを実行してみよう。Rで二項分布に従う乱数を生成する関数は、rbinom() である。この関数で指定しなければならない引数（ひきすう）は以下の3つである。

1. n：実験の回数
2. size: 試行回数（0以上の（非負の）整数）
3. prob: 成功確率（コイン投げで表が出る確率。0以上1以下）

たとえば、rbinom(n = 8, size = 10, prob = 0.4 ) とすると、「表が出る確率が0.4のコインを10回投げる」という実験を8回実行する。

rbinom(n = 8, size = 10, prob = 0.4)

## [1] 5 5 4 7 2 2 5 6

結果として、8つの数字が表示されるが、それぞれの数字が、1回毎の実験（コインを10回投げる）で表が何回出たかを表している。この数は乱数 (random numbers) 、つまり、Rによってランダムに生み出された数字なので、人によって異なる数字が得られるはずで、このWebページと同じ数字が出るとは限らない。 試しにまったく同じ関数をもう一度実行すると、違う数が得られる。

rbinom(n = 8, size = 10, prob = 0.4)

## [1] 4 3 4 4 7 4 2 5

表が出る確率が0.5のコインを使って、1回のコイン投げ実験を1回だけ実行するには、

rbinom(n = 1, size = 1, prob = 0.5)

## [1] 1

とする。0と出れば表が0回出た（つまり、裏が出た）ということであり、1ならば表が出たということである。何度か試してみよう。

rbinom(n = 1, size = 1, prob = 0.5)

## [1] 1

rbinom(n = 1, size = 1, prob = 0.5)

## [1] 1

rbinom(n = 1, size = 1, prob = 0.5)

## [1] 1

rbinom(n = 1, size = 1, prob = 0.5)

## [1] 0

rbinom(n = 1, size = 1, prob = 0.5)

## [1] 0

表が出る確率が0.5のコインを使って、1回のコイン投げ実験を10回まとめて実行するには、次のようにする。

rbinom(n = 10, size = 1, prob = 0.5)

##  [1] 0 1 1 1 1 0 0 0 1 0

実習課題 (2)：

• 3つの引数の値を変えて、様々な条件で二項分布からの乱数生成を試してみよう！

正規分布

正規分布 (normal distribution) からの乱数生成には、rnorm() を使う。この関数では、以下の3つの引数（ひきすう）を指定する。

1. n: 生成する乱数の個数（必ず指定する）
2. mean: 正規分布の平均値（指定しないと 0 に設定される）
3. sd: 正規分布の標準偏差（指定しないと 1 に設定される）

平均と標準偏差を指定しないと、平均が0で標準偏差が1の正規分布、すなわち標準正規分布からの乱数生成が行われる。

標準正規分布から100個乱数を生成してみよう。

b1 <- rnorm(n = 100)

作った乱数をヒストグラムにしてみよう。

df_n1 <- tibble(b1)
h_n1 <- ggplot(df_n1, aes(x = b1)) +
  geom_histogram(color = "black", fill = "dodgerblue") +
  labs(x = "標準正規分布からの乱数", y = "度数")
plot(h_n1)

正規分布に見えるだろうか？ 乱数の平均値と標準偏差を計算してみよう。

mean(b1)

## [1] 0.07602457

sd(b1)

## [1] 0.9794381

元の分布の平均値と標準偏差と比べると、どんなことが言えるだろうか？

乱数の数を増やして同じことしてみよう。

b2 <- rnorm(n = 10000)
df_n2 <- tibble(b2)
h_n2 <- ggplot(df_n2, aes(x = b2)) +
  geom_histogram(color = "black", fill = "dodgerblue") +
  labs(x = "標準正規分布からの乱数", y = "度数")
plot(h_n2)

今度はどうだろうか？

乱数の平均値と標準偏差を計算してみよう。

mean(b2)

## [1] 0.01088606

sd(b2)

## [1] 1.000074

元の分布の平均値と標準偏差と比べると、どんなことが言えるだろうか？

実習課題 (3)：

• 平均が10, 標準偏差が4の正規分布から乱数を1,000個生成し、結果をヒストグラムにしてみよう。また、乱数の平均値と標準偏差を計算し、元になった分布と比べてみよう。

中心極限定理

中心極限定理 (Central Limit Theorem) は、「標本サイズ（サンプルサイズ）が大きくなれば、標本平均は正規分布で近似できる」という定理である（より正確には、「標本サイズを無限大にすると標本平均を標準化したものが標準正規分布に収束する」という定理である）。厳密に理解する（証明する）のはこの授業の範囲を超える（つまり、難しい）。しかし、私たちが統計的推定や検定を行うときに正規分布（特に標準正規分布）ばかり使ってもいいのは、この定理のおかげである。したがって、実際に正規分布ではない分布から正規分布ができるということを理解する必要がある。

そこで、Rでシミュレーションを実行することを通じて中心極限定理を理解しよう。

例として、10個のボールが入った袋を考える。ボールにはそれぞれ0から9までの整数が書いてあるとする。この袋から、ランダムに1つボールを選ぶとすると、選んだボールに書かれた数は、0から9までの整数のどれかで、0から9が選ばれる確率は等しく10分の1（0.1）である。図にすると、以下のようになる。

つまり、ランダムにボールを1つ選ぶというのは、（離散）一様分布からの乱数生成である。 1つボールを選ぶとき、ボールに書かれている数の平均値（期待値）は、\((9 - 0) / 2 = 4.5\) である。

ここで、私たちはボールにどんな数が書かれているか知らない（0から9ではなく、100から109かもしれないし、すべての数が奇数かもしれない）と仮定しよう。この状態で、この袋からボールを引く作業を繰り返し（復元抽出する）、ボールに書かれている数の平均値を当てたい（正解は4.5）。

もっとも単純な方法は、ボールを \(n\) 回引いて、その平均値を使うという方法である。

まず、ボールを2回だけ引いて平均値を当てるという実験を行ってみよう。この実験をすると、1回目のボールの選び方は10通り、2回目のボールの選び方も10通りあるので、全部で100通りの選び方がある。しかし、2つのボールに書かれている数は0から9までの整数なので、可能な合計値は0から18までの19通りであり、平均値は「合計/2」 なので、平均値も19通りしかない。

理論的には、次の図のような確率で、それぞれの平均値が得られる。

この実験を1回だけ行うとき、正解である4.5が選ばれる確率は、0.1 である。 試しに1回やってみよう。

bag <- 0:9  # 袋の中身を定義する
exp_2 <- sample(bag, size = 2, replace = TRUE)  # 復元抽出
mean(exp_2)   # 平均値を求める

## [1] 6.5

今回は、たまたま6.5になった。

では、この実験を1,000回繰り返すと「それぞれの回での平均値の分布」はどんな形になるだろうか。 実際にやってみよう。

Rで同じ作業を繰り返し行う簡単な方法は、for ループを使うことである。for の直後の丸カッコで繰り返し回数を指定し、ループさせる内容を中括弧 { } で囲む。 たとえば、0からスタートして「1を足す」という作業を5回繰り返すには、次のようにする。

まず、スタート時の数である0を保存する。

a <- 0  # a に0を入れる

次に、結果を保存するために、result という名前の入れ物を用意する。NA は「空」の状態を表す（これを欠測値と呼ぶ）。また、 length.out で保存場所が何個分必要か指定する。

result <- rep(NA, length.out = 5) 

この時点で result の中身を確認してみよう。

result  # result の中身を確認する

## [1] NA NA NA NA NA

すべて NA になっている。

forループを利用して、数を1だけ加える作業を5回繰り返す。

for (i in 1:5) {    # i が1から5までの繰り返し
  a <- a + 1        # a に1を足す
  result[i] <- a    # i番目の足し算の結果を result[i] に入れる    
}
result  # resultの中身を確認

## [1] 1 2 3 4 5

このループを利用して、ボールを2回引く作業を1,000回繰り返してみよう。

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

結果をヒストグラムにしてみよう。

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 = "度数") +
  scale_x_continuous(breaks = 0:9)
plot(h_sim1)

この分布は正規分布に見えるだろうか？

次に、サンプルサイズ \(N\) を2から5に増やし、同様の実験をしてみよう。

N <- 5                                # サンプルサイズ
trials <- 1000                        # 実験の繰り返し回数
sim2 <- rep(NA, length.out = trials)  # 結果の保存容器
for (i in 1:trials) {
  experiment <- sample(bag, size = N, replace = TRUE)  # 復元抽出
  sim2[i] <- mean(experiment)         # i 回目の平均値を保存
}

結果をヒストグラムにしてみよう。

df_sim2 <- tibble(avg = sim2)
h_sim2 <- ggplot(df_sim2, aes(x = avg)) +
  geom_histogram(binwidth = 0.5, 
                 boundary = 0.5,
                 color = "black") +
  labs(x = "5個のボールの平均値", 
       y = "度数") +
  scale_x_continuous(breaks = 0:9)
plot(h_sim2)

この分布は正規分布に見えるだろうか？

サンプルサイズ \(N\) を10に増やしてみよう。

N <- 10                              # サンプルサイズ
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.5,
                 color = "black") +
  labs(x = "10個のボールの平均値", y = "度数")
plot(h_sim3)

この分布は正規分布に見えるだろうか？

サンプルサイズ \(N\) を100に増やしてみよう。

N <- 100                             # サンプルサイズ
trials <- 1000                       # 実験の繰り返し回数
sim4 <- rep(NA, length.out = trials) # 結果の保存容器
for (i in 1:trials) {
  experiment <- sample(bag, size = N, replace = TRUE)  # 復元抽出
  sim4[i] <- mean(experiment)        # i 回目の平均値を保存
}

結果をヒストグラムにしてみよう。

df_sim4 <- tibble(avg = sim4)
h_sim4 <- ggplot(df_sim4, aes(x = avg)) +
  geom_histogram(binwidth = 0.125, 
                 color = "black") +
  labs(x = "100個のボールの平均値", y = "度数")
plot(h_sim4)

このように、元の分布は一様分布でも、サンプルサイズ \(N\) を増やすと、「平均値の分布」は正規分布に近づく。よって、サンプルサイズ \(N\) が十分大きい（大まかな目安は \(N > 100\)）とき、正規分布を使って統計的推定や検定を行うことが許される。

実習課題 (4)：

• \(N=5\) と \(N=500\) のそれぞれの場合について、同様の実験を行ってみよう。

授業の内容に戻る