ベイズ統計R演習1

演習の概要

コイントス、サイコロ投げのシミュレーション
一様分布、二項分布、ベータ分布、ポアソン分布、正規分布

注意事項

慶應義塾大学SFCで開講している「ベイズ統計」の授業履修者向けの演習用ページです。
必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.5.3を使っています。
Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

Rの利用方法

フリー統計ソフトRを利用するには、Rをダウンロードする、Webブラウザ上で利用する、などの方法があります。
RはCRANからダウンロードし、インストール可能です。
RStudioには、より操作が容易なGUIが用意されています。
RStudio Cloudを使うと、RをダウンロードすることなくWebブラウザ上でRを操作できます。
RStudioは日本語が使えない場合があります。起動時にLC_CTYPE failedなどのエラーメッセージが出る場合は、以下の一行を実行しRStudioを再起動てください。

system("defaults write org.R-project.R force.LANG en_US.UTF-8")

コイントスとサイコロ投げ

コイントス

100回のコイントスを実行します
乱数はsample()で生成します
シミュレーションなので、結果は各自で異なります

# コイントス回数
N <- 100
# 確率0.5で表と裏がでる確率をN回ランダムに行う
cointoss1 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
# 全てのデータを見てみる
cointoss1

##   [1] 1 1 1 0 1 1 1 1 1 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 1 0 1 1 0 1
##  [36] 0 1 1 1 0 0 1 0 0 1 0 1 0 0 1 0 1 1 0 1 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0
##  [71] 1 0 1 0 0 0 0 1 1 1 1 1 0 0 1 0 1 1 1 1 1 0 0 0 1 0 1 0 0 0

# 最初の20回分をみてみる
head(cointoss1, 20)

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

# 最後の20回分をみてみる
tail(cointoss1, 20)

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

# 図をプロットする際の日本語フォントを指定
par(family="HiraMaruProN-W4")
# plot(x, y, ...)：図をプロット
# cumsum()：累積和
plot((1:N), cumsum(cointoss1)/(1:N) , 
# 点をoで描画
type="o",
# x軸とy軸の範囲
xlim=c(1,N) , ylim=c(0.0,1.0),
# 図タイトルおよびx軸とy軸の名前
main="コイントスで表が出る確率 (N=100)" , xlab="コイントス回数" , ylab="表が出る確率",
# xy軸、xy軸名前、図タイトルの大きさの指定
cex.axis=1.2, cex.lab=1.2, cex.main=1.5) 
# 確率0.5の線の描画
# col=""：色を指定
lines(c(1,N),c(0.5,0.5), col="red", cex=1.2)

# 最初の20回の結果を表示
head(cumsum(cointoss1)/(1:N), 20)

##  [1] 1.0000000 1.0000000 1.0000000 0.7500000 0.8000000 0.8333333 0.8571429
##  [8] 0.8750000 0.8888889 0.8000000 0.8181818 0.7500000 0.6923077 0.6428571
## [15] 0.6000000 0.6250000 0.5882353 0.5555556 0.5789474 0.5500000

# 最後の20回の結果を表示
tail(cumsum(cointoss1)/(1:N), 20)

##  [1] 0.4691358 0.4756098 0.4698795 0.4642857 0.4705882 0.4651163 0.4712644
##  [8] 0.4772727 0.4831461 0.4888889 0.4945055 0.4891304 0.4838710 0.4787234
## [15] 0.4842105 0.4791667 0.4845361 0.4795918 0.4747475 0.4700000

10000回のコイントス

N <- 10000
cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
par(family="HiraMaruProN-W4")
plot((1:N), cumsum(cointoss2)/(1:N) ,  type="o", log="x", xlim=c(1,N) , ylim=c(0.0,1.0),
main="コイントスで表が出る確率 (N=10000)" , xlab="コイントス回数" , ylab="表が出る確率",
cex.axis=1.2, cex.lab=1.2, cex.main=1.5) 
lines(c(1,N),c(0.5,0.5), col="red", cex=1.2)

# 最初の20回の結果を表示
head(cumsum(cointoss2)/(1:N), 20)

