• 数学B（数研出版）「2章 第1節 確率分布」をRを用いて解説します。
    問題については当該の教科書を参照してください。

確率変数は値を出力する機械であり、その出力値は固有の分布を持つ。
当サイト確率変数と確率分布　では確率変数を連続量から導入している。
一方、数Bの教科書では離散量から導入されている。ここでは教科書に従い離散量から導入する。

当サイトでは「確率変数は自販機に喩えられる」と説明している。
（ 確率変数と自販機 参照）

自販機は値を出力する乱数という関数である。 実際にRで自販機をつくり、シミュレーションしてみよう。

1 例題１（p55）

1.1 母集団

母集団から１つの標本Xを抽出するとき、Xは母集団と同じ分布を持つ確率変数である。
Rでは関数sampleを用いて標本を抽出する。
母集団と抽出による確率変数を構成し、確率変数の出力が母集団分布に一致することを確かめよう。

A=c(rep(10000,5),rep(1000,20),rep(200,75),rep(0,900))#母集団
table(A)#度数表示
## A
##     0   200  1000 10000 
##   900    75    20     5
prop.table(table(A))#確率表示
## A
##     0   200  1000 10000 
## 0.900 0.075 0.020 0.005
X <- sample(A, size = 1, replace = TRUE)#標本

Xは1000本のクジから1本を引くときの出力値であり、このXが確率変数の働きをする。 prop.table を用いると確率表示され、これが母集団の確率分布となる。

1.2 標本

自販機から多くのデータを出力し集計すれば、母集団の確率分布に近くなる。

1.2.1 400人がクジを引く

sizeを400とする。

X400 <- sample(A, size = 400, replace = TRUE)#標本
prop.table(table(X400))
## X400
##      0    200   1000  10000 
## 0.9050 0.0625 0.0200 0.0125
barplot(prop.table(table(X400)))

plot(prop.table(table(X400)))

plotを用いると、横軸は間隔尺度になる。（barplotの横軸は順序尺度）

［参考サイト］

1.2.2 4000人がクジを引く

sizeを4000とする。

X4000 <- sample(A, size = 4000, replace = TRUE)#標本
round(prop.table(table(X4000)),3)
## X4000
##     0   200  1000 10000 
## 0.907 0.072 0.016 0.006
barplot(prop.table(table(X4000)))

plot(prop.table(table(X4000)))

大標本でも完全に母集団分布と一致することはないことに注意。

2 例題２（p57）

2.1 母集団

A=c(rep(1,3),rep(2,6),rep(3,1))#母集団
A
##  [1] 1 1 1 2 2 2 2 2 2 3
table(A)#度数表示
## A
## 1 2 3 
## 3 6 1
prop.table(table(A))#確率表示
## A
##   1   2   3 
## 0.3 0.6 0.1
X <- sample(A, size = 1, replace = TRUE)#標本

2.2 標本

sizeを400とする。

X400 <- sample(A, size = 400, replace = T)#標本
prop.table(table(X400))
## X400
##      1      2      3 
## 0.2950 0.6075 0.0975
mean(X400)
## [1] 1.8025
mean(A)
## [1] 1.8

期待値と400個の平均はほぼ一致している。

2.3 研究

確率分布の表のように、各値の出現頻度を指定して抽出することもできる。

A=c(1,2,3)
X400 <- sample(A, size = 4000, replace = TRUE,prob=c(3,6,1))#標本
prop.table(table(X400))
## X400
##       1       2       3 
## 0.28850 0.60525 0.10625
mean(X400)
## [1] 1.81775

2.4 発展

2.4.1 超幾何分布

Rの超幾何分布乱数を用いて出力してみよう。

X <- rhyper(1000, 3, 2, 3)
prop.table(table(X))
## X
##     1     2     3 
## 0.310 0.585 0.105
hist(X,prob=T,xlim=c(0,4))

\(E(X)=n\frac{M}{N}\)

\(V(X)=n\left(\frac{M}{N}\right)\left(1-\frac{M}{N}\right)\left(\frac{N-n}{N-1}\right)\)
期待値は二項分布と同じ、分散は二項分布の分散に有限修正を施したものになる。

