for文を使いこなす

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()
}