##  [1] 1.0000000 1.0000000 0.6666667 0.5000000 0.4000000 0.3333333 0.2857143
##  [8] 0.2500000 0.3333333 0.4000000 0.3636364 0.4166667 0.4615385 0.5000000
## [15] 0.5333333 0.5000000 0.5294118 0.5000000 0.4736842 0.4500000

# 最後の20回の結果を表示
tail(cumsum(cointoss2)/(1:N), 20)

##  [1] 0.5017533 0.5018032 0.5017530 0.5018029 0.5018528 0.5018025 0.5017523
##  [8] 0.5017020 0.5017519 0.5018018 0.5018517 0.5019015 0.5019514 0.5020012
## [15] 0.5019510 0.5019008 0.5019506 0.5020004 0.5019502 0.5020000

10000回のコイントスを５回行う

N <- 10000
cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
par(family="HiraMaruProN-W4")
plot((1:N), cumsum(cointoss2)/(1:N) ,  type="l", xlim=c(1,N) , ylim=c(0.0,1.0),
main="コイントスで表が出る確率 (N=10000)" , xlab="コイントス回数" , ylab="表が出る確率",
cex.axis=1.2, cex.lab=1.2, cex.main=1.5, col=1, lwd=5) 
cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
lines((1:N), (cumsum(cointoss2)/(1:N)), col=2, lwd=3)
cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
lines((1:N), (cumsum(cointoss2)/(1:N)), col=3, lwd=1.5)
cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
lines((1:N), (cumsum(cointoss2)/(1:N)), col=4, lwd=1)
cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
lines((1:N), (cumsum(cointoss2)/(1:N)), col=5, lwd=1)

サイコロ投げ

100回のサイコロ投げ

N=100
# seq()：連続した整数
# rep(x,n): xをn回繰り返す
diceroll <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE)
head(diceroll,20) # 最初の20個を表示する

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

diceroll1 <- diceroll[diceroll==1] # 1の目が出たものを取り出す
diceroll1 # すべて表示してみる

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

length(diceroll1) # 個数を数える

## [1] 15

table(diceroll) # クロス集計

## diceroll
##  1  2  3  4  5  6 
## 15 20 15 16 17 17

prop.table(table(diceroll)) # 割合を計算

## diceroll
##    1    2    3    4    5    6 
## 0.15 0.20 0.15 0.16 0.17 0.17

# ヒストグラムで度数を表示
par(family="HiraMaruProN-W4")
hist(diceroll, breaks=seq(0,6), freq=TRUE, ylim=c(0,50), 
main="サイコロ投げ 100回", xlab="サイコロの目", ylab="度数")
lines(c(0,6),c(1/6*N,1/6*N), col="red", cex=1.2)
text(5.5,1/6*N+1, "理論確率", col="red")

# ヒストグラムで確率を表示
par(family="HiraMaruProN-W4")
hist(diceroll, breaks=seq(0,6), freq=FALSE, ylim=c(0,0.4), 
main="サイコロ投げ 100回", xlab="サイコロの目", ylab="確率")
lines(c(0,6),c(1/6,1/6), col="red", cex=1.2)
text(5.5,0.18, "理論確率", col="red")

10000回のサイコロ投げ

N=10000
diceroll <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE)
# head(diceroll,20) # 最初の20個を表示する
table(diceroll) # クロス集計

## diceroll
##    1    2    3    4    5    6 
## 1662 1621 1661 1658 1725 1673

prop.table(table(diceroll)) # 割合を計算

## diceroll
##      1      2      3      4      5      6 
## 0.1662 0.1621 0.1661 0.1658 0.1725 0.1673

# ヒストグラムで度数を表示
par(family="HiraMaruProN-W4")
hist(diceroll, breaks=seq(0,6), freq=TRUE, ylim=c(0,(1/6*N+1000)), 
main="サイコロ投げ 10000回", xlab="サイコロの目", ylab="度数")
lines(c(0,6),c(1/6*N,1/6*N), col="red", cex=1.2)

