サイコロの出目の合計の分布

「\(n\)個のサイコロの出目の合計の分布は\(n\)が大きくなると正規分布に近づく」ことを、Rを使って確認しましょう。

先日のゼミで「\(n\)個のサイコロの出目の合計の分布は\(n\)が大きくなると正規分布に近づく」という話をしました。 このことをRを使って確認しましょう。

サイコロの確率分布

サイコロの目は1〜6までの6種類あります。

d1 <- 1:6
d1

[1] 1 2 3 4 5 6

これを集計して棒グラフにしてみます。 データの個数（この場合は6）で割って、確率に直しています。

# n = 1 のとき
(table(d1) / length(d1)) |> barplot(main = "n = 1", col = "lightblue")

サイコロが2個のとき、出目の和の分布はどうなるでしょうか？ outer関数を使って、すべての組みわせについて出目の和を計算してみます。

d2 <- outer(d1, d1, `+`)
d2

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    2    3    4    5    6    7
[2,]    3    4    5    6    7    8
[3,]    4    5    6    7    8    9
[4,]    5    6    7    8    9   10
[5,]    6    7    8    9   10   11
[6,]    7    8    9   10   11   12

\(6\times 6=36\)通りのすべての組み合わせについて出目の和を計算すると、2〜12までの数字になります。 これを集計して棒グラフにすると、下の図のように三角形の分布になることがわかります。

(table(d2) / length(d2)) |> barplot(main = "n = 2", col = "lightblue")

サイコロが3つの場合についても、同じようにやってみましょう。 この場合は、\(6^3=216\)通りの組み合わせについて、出目の合計を計算すればいいですね。

d3 <- outer(d2, d1, `+`)
d3

, , 1

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    3    4    5    6    7    8
[2,]    4    5    6    7    8    9
[3,]    5    6    7    8    9   10
[4,]    6    7    8    9   10   11
[5,]    7    8    9   10   11   12
[6,]    8    9   10   11   12   13

, , 2

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    4    5    6    7    8    9
[2,]    5    6    7    8    9   10
[3,]    6    7    8    9   10   11
[4,]    7    8    9   10   11   12
[5,]    8    9   10   11   12   13
[6,]    9   10   11   12   13   14

, , 3

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    5    6    7    8    9   10
[2,]    6    7    8    9   10   11
[3,]    7    8    9   10   11   12
[4,]    8    9   10   11   12   13
[5,]    9   10   11   12   13   14
[6,]   10   11   12   13   14   15

, , 4

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    6    7    8    9   10   11
[2,]    7    8    9   10   11   12
[3,]    8    9   10   11   12   13
[4,]    9   10   11   12   13   14
[5,]   10   11   12   13   14   15
[6,]   11   12   13   14   15   16

, , 5

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    7    8    9   10   11   12
[2,]    8    9   10   11   12   13
[3,]    9   10   11   12   13   14
[4,]   10   11   12   13   14   15
[5,]   11   12   13   14   15   16
[6,]   12   13   14   15   16   17

, , 6

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    8    9   10   11   12   13
[2,]    9   10   11   12   13   14
[3,]   10   11   12   13   14   15
[4,]   11   12   13   14   15   16
[5,]   12   13   14   15   16   17
[6,]   13   14   15   16   17   18

(table(d3) / length(d3)) |> barplot(main = "n = 3", col = "lightblue")

このように計算していくと、\(n=10\)のときには\(6^{10}\)（およそ6千万）の数字を計算しなくてはならなくなってしまいます。 これは効率的ではないので、以降では畳み込み（convolve関数）を使って計算した結果を示します。

サイコロ1個のとき、すべての出目の確率は等しく\(\frac{1}{6}\)です。

p <- rep(1/6, 6)

棒グラフにしてみます。 先程と同じグラフになっていることを確認してください。

# n = 1 のとき
dist <- p
names(dist) <- 1:6
dist |> barplot(main = "n = 1", col = "lightblue")

サイコロ2個の場合には、先程の分布にサイコロの出目の分布を畳み込みます。

# n = 2 のとき
dist <- convolve(dist, p, type = "open")
names(dist) <- 2:12
barplot(dist, main = "n = 2", col = "lightblue")

サイコロが3個に増えたら、もう一度畳み込みを行います。

# n = 3 のとき
dist <- convolve(dist, p, type = "open")
names(dist) <- 3:18
barplot(dist, main = "n = 3", col = "lightblue")

サイコロが増えただけ、畳み込みを繰り返せばいいです。 10個のサイコロを考えたいので、あと7回繰り返します。

for(i in 4:10)dist <- convolve(dist, p, type = "open")
names(dist) <- 10:60
barplot(dist, main = "n = 10", col = "lightblue")

dnorm関数を使って、この棒グラフに理論的な正規分布のグラフを重ねてみましょう。

1つのサイコロを考えたとき、その出目の平均\(\mu\)は

\[\mu = \frac{1+2+3+4+5+6}{6}= \frac{7}{2}\]、

分散\(\sigma^2\)は

\[\sigma^2 = \frac{(1 - 7/2)^2+(2 - 7/2)^2+(3 - 7/2)^2+(4 - 7/2)^2+(5 - 7/2)^2+(6 - 7/2)^2}{6}=\frac{35}{12}\]

となります。 これを使えば、サイコロが10個の時の平均値と分散も簡単に計算することができます。

x <- barplot(dist, main = "n = 10", col = "lightblue")
mu <- 7 / 2 * 10
sd <- sqrt(35 / 12 * 10)
lines(x, dnorm(10:60, mean = mu, sd = sd), col = "red", lwd = 2)

正規分布にかなり近い分布になっていることが視覚的に明らかになりました。