for文を使えるようになると，いろいろ便利です．以下のような図も作れます．

◆準備

setwd()でワーキングディレクトリを指定し，演習用フォルダを指定しておいてください．参考：地理空間情報の可視化

◆二項分布に親しむ

表が出る確率がpのコインをn回投げて，表がk回出る確率の分布．

― 二項分布の数式 ―

問1 パラメータがそれぞれ \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において，表が出る回数を \(k\) と表す．このとき， \(k=20\)となる確率を求めよ．

ヒント：dbinom(k,n,p)という関数を使うことで，\(B(n,p)\) において\(k\)が出る確率を求めることができる．

dbinom(k,n,p)

問2 パラメータがそれぞれ \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) から，20個の乱数をサンプリングせよ．

ヒント：rbinom(size,n,p)という関数を使うことで，\(B(n,p)\) からsize個の乱数をサンプリングすることが出来る．

rbinom(size,n,p)

◆for文を使ってみる！

for文は，以下のような構文である．Rのfor文では，()内のiに，vectorの値をひとつずつ入れて， {} 内の処理をおこなう．たとえば，for(i in 1:5)であれば，i=1, i=2, i=3, i=4, i=5をそれぞれ代入し，{}内のprint(i)を実行する．print関数は()内のオブジェクトをコンソールに吐き出す関数である．

printをしないでiだけ書くとどうなるかも試してみよう．

for (i in vector) {
  print(i)
}

問3 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において、表が出る回数\( k \)が， \( 0,1,2,…50 \) となる確率をすべて計算せよ．

ヒント：もちろん，1行ずつkを変えて…となると大変．for文を使おう．????に繰り返したい値をベクトルに入れて，指定する．i,n,pもそれぞれ指定する．print関数をつけないとどうなるか…？

for (i in ????) {
  dbinom(i,n,p)
}

問4 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) から、20個の乱数をサンプリングする，という操作を100回繰り返し，その結果を表示せよ．

ヒント： ????に繰り返したい値をベクトルに入れて，指定する．i,n,pもそれぞれ指定する．print関数を使って表示する．

for (j in ????) {
  rbinom(size,n,p)
}

◆for文を使って可視化する

for文を使うと，条件を少しずつ変えたグラフなども簡単に作れる．

問5 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において、表が出る回数\( k \)が， \( 0,1,2,…50 \) となる確率をすべて計算し，x軸にkを，y軸にkが出る確率をとった二次元のグラフを作成してみよう．

ヒント1：これもfor文を使うこと（for文を使わない方法もある）．

ヒント2：グラフを出力する関数はplot()である．plot(x,y)で，横軸がx，縦軸がyのグラフが描かれる（xとyのベクトルの長さは同じでなければならない）．

ヒント3：par(new=T)をplot関数の前に指定することで，その前に書いたグラフに重ね書きができる．

for (i in ????) {
  plot(i,dbinom(i,n,p))
  par(new=T)
}

ヒント4：x軸とy軸がぐちゃぐちゃになったのではないでしょうか？そういう時は，以下のようにxlim，ylim というパラメータを指定することで，毎回同じ範囲の軸を描くことが出来る．1回目だけ描いて2回目以降は描かないように指定することで，重ね書きによって文字が分厚くなることも防げそうだが，今回はそれはやらない．

xlim = c(0,50),ylim = c(0,0.2)

問6　 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) から、20個の乱数をサンプリングする，という操作を100回繰り返し，その結果をヒストグラムに書き表してみよう．

ヒント：breaksというパラメータを使って，ヒストグラムの範囲と幅を指定することが出来る．たとえば，以下のようにパラメータ指定すると，0から50の値を幅1で表示する，というものになる．

breaks=seq(0,50,1)

ヒント2：黒々として見にくければ，colというパラメータで色を分けて指定することもできる．col = i のように指定すれば，毎回色を変えて可視化することができる．

問7　 \(n=50, p=0.3\) の 二項分布\(B(50,0.3)\) において、表が出る回数\( k \)が， \( 0,1,2,…50 \) となる確率を計算し，その確率に各々50をかけた値をy軸にして表示せよ．また，これを問6のグラフに重ね合わせてみよう．

◆ファイルの出力