# ヒストグラムで確率を表示
par(family="HiraMaruProN-W4")
hist(diceroll, breaks=seq(0,6), freq=FALSE, ylim=c(0,0.3), 
main="サイコロ投げ 10000回", xlab="サイコロの目", ylab="確率")
lines(c(0,6),c(1/6,1/6), col="red", cex=1.2)
text(5.5,0.18, "理論確率", col="red")

ガリレオのサイコロ問題をシミュレーションで確かめてみる

N=10000
dice.G <- matrix(0, ncol=4, nrow=N)
colnames(dice.G)=c("dice1", "dice2", "dice3", "Total")
dice.G[,1] <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE) #1個目のサイコロ
dice.G[,2] <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE) #2個目のサイコロ
dice.G[,3] <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE) #3個目のサイコロ
dice.G[,4] <- rowSums(dice.G[,1:3]) # 行和
table(dice.G[,4]) # クロス集計

## 
##    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17 
##   45  167  260  468  640  992 1120 1255 1217 1236  996  652  468  290  141 
##   18 
##   53

prop.table(table(dice.G[,4])) # 割合を計算

## 
##      3      4      5      6      7      8      9     10     11     12 
## 0.0045 0.0167 0.0260 0.0468 0.0640 0.0992 0.1120 0.1255 0.1217 0.1236 
##     13     14     15     16     17     18 
## 0.0996 0.0652 0.0468 0.0290 0.0141 0.0053

基本的な確率分布

一様分布

区間[0, 1]で1となる一様分布を描画する
一様分布の確率密度関数はdunif()で計算できます

x <- seq(-1,2,length.out=1000)
par(family="HiraMaruProN-W4")
plot(x, dunif(x, min=0,max=1), xlim=c(-1,2), ylim=c(0,2), type="l", main="一様分布")

区間[0, 0.5]で2となる一様分布

x <- seq(-1,2,length.out=1000)
par(family="HiraMaruProN-W4")
plot(x, dunif(x, min=0,max=0.5), xlim=c(-1,2), ylim=c(0,2), type="l", main="一様分布")

二項分布

n=10回の試行でs=1 (r=0.1)となるベルヌーイ試行を描画する
二項分布の確率密度関数はdbinom()で計算できます

par(family="HiraMaruProN-W4")
plot(dbinom(seq(0,10,1), 10, 0.1), type="l", ylab="確率")
text(3, 0.37, labels="B(10,0.1)")

n=10回の試行でs=3 (r=0.3)となるベルヌーイ試行

par(family="HiraMaruProN-W4")
plot(dbinom(seq(0,10,1), 10, 0.3), type="l", ylab="確率")
# lines(dbinom(seq(0,10,1), 10, 0.3))
text(6, 0.2, labels="B(10,0.3)")

n=10回の試行でs=5 (r=0.5)となるベルヌーイ試行

par(family="HiraMaruProN-W4")
plot(dbinom(seq(0,10,1), 10, 0.5), type="l", ylab="確率")
# lines(dbinom(seq(0,10,1), 10, 0.5))
text(8, 0.2, labels="B(10,0.5)")

n=10回の試行でs=7 (r=0.7)となるベルヌーイ試行

par(family="HiraMaruProN-W4")
plot(dbinom(seq(0,10,1), 10, 0.7), type="l", ylab="確率")
# lines(dbinom(seq(0,10,1), 10, 0.7)) 
text(6, 0.2, labels="B(10,0.7)")

ベータ分布

様々なベータ分布を描画する
ベータ分布の確率密度関数はdbeta()で計算できます

