Rを使ったプログラミング演習2：　シミュレーション †

前回からRを使ったプログラミングについて講義を行っている。前回、プログラミングの基本である、繰り返しと代入について解説したので、今回は残りの条件分岐についての解説を行う。さらに、この３つの基本技を使って、シミュレーションに挑戦する。

なお、アンケート調査によると、プログラミングが難しいという意見がたくさんみられた。初めてプログラミングに触れる人は、いろんな命令文がでてくるので、なかなか「理解した」という実感がえられないかもしれない。でも、まずは、自分で作った一連の命令（つまりプログラム）を、コンピュータに与えらることができるということを実感してほしい。また、繰り返し・代入・条件分岐というプログラムの構造は理解してほしい。その他の個々の命令（matrixやplotなど）は、必要に応じて説明文書などを参照して、使い方を調べることができるようになればよしとしよう。

↑

第10回授業の獲得目標： [worried] †

↑

条件分岐 [smile] †

プログラミングの基本技の最後は条件分岐。if()を使って表現する。

i=4                       #iに4を代入
if(i==3){                 #iの値が3ならば
  print("三")            #   三 と表示する
} else {                  #iの値が3では無ければ
  print("三以外")   #   三以外と表示する
}

if命令の()の中の式を条件式といい、iの値を評価している。評価に使われるのは比較演算子というもので、

==　    等しければ
!=      等しく無ければ
>, >=   左辺が右辺より大きいなら、左辺が右辺以上ならば
<, <=   左辺が右辺より小さいなら、左辺が右辺以下ならば

ということを意味している

''比較演算子''である''=='' に対して、代入の時に使った ''='' は''代入演算子''といいう。

上のif命令では、次のような条件分岐が行われている。

if(i==3){
   print("三")  # ・もしiの値が３ならば、画面に三を表示。
} else {                  #もしiの値が3では無ければ
  print("三以外")   #   三以外と表示する
}
↑

乱数を使った条件分岐の演習 †

上の例では、iの値が最初から決まっているので、出てくる結果は１通りしかなく、プログラムとしては面白く無い。そこで、乱数を１つ発生させ、値が0.5以上なら「0.5以上」、0.5より小さいなら「0.5より小さい」と表示させるプログラムを考えてみよう。
乱数を1つ発生させるにはrunif()という関数を使う。試しに

runif(5)

と入力してみよう。0から1までの乱数が5個表示されたはずだ。乱数はシミュレーションではよく使うので、もう少し慣れ親しんでおこう。 例えば、runif()を使って乱数を10個発生させ、ヒストグラムを書いてみよう。

hist(runif(10))

同じことを、100個、1000個、10000個、100000個についてもやってみよう。バラツキがだんだんと無くなってきたのがわかるだろうか。

#穴埋め式練習問題
if (runif(1) >= _ ) {   #条件式0.5以上ならば
  print("0.5以上")
} _ { 
  print("0.5より小さい")
}

これでは、どの乱数が発生したのか分からないので、一度別の変数に代入しておいて、判定基準と一緒に表示させてみよう。

#穴埋め式練習問題: _の部分に適当な文字を入力して実行しよう
x=runif(1)      #乱数を1つ発生させてxに代入
kekka=c()     #結果をいれておく空ベクトルを作る
if (_>= 0.5 ) {   #条件式0.5以上ならば
  kekka=c(kekka, x)
  kekka=c(kekka, "は0.5以上")
} else { 
  kekka=c(kekka, x)
  kekka=c(kekka, "は0.5より小さい")
}
  print(kekka)
↑

繰り返しと条件分岐 †

　では、上の繰り返しと条件分岐を組み合わせて使ってみよう

練習問題：for命令を使って1から10までの数字を表示させなさい、5だけは"five"と英単語で表示させなさい
↑

関数の定義 [smile] †

「関数の定義」なんていうと難しそうだが、ようするに、先ほどまでに作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにしようというものだ。こういうときに関数を定義する関数、function()を使う。 例えば、円の面積を計算する関数を作るなら、

 menseki=function(r){r*r*pi}