2.4.2 多項分布

結果の出力回数に注目したのが多項分布である。
100回の多項分布の回数分布と超幾何分布をtableで整理したものはほぼ一致している。

#多項分布
d <- rmultinom(n=1, size = 100, prob = c(0.3, 0.6, 0.1))  # 乱数を発生
print(t(d))
##      [,1] [,2] [,3]
## [1,]   33   59    8
#超幾何分布
X <- rhyper(100, 3, 2, 3)
table(X)
## X
##  1  2  3 
## 29 58 13

3 例題３（p64）

3.1 母集団

当りを１、はずれを０で表す。
母集団の度数分布はtableで、確率分布はprop.tableを用いる。

U=c(rep(1,2),rep(0,8))#母集団
U
##  [1] 1 1 0 0 0 0 0 0 0 0
table(U)#度数表示
## U
## 0 1 
## 8 2
prop.table(table(U))#確率表示
## U
##   0   1 
## 0.8 0.2

3.2 非復元抽出

Aさん、Bさんの抽出を非復元にするには、replace = FALSEとする。
抽出を45回行い、当たりくじの数を調べる。
当たりくじの数は、１の数をsumで合計すればよい。

X=c()
Y=c()
for (i in 1:1000){
  s <- sample(U, size = 3, replace = FALSE)
  X <- c(X,sum(s[1:2]))
  Y <- c(Y,sum(s[3]))
}
head(X)
## [1] 0 1 0 1 2 1
head(Y)
## [1] 0 1 0 0 0 0
data01 <- data.frame(X, Y)
head(data01)
X Y
0 0
1 1
0 0
1 0
2 0
1 0 
#t01 <- table(data01$X, data01$Y)
#t01
t01 <- table(data01)
t01
##    Y
## X     0   1
##   0 475 152
##   1 307  38
##   2  28   0
pt01 <- prop.table(t01)
pt01
##    Y
## X       0     1
##   0 0.475 0.152
##   1 0.307 0.038
##   2 0.028 0.000

このようにRを用いた思考実験によるシミュレーションが可能になる。
シミュレーションが威力を発揮するのは近似や収束についてである。

3.3 研究

「2.3_研究」で行ったように、各セルの出現頻度を指定して抽出することもできる。
この場合“dplyr”というパッケージが必要になる。

library(dplyr)
## 
##  次のパッケージを付け加えます: 'dplyr'
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union
data <- data.frame(X=c(0,0,0,1,1,1),Y=c(0,1,2,0,1,2))
data
X Y
0 0
0 1
0 2
1 0
1 1
1 2 
s <- sample_n(tbl = data, size = 1000,replace=T,weight=c(21,14,1,7,2,0))
table(s)
##    Y
## X     0   1   2
##   0 497 299  16
##   1 144  44   0

比率をweightで指定している。
点をずらして散布図をかくと

plot(jitter(X),jitter(Y))

4 例題７（p86）

4.1 二項分布の正規近似（p85）

サイコロを720回投げて１の目が出る回数を調べるという試行を10000セット繰り返し、10000個のデータでヒストグラムを作る。

X <- rbinom(10000, 720, 1/6)
hist(X,prob=T)
curve(dnorm(x,120,10),add=TRUE,col="blue")

理論上、
　　期待値 \(m=\dfrac{720}{6}=120\) 、分散 \(\sigma^2=720\dfrac{1}{6}\cdot \dfrac{5}{6}=100\)、 標準偏差　 \(\sigma=10\)
であり、青のグラフは正規分布 \(N(120,10^2)\) である。

標準化すると標準正規分布にほぼ一致する。

hist(scale(X),prob=T,xlim=c(-5,5))
curve(dnorm(x,0,1),add=TRUE,col="blue")

Z値に変換するにはscale関数を用いる。
青は標準正規分布のグラフである。

720回のうち１が出る回数の確率分布は次のようになる。

x<-0:720
prob<-dbinom(x,720, 1/6)
plot(prob,col=1,type="l")