問8　ファイルに出力する．

まずは，setwdで保存したいフォルダのパスを指定する．

setwd("保存したいフォルダのパス")
getwd()

ファイルに画像データを保存する場合は， plot()の前 にpng(“パス/ファイル名”)を実行し， plot()の実行後，dev.off()を実行することで，plot()で描いたグラフが “パス/ファイル名” に保存される．

ファイル名が同じだと，毎回上書きされてしまう．これでは困る．ファイル名を各繰り返しの度に変えたい．

そこで，paste()という関数を使う．これは，paste(“Hides”,”hima”)のように実行すると”Hideshima”のように，文字列をくっつけてくれる関数である．これを使って，毎回異なるファイル名を指定できるようにしてみよう．

sep = “”　というパラメータをつけないと， “Hides hima” のようにスペースが入る．

paste("Hides","hima",sep = "")

実行後，指定したフォルダ内を確認してみよう．

for (i in 0:50) {
  png(paste("binom",i,".png"))
  hist(rbinom(20,50,0.3),xlim = c(0,50),ylim = c(0,20),breaks=seq(0,50,1),col = i)
  dev.off()
}

◆動画の作成

「 png　複数　gif化 」等で各自ネットで調べてください．

◆演習の解答

以下は解答例です．もっといいやり方もあると思うので，各自模索してみてください．

下記コードのダウンロードはこちらから．

# 2020.6.3
# coded by F.Nakai

#二項分布
### 肩慣らし
#Q1: p=0.3, N=50の二項分布において、k=20となる確率
dbinom(20,50,0.3)

#Q2: p=0.3, N=50の二項分布から、20個の乱数をサンプリングする
rbinom(20,50,0.3)

### for文を使ってみる！
#Q3: p=0.3, N=50の二項分布において、K=0,1,2,...50となる確率を計算する

for (i in 0:50) {
  print(dbinom(i,50,0.3))
}

#Q4:p=0.3, N=5の二項分布から、20個の乱数をサンプリングするのを100回繰り返す

for (j in 1:100) {
  print(rbinom(20,50,0.3))
}

### for文を使って可視化してみる
#Q5: p=0.3, N=50の二項分布において、K=0,1,2,...50となる確率を計算して可視化する
for (i in 0:50) {
  print(dbinom(i,50,0.3))
  plot(i,dbinom(i,50,0.3),xlim = c(0,50),ylim = c(0,0.2))
  par(new=T)
}


#Q6:p=0.3, N=50の二項分布から、20個の乱数をサンプリングするのを100回繰り返す

par(new=F)
for (j in 1:100) {
  print(rbinom(20,50,0.3))
  hist(rbinom(20,50,0.3),xlim = c(0,50),ylim = c(0,20),breaks=seq(0,50,1),col = j)
  par(new=T)
}

#Q7: p=0.3, N=50の二項分布において、K=0,1,2,...50となる確率を計算し、回数に変換して可視化する

par(new=T)
for (i in 0:50) {
  print(dbinom(i,50,0.3))
  plot(i,dbinom(i,50,0.3)*50,xlim = c(0,50),ylim = c(0,20))
  par(new=T)
}

#Q8: ファイルにつくったグラフを保存する．

#出力するフォルダを指定する．
setwd("保存したいフォルダのパスを設定")
getwd()

paste("hides","hima",sep = "")

#for文で1つずつグラフを書いていく
par(new=T)
for (i in 0:100) {
  png(paste("binom",i,".png"))
  hist(rbinom(i*10,50,0.3),xlim = c(0,50),ylim = c(0,300),breaks=seq(0,50,1),col = "gray")
  dev.off()
}

#maptoolsライブラリを使ってみる
#グラフのプロットにラベルを付ける
library(maptools)
mes <- c("今","日","は","お","つ","か","れ","さ","ま","で","し","た","!")
mes1 <- length(rep(mes,5))

for (i in seq(1,50,1)) {
  png(paste("binom",i,".png"))
  plot(i,dbinom(i,50,0.3),xlim = c(0,50),ylim = c(0,0.2),col = "white",pch = 19,cex = 10)
  pointLabel(x=i, y=dbinom(i,50,0.3), labels=mes[i])
  dev.off()
}