とすれば良い。これで自分独自の関数menseki()ができた。 では、実行してみよう。

menseki(10)  #好きな数字を入れて円の面積を計算

演習：1から入力した数値までの全てを横一列に表示させるプログラムを作りdisplayという名前の関数として定義する

display=function(a){       #関数定義の始まり
 kekka=c()                 #kekkaに空ベクトルを代入して初期化
 for(i in 1:a){              #a回（iの値を1からaまで変化させる間）繰り返し
   kekka=c(kekka,i)     #kekkaというベクトルにiを要素として代入
 }
 print(kekka)              #kekkaの内容を表示
}                          #関数定義の終わり
display(10)

↑を使って何回も実行してみると、結果がいろいろ変わるのがわかる。

演習上の関数の定義方法に従って、入力した数までの合計値を計算するsumupという名前の関数を作成する。下の囲みの中の_の部分（１文字に対応するとは限らない）を埋めて、プログラムを完成しなさい。

 sumup = ______       #関数sumpuを定義
   kotae = _          #kotaeを初期化
   _ (i in ___){      #1からaまで繰り返し
     kotae = _____    #kotaeにiの値を足したものとkotaeに代入
   }                  #繰り返し終了
  _____               #kotaeを表示
  }                   #関数定義の終了
↑

繰り返し・代入・条件分岐のおさらい †

次のそれぞれのプログラムの＿部分に何を入れるべきかを選択肢から選びなさい。また、Rを使って実行してみましょう。

↑

繰り返し †

for(i in 1:__){    #{と}の間を10回繰り返す。i の値は毎回1ずつ増える
  print(i)
}

選択肢： a) 1    b)  5     c)  10
↑

代入 †

goukei=0           #goukeiという変数に初期値0を代入
for(i in 1:10) {   
  goukei=_____+i  #goukeiに前回までのgoukeiの値にiの値を足したもの代入
} 
print(goukei)

選択肢： a) i    b)  goukei   c)  print
↑

条件分岐 †

goukei=0           #goukeiという変数に初期値0を代入
 for(i in 1:10) {
   if(i %% 2 __ 0) { # %% は割り算をした余りを返す演算子。3 %% 2 の計算結果は 1 になる
     goukei=goukei+i # if文の内容は：　「もし i を 2 で割った余りが 0 ならば、{}内を実行
   }               #それ以外は何もしない
 }              #for 文のカッコ閉じるに注意
 print(goukei)

選択肢： a) =    b)  ==     c) <=
↑

授業で使ったRの基本関数 †

授業で使うR関数一覧 †

↑

演算子 †

演算子意味使用例
+足し算> 4+3
[1] 7
-引き算> 4-3
[1] 1
*掛け算> 4*3
[1] 12
/割り算> 4/3
[1] 1.333333
^累乗> 4^3
[1] 64
%/%整数商> 7%/%3
[1] 2
%%剰余> 7%%3
[1] 1
↑

定数（Rの基本システムで決まっている値):　次の５つ †

