-数学B（数研出版）「第２章 統計的な推測」を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.8575 0.1275 0.0125 0.0025
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.902 0.071 0.022 0.004
barplot(prop.table(table(X4000)))

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

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

2 例題２（p59）

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 = TRUE)#標本
prop.table(table(X400))
## X400
##      1      2      3 
## 0.3175 0.6075 0.0750
mean(X400)
## [1] 1.7575

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

2.3 超幾何分布

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

X <- rhyper(1000, 3, 2, 3)
prop.table(table(X))
## X
##     1     2     3 
## 0.294 0.604 0.102
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 多項分布

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

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

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] 1 0 0 1 0 0
head(Y)
## [1] 0 0 0 0 0 0
data01 <- data.frame(X, Y)
head(data01)
X Y
1 0
0 0
0 0
1 0
0 0
0 0 
#t01 <- table(data01$X, data01$Y)
#t01
t01 <- table(data01)
t01
##    Y
## X     0   1
##   0 462 159
##   1 316  45
##   2  18   0
pt01 <- prop.table(t01)
pt01
##    Y
## X       0     1
##   0 0.462 0.159
##   1 0.316 0.045
##   2 0.018 0.000

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

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

X <- rbinom(100, 10, 1/6)
hist(X,prob=T,xlim=c(0,20))

hist(scale(X),prob=T,xlim=c(-5,5))

#plot(x,dbinom(x, 10, 1/6),)
p<-1/6
x<-0:20
prob<-dbinom(x, length(x)-1, p)
#names(prob)<-x
plot(prob,col=1,type="b")

5 例題７（p86）

5.1 母集団

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