# ベータ分布Β(0.5, 0.5)
x <- seq(0,1,length=100)
par(family="HiraMaruProN-W4")
curve(dbeta(x, 0.5,0.5), type="l", xlim=c(0,1), ylim=c(0, 8), ylab="確率", lwd=2)
#  Β(0.5, 1) ⇒ 単調減少分布となる
curve(dbeta(x, 0.5, 1), type="l", xlim=c(0,1), col="black", add=TRUE, lty=2, lwd=2)
# Β(1, 2) ⇒ 単調増加分布となる
curve(dbeta(x, 1, 2), type="l", xlim=c(0,1), col="grey", add=TRUE, lwd=2)
#  Β(1, 1) ⇒ 一様分布となる
curve(dbeta(x, 1, 1), type="l", xlim=c(0,1), col="red", add=TRUE, lwd=2)
# Β(5, 5) ⇒ 0.5をピークとする釣鐘型の分布となる
curve(dbeta(x, 5, 5), type="l", xlim=c(0,1), col="blue", add=TRUE, lwd=2)
# Β(10, 2)
curve(dbeta(x, 10, 2), type="l", xlim=c(0,1), col="green", add=TRUE, lwd=2)
# Β(10, 12)
curve(dbeta(x, 10, 12), type="l", xlim=c(0,1), col="yellow",  add=TRUE, lwd=2)
# 凡例
legend(.5,8,c("Β(0.5, 0.5)","Β(0.5, 1)","Β(1, 2)","Β(1, 1)","Β(5, 5)","Β(10, 2)","Β(10, 12)"), 
col=c("black","black","grey","red","blue","green","yellow"),lwd=c(2,2,2,2,2,2,2),  lty=c(1,2,1,1,1,1,1))

ポアソン分布

ポアソン分布Po(1), Po(3), Po(5), Po(10)を描画する
ポアソン分布の確率密度関数はdpois()で計算できます

# ポアソン分布Po(1)
par(family="HiraMaruProN-W4")
plot(dpois(0:50, lambda=1), type="b", xlim=c(0,20), lwd=2, ylab="確率", xlab="x")
text(3.5, 0.35, "Po(1)")
# Po(3)
lines(dpois(0:50, lambda=3), type="b", lwd=2, col=2)
text(5, 0.24, "Po(3)", col=2)
# Po(5)
lines(dpois(0:50, lambda=5), type="b", lwd=2, col=3)
text(7, 0.19, "Po(5)", col=3)
# Po(10)
lines(dpois(0:50, lambda=10), type="b", lwd=2, col=4)
text(11, 0.14, "Po(10)", col=4)

正規分布

様々な正規分布を描画する
正規分布の確率密度関数はdnorm()で計算できます

# 区間[-100, 100]で平均0, 標準偏差20の正規分布
x <- seq(-100,100,length=2000)
plot(x, dnorm(x, 0,20), type="l", ylim=c(0,0.02), xlim=c(-100,100), lwd=2)
# 平均10, 標準偏差20の正規分布
curve(dnorm(x, 10,20), type="l", ylim=c(0,0.02), xlim=c(-100,100), lwd=2, col="red", add=TRUE)
# 平均0, 標準偏差40の正規分布
curve(dnorm(x, 0,40), type="l", ylim=c(0,0.02), xlim=c(-100,100), lwd=2, col="blue", add=TRUE)
# 凡例
legend(30,0.02,c("N(0,20^2)","N(10,20^2)","N(0,40^2)"), col=c("black","red","blue"),lwd=c(2,2,2))

ガンマ分布

Ga(1, 1)とGa(1, 3)を描画する
ガンマ分布の確率密度分布はdgamma()で計算できます

# Ga(1, 1)
# Ga(α, λ) ⇒ dgamma(x, shape, rate, scale=1/rate)
# α⇒shape, λ=rate
x <- seq(0,10,by=0.01)
# α=1, λ=1⇒指数関数と同じになる
par(family="HiraMaruProN-W4")
plot(x, dgamma(x, 1, 1), type="l", lwd=2, ylab="確率")

# Ga(1, 3)
curve(dgamma(x, 1, 3), type="l", lwd=2)

Ga(5, 1)

x <- seq(0,20,by=0.01)
plot(x, dgamma(x, 5, 1), type="l", lwd=2)

Ga(40, 1)

plot(x, dgamma(x, 40,1), type="l",  lwd=2)

Ga(2, 0.5)

x <- seq(0,20,by=0.01)
plot(x, dgamma(x, 2, 0.5), type="l", lwd=2)