> LETTERS
  [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U"
 [22] "V" "W" "X" "Y" "Z"
> letters
 [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s" "t" "u"
 [22] "v" "w" "x" "y" "z"
> month.abb
  [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"
> month.name
  [1] "January"   "February"  "March"     "April"     "May"       "June"      "July"     
  [8] "August"    "September" "October"   "November"  "December" 
> pi
 [1] 3.141593
↑

関数 †

Rの関数意味使用例
c()ベクトルの作成x = c(1, 2, 3, 4)
mean()ベクトルの平均mean(1, 2, 3, 4, 5)
scan()コピー・ペーストや手入力でデータをベクトルとして読み込む。x=scan()
hist()ヒストグラムを表示x=c(1,1, 3, 1, 5, 3, 4, 4); hist(x)
summary()要約統計量の表示x=c(1,1, 3, 1, 5, 3, 4, 4); summary(x)
plot()散布図の表示x=c(1,1, 3, 1, 5, 3, 4, 4); plot(x)
print()オブジェクトを表示x=c(1,1, 3, 1, 5, 3, 4, 4); print(x)
for(ループ変数 in リスト) {式}ループ変数をリストの内容に従って、１つずつ変化させる。変化の1回ごとに、{ }の式を実行。for(i in 1:10) { print(i) } #注：1:10は 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 を表す
matrix(ベクトル, 行数, 列数)ベクトルを行列に変換して表示matrix(1:6, nrow=2, ncol=3)
if (条件式) {式1} else {式2}条件式が真の時は 式1を実行、偽のときは式2を実行if (runif(1) < 0.5) { print("0.5より小さい")} else { print("0.5以上")}
runif()0-1の乱数を指定した個数発生runif(5)
function(引数1, 引数2,....){式 }自分で関数を定義するmenseki=function(r){r*r*pi}
plot.new()図を新しく書き直すという関数
↑

Rを用いたシミュレーション [smile] †

これまでに学んだ繰り返し・代入・条件分岐の３つの命令それぞれは単純だが、組み合わせれば、かなり複雑なこと表現できる。そこで、今回は、これらの命令を組み合わせて、Rを用いたシミュレーションに挑戦してみよう。
シミュレーションという言葉は「株式シミュレーション」とか「シミュレーションゲーム」とかで聞いたことがあるだろう。日本語にすると、「模擬実験」で、自分では実行が難しい繰り返し実験を、コンピュータに計算させることで行うものだ。コンピュータは単純な繰り返し計算は得意なので、砂粒をばらまいて、ある図形の中に入った数を数えるなんていう作業を、10万回でも、100万回でも繰り返すことができる。

↑

円周率(π）を「繰り返し」・「乱数」・「条件分岐」で求める †

#10_1.jpg
左に１辺の長さが1の正方形と、それに内接する扇形がある。このとき、

正方形の面積は　 1　　扇形の面積は       π/4

この図形の上に、砂粒をばらまいてみる。 砂粒が図形の上にランダムに散らばるとすると、砂粒をすごく沢山まいて数を数えれば、

正方形の中の砂粒の数(n) ： 扇形の中砂粒の数(m)  ＝　
正方形の面積 ：扇形の面積 ＝
１　：　π/4

になると考えられる。つまり

   n : m = 1 : π/4　　　だから、   π = 4m/n

もしあなたが、砂粒をたくさんランダムにまいて、正方形の中の砂粒と扇形の中の砂粒の数を数えることができれば、πの値がわかるということになる。

↑

シミュレーションの考え方 †

さて、砂粒を沢山まいて、それが円の中か外かを判定する実験を、コンピュータにやらせるにはどうすればよいだろうか。 まず、点（砂粒）を大量に発生させる必要がある。さて、どうすればいいだろうか？
もうお分かりのように、「繰り返し」命令（for文）を使う。

#とりあえず、1,000つぶの砂をまくことをfor文で表現する
for (i in 1:1000){
 ＜砂を１粒図形の上にまく＞
}

次に、砂を１粒図形の上にまいて、それが扇形の中か外かを判定することを考える。

砂を１粒まくことは、無作為に（ランダムに）図形の上に点を打つことと同じことだ。「ランダム」というと、我々はすでに乱数を発生させる方法をしっている。図形の上に１つの点を打つということは、２つの乱数を発生させることで表現できるだろう。

#10_2.jpg
例えば、0-1の範囲の乱数を２つ発生させたら、　(0.230 , 0.782) という２つの数字が得られたとする。この２つの数字が座標上の点を表していると考えると、上の図のようになる。
原点からこの点までの距離が1より小さいとき、この点は扇形の中に入っているということができる。この例では、

0.230^2+0.782^2 < 1

なので、ランダムに発生させた（図形の上に偶然落ちた）1つの点は、扇形の中だったということになる。

point=runif(2)  #この命令で、0-1の範囲の乱数を２つ発生させ、pointというベクトルに入れる

このとき、pointに入ったそれぞれの値を指定するには、

point[1]
point[2]

のようにカギ括弧［］を使う。

例：　プロンプト以降を入力してみよう
        > x=c(0.230 , 0.782) 
        > x[1]
        > x[2]

砂粒（点）が扇形の中かどうかを判断するにはどうすれば良いか？これは条件分岐で表現できる。つまり、runif(2)で発生させた２つの数値のそれぞれを2乗して足しあわせたものの平方根が1より小さいか否かを判断する。

#扇形の中かどうかを判断する条件文
 distant=sqrt(point[1]^2+point[2]^2)   #発生させた２つの乱数の２乗の和の平方根
   if(distant<1){　　　　　　#原点からの距離が1より小さければ扇形の中
       ＜個数を数える＞
  }

では、個数を数えるにはどうすればよいか？ここでは代入文を使う

 #個数を数える部分
 kaisu=0   #回数の初期値を0にする
   if(distant<1){　　　　　　#原点からの距離が1より小さければ扇形の中
       kaisu = kaisu + 1   #条件にあったとき、kaisuの値を1増やす
  }

以上のように、これまでに学習した、繰り返し、条件分岐、代入で、図形の上に砂粒を1000個まく実験が表現できた。πの値は砂粒の個数を4倍して、全体の個数で割ったものだから、以上をまとめると次のようになる。

	kaisu=0;  #個数の初期値を0にする
	for (i in 1:1000) {
		point=runif(2)
		distant=sqrt(point[1]^2+point[2]^2)
		if(distant<1){
			kaisu=kaisu+1
		}
	}
	answer=4*kaisu/1000
	print(answer)

さらに、functionを使って、好きな回数繰り返せる関数を定義してみよう

simpai=function(x){
	kaisu=0;  #個数の初期値を0にする
	for (i in 1:x) {
		point=runif(2)
		distant=sqrt(point[1]^2+point[2]^2)
		if(distant<1){
			kaisu=kaisu+1
		}
	}
	answer=4*kaisu/x
	print(answer)
}
> simpai(10)
[1] 2.8
> simpai(100)
[1] 3.16
> simpai(1000)
[1] 3.156
> simpai(10000)
[1] 3.1624
> simpai(100000)
[1] 3.15652
> simpai(1000000)   #ここから先は時間がかかるので、授業中はやらない方がいい
[1] 3.13994
> simpai(10000000)  #ここから先は時間がかかるので、授業中はやらない方がいい
[1] 3.140556
↑

せっかくだから、点をうつところを図で表現してみよう †

simpai=function(x){
	kaisu=0;  #回数の初期値を0にする
	for (i in 1:x) {
		point=runif(2)
               par(new=T)    #複数のグラフを重ね合わせる
               plot(point[1], point[2], ylim=c(0,1), xlim=c(0,1))  #0-1の範囲で乱数の値をプロット
		distant=sqrt(point[1]^2+point[2]^2)
		if(distant<1){
			kaisu=kaisu+1
		}
	}
	answer=4*kaisu/x
	print(answer)
}

simpai100000.jpg
10万回の実行結果。

↑

袋の中から赤玉・青玉を取り出す実験のシミュレーション †

赤玉と青玉が同じ数だけたくさん入った筒から、20個の玉を取り出して、その中に含まれる赤玉の割合に従って筒に含まれる赤玉の個数を変え、また20個取り出す。。。こういう操作を何度も繰り返すというモデルで、何度も実験すると、袋の中の玉の割合がどのように変化するかをシミュレーションで示してみよう。

質問：　この実験を何度も繰り返しておこなうとき、
筒の中の玉が赤か白かどちらか１色になってしまうのは、何回目だろうか？

Untitled-9.gif

「赤玉と青玉の同じ数入っている筒から無作為に20個の玉を取り出したら、赤玉が7個入っていた」
  -----------------------------------------------------
上でやった操作には、「20個取り出して赤玉かどうか判断する」というところに「繰り返し」が含まれている。
「繰り返し」がはっきりと分かるように説明すると、
 1. 筒からボールを一つ取り出す
 2. そのボールが赤玉だったら数える（１つ目だったら「一つ」、２つめは「２つ」..)。青玉だったら数えない
 3.  1 と 2 の操作を20回繰り返す

この部分を大雑把でいいので、プログラミングしてみる。繰り返しはfor文で表現できる

 カウンタの値は最初は0にしておく
 for ( 20回繰り返す  ) {
   if (取り出したのが赤玉だったら) {
     カウンタの値を１増やす
   }
 }

カウンタのところは代入文で表現できる。for文も、Rで実行できる形式に直す

 カウンタ = 0
 for ( i in 1:20 ) {
   if (取り出したのが赤玉だったら) {
     カウンタ = カウンタ + 1
   }
 }

「もし、取り出したのが赤玉だったら」というif文をどうプログラミングすれば良いかが、今回のシミュレーションの最重要ポイントだろう。この部分は乱数を使って表現できる。

Untitled-10.gif
つまり、0-1の範囲の乱数を１つ発生させ、それが0.5よりも小さい値だったら、赤玉が取り出されたとみなせばいい。 乱数を１つ発生させるには runif(1) を使う。

 カウンタ = 0
 for ( i in 1:20 ) {
   if (runif(1) < 0.5 ) {
     カウンタ = カウンタ + 1
   }
 }
print(カウンタ)       #この命令で'カウンタ'という変数の内容を表示

これで20個取り出したうちの赤玉の個数を数えるという部分は完成です

Untitled-11.gif
上の図で、最初の矢印の部分の操作を行い、筒から取り出した赤玉の数を数えるシミュレーションのプログラミングは終わった。次に、取り出した赤玉の割合になるように筒の中の赤玉の数を変更して、同様の操作を100回繰り返すという部分のプログラミングを考える。繰り返しが含まれているのでfor文を使って大雑把に書くと、

 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   ※20個取り出して赤玉の数を数える
 }

たったこれだけ。上で作った赤玉の数を数えるプログラムを※のところに入れてみると

 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   カウンタ = 0
   for ( i in 1:20 ) {
     if (runif(1) < 0.5 ) {
       カウンタ = カウンタ + 1
     }
   }
  print(カウンタ)       #この命令で'カウンタ'という変数の内容を表示
 }

ん？？これでは上の図の例で、20個取り出した中に含まれていた赤玉の数の割合に、筒の中の赤玉の割合を変更したことが反映されていない。取り出した赤玉の数は、内側のfor文が終わった後に、カウンタの値を20で割ることで求められるから、「赤玉の割合」という変数を一つ作って、毎回、新しい値をこの変数に入れるようにしよう。

 赤玉の割合=0.5         #初期値は0.5
 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   カウンタ = 0
   for ( i in 1:20 ) {
     if (runif(1) < 赤玉の割合 ) {
       カウンタ = カウンタ + 1
     }
   }
   赤玉の割合 = カウンタ/20  #20個の中の赤玉の割合を算出
   print(カウンタ)       #この命令で'カウンタ'という変数の内容を表示
 }

あとは結果をプリントする文を入れれば完成。何回目に1色になったかが分かるように、j の値も一緒にプリントしよう。

 赤玉の割合=0.5         #初期値は0.5
 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   カウンタ = 0
   for ( i in 1:20 ) {
     if (runif(1) < 赤玉の割合 ) {
       カウンタ = カウンタ + 1
     }
   }
   赤玉の割合 = カウンタ/20  #20個の中の赤玉の割合を算出
   結果=c(j, 赤玉の割合)
   print(結果)
 }

できる人は、結果をグラフに表示してみよう。

実は、今作成したシミュレーションは、生物学科の人は皆知っている遺伝的浮動のシミュレーションだ。この例で、それぞれの変数は次のように対応している。

赤玉の割合：　遺伝子頻度の初期値
繰り返し数（ここでは100回）：　世代数
筒の中の玉の数 (ここでは20）：　集団サイズ

昨年度までは、最初から遺伝的浮動のシミュレーションを目指して授業を行ったのだが、今年度は話を少し単純にして、筒の中の玉の個数の変化という実験にしてみた。遺伝的浮動に興味がある人は、昨年度の授業ページを参照して欲しい。

↑

最後に †

一昨年、同様の授業を行ったところ、受講生から次のような質問があった。

「こういうシミュレーションのプログラミングって、これから生物学を学ぶのに、どのくらい必要なんですか？」

これまでプログラミングなどしたことの無い人にとって、とても素朴な疑問だろうと思う。そこで、次のように答えておいた。

画一的な答えは難しいです。シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。
ただ、私は、プログラミングは、生物学における発展的な情報処理をする上での基本技だと考えています。
千葉大の研究室に限って言えば、プログラミング技術が無いと卒業研究が終わらないというところは
現時点では（たぶん）無いでしょう。但し、プログラミングができると、卒業研究や修士の研究でより
進んだことができたり、理解が深まったり、研究が楽になるようなところはいくつもあります。
将来、数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、
敢えて、不可欠な技術だと言っておきます。
↑

第11回授業の課題 †

↑

課題１．意見調査 †

  1. http://bean.bio.chiba-u.jp/joho/index.php?joho20 に、「自分のID」/11　という新しいページを作成し、下の囲みの中にあるアンケートをコピー・ペーストして、「回答：」の後に答えを書き込むこと。
↑

課題2：（基本課題）コインを投げる実験のシミュレーション †

下のプログラムは、正常なコイン（表、裏の出る確率はそれぞれ0.5）を100回投げて、表の出た回数を記録するという実験を、1000回繰り返すシミュレーションのプログラムです。

omote=0                 #表の出る回数を初期化
kekka=c()
for(i in 1:【1】){
	for (i in 1:【2】){
   		if(runif(1)<【3】){
			omote=【4】+1
		}	
	}
	kekka=c(kekka,omote)
	omote=0                #表の出る回数を初期化
}
print(kekka)
↑

課題3:（発展課題）遺伝的浮動のシミュレーション †

下のプログラムは遺伝的浮動をシミュレートするRのプログラムです。このシミュレーションで観察する集団の数(num_repeats)、１集団あたり観察する世代数(num_generations)、集団サイズ（遺伝子数）(size_population)、０世代目における対立遺伝子aの数(num_a_allele)を与えて実行すれば、、ランダムに選ばれる遺伝子が、次世代の集団の遺伝子頻度にどのように影響するかを、シミュレートする集団数（観察集団数）の数だけ色つきの折れ線グラフで表示することができます。下のプログラムを見て以下の問に答えなさい。

drift= function(num_repeats,num_generations,size_population, num_a_allele){
	 results=c()
        a=num_a_allele
	 for(i in 1:【1】){
		  for(j in 1:【2】){
		      count_a=0
		       for(k in 1:size_population){
			    if ( runif(1) < a/【3】 ){
			     count_a=count_a+1
			    }
		       }
		       a=count_a
		       results=append(results, a/size_population)
		  }
		  a=num_a_allele
	 }
	 rmatrix=matrix(results, nrow=num_generations, ncol=num_repeats)
	 return(matplot(rmatrix, type="l"))
}
↑

参考：　円周率の計算を乱数を使わずに行うプログラム †

simpai=function(x){
	kaisu=0;  #個数の初期値を0にする
	for (i in 1:x) {
		for (j in 1:x) {
			point=c(i/x,j/x)
			distant=sqrt(point[1]^2+point[2]^2)
			if(distant<1){
			kaisu=kaisu+1
			}
		}
	}
	answer=4*kaisu/x^2
	print(answer)
}