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

はじめに

Rでプログラミングの授業を始めたころ、受講生から次のような質問があった。

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

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

シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。
プログラミングは、生物学における発展的な情報処理として、とても重要なものだと思います。
将来、数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、
必須技術といっても良いでしょう。

第11回授業の獲得目標: [worried]

  • 1. プログラミングの基本:繰り返し・代入・条件分岐に習熟する
  • 2. ユーザー定義関数を理解する
  • 3. シミュレーションの考え方を習得する
  • 4. シミュレーションのプログラミングを体験する
  • 5. 遺伝的浮動のシミュレーションに挑戦する!

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

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

繰り返し

  • 演習: 1から10までの数字を表示させるプログラムを作りなさい
    for(i in 1:__){    #{と}の間を10回繰り返す。i の値は毎回1ずつ増える
      print(i)
    }
    
    選択肢: a) 1    b)  5     c)  10

代入

  • 演習: 1から10までの数字の合計を計算するプログラムを、for文を使って作りなさい
    goukei=0           #goukeiという変数に初期値0を代入
    for(i in 1:10) {   
      goukei=_____+i  #goukeiに前回までのgoukeiの値にiの値を足したもの代入
    } 
    print(goukei)
    
    選択肢: a) i    b)  goukei   c)  print

条件分岐

  • 演習: 1から10までの数字のうち、2の倍数だけを合計するプログラムを、for文とif文を使って作りなさい
    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) <=

プログラミングTips

・1行に1つの命令を書いたら、必ず改行
・{  }  のそれぞれの後で改行すると見やすい(このページの例を参照)

おさらい: プログラミングは手順書き! 子供にお買い物を言葉で教えるようなもの

前の授業で、プログラミングは数学というよりも、コンピュータにことばで手順を伝えることだという説明をした。
例えば、おつかいに行く子供に、次のような買い物の手順を説明してみよう。

お財布の中に500円入っているから、これを持って買い物に行ってね。
この通りをまっすぐに進むとお店が3軒あるから、順番に買い物をするんだよ。
最初のお店で、100円のゴーヤを1本もらって、100円を渡すんだよ。
次のお店で100円のお豆腐を一丁もらって、100円を渡すんだよ。
次のお店で100g 100円の豚肉を、200gもらって、200円を渡すんだよ。
残ったお金はいくらだった?あとでお菓子をかってもいいわよ。

このお買い物の手順を、先週学んだ「繰り返し・代入・条件分岐」 を使って書いてみると、例えば、下のようなプログラムになる。

お財布の中 = 500
for (i in 1:3) {  #3回買い物を繰り返す
    if (i != 3) {    #3件目以外は100円ずつ使う
        お財布の中 = お財布の中 - 100
    } else {
         お財布の中 = お財布の中 - 200
    }
}
print(お財布の中)

プログラミングはこうやって、コンピュータに「手順」を伝えるものなのです。

  

おさらい:ユーザー定義関数 [smile]

「ユーザー定義関数」というのは、関数を定義する関数function()を使って、自分で作成したプログラムに名前を付けるようなものだった。オサライのため、1つ演習をしておこう。

  • 演習: #穴埋め式練習問題: _の部分に適当な文字を入力して実行しよう
    #入力した回数だけ乱数を1つ発生させて、それが0.5より小さいものだけを抜き出して表示させる
    ransu=___(a){   #ユーザー定義関数を定義
     kekka=c()     #結果をいれておく空ベクトルを作る
     for (i in 1:_) {  # a 回繰り返す
      x = __(1)           #乱数を1つ発生させる
      if (x  < 0.5 ) {   #発生させた乱数が0.5より小さいならば
         kekka=c(kekka, x)
      }
     }
     print(kekka)
    }

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

授業/H20/情報処理/R関数一覧  これまでに学んだRのいろんな関数を一覧できる。

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

これまでに学んだ繰り返し・代入・条件分岐の3つの命令それぞれは単純だが、組み合わせれば、かなり複雑なことでも表現できる。そこで、今回は、これらの命令を組み合わせて、Rを用いたシミュレーションに挑戦してみよう。
シミュレーションという言葉は「株式シミュレーション」とか「シミュレーションゲーム」とかで聞いたことがあるだろう。日本語にすると、「模擬実験」で、自分では実行が難しい繰り返し実験を、コンピュータに計算させることで行うものだ。コンピュータは単純な繰り返し計算は得意なので、砂粒をばらまいて、ある図形の中に入った数を数えるなんていう作業を、10万回でも、100万回でも繰り返すことができる。下の例では、そういうシミュレーションの方法を用いて、円周率(π)の値を求める。このように、シミュレーションをうまく使えば、たとえ数学的に解くのが難しい問題についても、近似解を求めることができる。また、ある現象について作成したモデルを用い、パラメーターの変化がどのような状態の変化をもたらすかを予測することができるため、進化生物学や生態学でもよく使われる。

今回は、乱数を用いたシミュレーションを行う。乱数とは、前回も勉強したが、ランダム(無作為)に得られる数値のことで、次にどんな数字が現れるのか予測できないような数字の集合のことを言う。Rでは、

runif()

という関数を使って、好きな数だけ乱数を得ることができる。
例えば、

runif(5)

と入力すれば、乱数が5個得られる。せっかくなので、得られた乱数をグラフにすることで、視覚的に乱数の性質を理解してみよう。

plot(runif(1000))   #枠内に、1000個の点がランダムにプロットされる
hist(runif(1000)    #発生させた1000個の乱数のヒストグラム; 1000個ぐらいだと、けっこうばらつく
 #数を増やすと、バラツキが小さくなるのが分かるだろうか?(これを大数の法則という)

前回のおさらい:ランダムウォークのシミュレーション

このシミュレーションは、「ランダムウォーク」で、次に移動する位置が確率的に無作為に決まるというものだ。
例えば、葉っぱの真ん中にポトリと落ちた虫が、全く無作為に葉っぱを食べ進むとき、食べ後はどのようなパターンを描くだろうか?
自分がすでに通った場所にも無作為に戻ると仮定して、ランダムウォークのシミュレーションを作ると、次のようになる。
ここで想定しているのは、x軸が-30から30、y軸が-30から30の正方形の空間で、原点からランダムウォークがスタートするとしたものだ。いくつか目新しい関数が入っているが、とりあえず無視して、Rにコピー・ペーストして、実行してみよう。

rwalk=function(x){
       frame() #グラフをクリア
	kaisu=0  #回数の初期値を0にする
	point=c(0,0) #原点を示すベクトルをpointという名前で作成
	move=floor(runif(x, 0, 4)) #0-3の範囲でx個の整数乱数を発生させる
       #floor(x)はx以上では無い最大の整数を返す関数
	for (i in 1:x) {
		if (move[i] == 0) { point[1] = point[1]+1 } #0なら右
		if (move[i] == 1) { point[2] = point[2]-1 } #1なら下
		if (move[i] == 2)	{ point[1] = point[1]-1 } #2なら左
		if (move[i] == 3)	{ point[2] = point[2]+1 } #3なら上
		if ((point[1] > 30) || (point[2] > 30) || (point[1] < -30) || (point[2] < -30)){
			break  #描画範囲の上限を超えたら止める
		}
              par(new=T)    #複数のグラフを重ね合わせる
              plot(point[1], point[2], ylim=c(-30,30), xlim=c(-30,30)) 
                                      #次の移動位置(pointの場所)をプロット
	}
}

実行するには、

> rwalk(100)

と入力する。動かす回数は()内の数値で指定する. 60 x 60の正方形の空間をはみ出すには、何回ぐらいかかるかを予想してからやってみると面白い。
また、ここにあげたのは非常にシンプルな条件でのランダムウォークのシミュレーションだが、虫が実際に葉っぱをたべる様子をシミュレートする場合、葉の縁から食べ始めるとか、一度たべたところには戻らないとか、自分が落ちてしまわないように食べるなどの条件を付け加える。シミュレーションのプログラミングは、単純なモデルからはじめ、徐々に複雑にして行けば良い。

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

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

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

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

正方形の中の砂粒の数(n) : 扇形の中砂粒の数(m)  = 
正方形の面積 :扇形の面積 =
1 : π/4

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

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

もしあなたが、砂粒をたくさんランダムにまいて、正方形の中の砂粒と扇形の中の砂粒の数を数えることができれば、πの値がわかるということになる。自分で実際に実験することを想像すると、恐ろしく大変だが、そういう面倒な実験はコンピュータにやらせてしまえばいい!

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

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

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

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

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

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

sqrt(0.230^2+0.782^2) < 1

となり、ランダムに発生させた(図形の上に偶然落ちた)1つの点は、扇形の中だったということになる。
実際に2つの乱数を発生させて、それを、pointというオブジェクトに代入するには、次のようにする。

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

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

point[1]
point[2]

のようにカギ括弧[]を使う。

例: []の使い方を、試してみよう
        > x=runif(2)
        > x
        > x[1]
        > x[2]

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

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

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

 #個数を数える部分
 kaisu=0   #回数を数えた結果を入れる場所なので、最初の値(初期値)は0にしておく
   if(distance<1){      #原点からの距離が1より小さければ扇形の中
       kaisu = kaisu + 1   #条件にあったとき、kaisuの値を1増やす
  }

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

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

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

simpai=function(x){
	kaisu=0;  #個数の初期値を0にする
	for (i in 1:x) {
		point=runif(2)
		distance=sqrt(point[1]^2+point[2]^2)
		if(distance<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

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

plot()関数でylimというオプションを使うと、グラフの軸の範囲を指定できる(詳しくはマニュアルを参照)

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の範囲で乱数の値をプロット
		distance=sqrt(point[1]^2+point[2]^2)
		if(distance<1){
			kaisu=kaisu+1
		}
	}
	answer=4*kaisu/x
	print(answer)
}

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

発展演習: 円周率の計算を乱数を使わずに行うプログラム

  • これまでの説明が分かった人は、次のような考え方でπを求めてみよう。
  • 上のプログラムとよく似ているが、乱数を使わず、1x1の正方形をマス目に分割して、扇形の中に入っている交点の数を数えるもの。穴埋め問題にしておくので、プログラムを完成させてほしい。得られる値は、100のとき3.1、1000のとき3.137524、10000のとき3.14119。100000のときは計算に時間がかかりすぎて途中でストップ。
simpai=function(x){
	kaisu=0;  #個数の初期値を0にする
	for (_ in 1:_) {
		for (_ in 1:_) {
			point=c(___,___)
			distance=sqrt(point[1]^2+point[2]^2)
			if(distance<1){
			 kaisu=kaisu+1
			}
		}
	}
	answer=4*kaisu/x^2
	print(answer)
}

遺伝的浮動のシミュレーションに挑戦する! [smile]

次に、生物学科の学生なら、入試のときに必ず勉強している遺伝的浮動のシミュレーションをやってみよう。
遺伝的浮動とは、

大きさの決まった集団で起きる、配偶子のランダムサンプリングによる遺伝子頻度の変動

のことだ。例えば、集団中の2つの対立遺伝子を、青玉と赤玉で表すと、袋(配偶子の集団: 配偶子の数はものすごく多く、青と赤の比率はそれぞれ50%)から取り出される青玉・赤玉の割合(集団中の遺伝子頻度)には、下の図に示したような変動が生じる(下の例で、青玉の遺伝子頻度は、0.8, 0.7, 1.0, と変動している)。もし、袋から取り出されたのが青玉ばっかり(あるいは、赤玉ばっかり)だと、それ以降、何度繰り返しても変化は生じない(下の例だと、「集団は青玉に固定した」などという)。
#11_1.gif

遺伝的浮動のシミュレーションをプログラミングして、赤玉と青玉の変動の度合いをグラフ化する

それでは、赤玉と青玉が同じ数だけたくさん入った袋から、10個の玉を取り出して、その中に含まれる青玉の割合に従って袋に含まれる青玉配偶子の個数を変え、また10個取り出す。。。こういう操作を何度も繰り返すと、取り出した玉のうちの青玉と赤玉の割合の変動がわかる。

自分で何度も実験を繰り返してもいいけれど、とても面倒!
そこで、この実験をコンピュータにやらせ、取り出した玉の割合がどのように変化するかをグラフにしよう!
  • でも、シミュレーションを始める前に、
    質問: この実験を何度も繰り返しておこなうとき、取り出した玉が赤か青かどちらか1色になってしまうのは、
    何回目だと思いますか?
    (袋の中の青玉と赤玉の数は同数からはじめ、 袋から取り出すのは10個)

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

もう皆さんは繰り返し、代入、条件分岐というプログラミングの基本技を修得しているので、今問題にしている実験を、この3つのワザを使ってどのように表現するかが、プログラミングの重要なポイントだ。

  1. まず袋から玉を10個、ランダムに取り出すところを考える

    #11_2.gif

    「赤玉と青玉の同じ数入っている袋から無作為に10個の玉を取り出したら、青玉が7個入っていた」
      -----------------------------------------------------
    上でやった操作には、「10個取り出して青玉かどうか判断する」というところに「繰り返し」が含まれている。
    「繰り返し」がはっきりと分かるように説明すると、
     1. 袋から玉を一つ取り出す
     2. その玉が青玉だったら数える(1つ目だったら「一つ」、2つめは「2つ」..)。赤玉だったら数えない
     3.  1 と 2 の操作を10回繰り返す
    この部分を大雑把でいいので、プログラミングしてみる。繰り返しはfor文で表現できる
     カウンタの値は最初は0にしておく
     for ( 10回繰り返す  ) {
       if (取り出したのが青玉だったら) {
         カウンタの値を1増やす
       }
     }
    カウンタのところは代入文で表現できる。for文も、Rで実行できる形式に直す
     カウンタ = 0
     for ( i in 1:10 ) {
       if (取り出したのが青玉だったら) {
         カウンタ = カウンタ + 1
       }
     }
  2. 次に、「もし、取り出したのが青玉だったら」という部分を、if文を使ってプログラミングする
    この点が、シミュレーションの最重要ポイントだ。
    ここでは、前回の予習課題で使った乱数であるrunif()という関数を使う。

    #11_3.gif
    つまり、

    0-1の範囲の乱数を1つ発生させ、それが0.5よりも小さい値だったら、
    青玉が取り出されたとみなす
    あるいは、こんな説明でも良いかもしれない
    例えば、袋の中に100個の玉が入っているとき、
    青には1-50の、赤には51-100の数字を書いておく
    袋から玉を一個とりだして、書いてある数字が50以下なら、青としてカウントする
     数字: 1 - 50 : 青  → 乱数: 0-0.5
     数字: 51-100:  赤  → 乱数: 0.5-1
    乱数を1つ発生させるには runif(1) を使って、上のプログラム「取り出したのが青玉だったら」の部分を書き換える
     カウンタ = 0
     for ( i in 1:10 ) {
       if (runif(1) < 0.5 ) {
         カウンタ = カウンタ + 1
       }
     }
    print(カウンタ)       #この命令で'カウンタ'という変数の内容を表示
    これで10個取り出したうちの青玉の個数を数えるという部分は完成!
    Rのスクリプトエディタにコピーして、実行してみよう。何度も実行すると、値は変化するのが分かる。
  3. 上の実験を100回繰り返して、遺伝子頻度の変動を見る

    #11_2.gif
    上の図で黄色く塗った部分:「袋から取り出した青玉の数を数えるシミュレーションのプログラミング」は終わった。次に、取り出した青玉の割合になるように袋の中の配偶子の割合を変更して、同様の操作を100回繰り返すという部分のプログラミングを考える。上の図で青く塗った囲みの部分がそれに相当する。

    繰り返しが含まれているのでfor文を使って書けることは分かるだろう

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

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

     for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
       カウンタ = 0
       for ( i in 1:10 ) {
         if (runif(1) < 0.5 ) {
           カウンタ = カウンタ + 1
         }
       }
      print(カウンタ)       #カウンタの内容を表示
     }
    実際に走らせてみると、、、、???ん?何度やっても0や10に固定しない。
     ??どこに問題があるか分かるだろうか??
    問題は、上の図のように、「10個取り出した中の青玉の数を数え、その割合に袋の中の青玉の割合を変更する」ということが反映されていないからだ。取り出した青玉の数は、内側のfor文が終わった後のカウンタの値に入っているので、その頻度を得るには、カウンタを10で割ることで求められる。
    そこで、「青玉の頻度」という変数を一つ作って、毎回、新しい値をこの変数に入れるようにしよう。
     青玉の頻度=0.5         #初期値は0.5
     for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
       カウンタ = 0
       for ( i in 1:10 ) {
         if (runif(1) < 青玉の頻度 ) {
           カウンタ = カウンタ + 1
         }
       }
       青玉の頻度 = カウンタ/10  #10個の中の青玉の割合を算出
       print(カウンタ)       #この命令で'カウンタ'という変数の内容を表示
     }
    これを実行すると、何回目かには0か1に固定することがわかる。
    これで、遺伝的浮動のシミュレーションは概ね完成した!!
    あとは結果をプリントする文を入れれば完成。何回目に1色になったかが分かるように、j の値も一緒にプリントしてみよう。
     青玉の割合=0.5         #初期値は0.5
     結果=c()      #結果を入れる空ベクトルを準備しておく
     for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
       カウンタ = 0
       for ( i in 1:10 ) {
         if (runif(1) < 青玉の割合 ) {
           カウンタ = カウンタ + 1
         }
       }
       青玉の割合 = カウンタ/10  #10個の中の青玉の割合を算出
       結果=c(j, 青玉の割合)  #何回目かという値はjに入っている。青玉の割合と一緒にベクトルに入れる
       print(結果)
     }
  4. では、最後に、結果をグラフにして表示してみよう
    グラフにするにはplot()という関数を前に使った。グラフを作るには、カウンタの値をベクトルに保存して
    (0.65, 0.7, 0.75, 0.85, ....)
    という形式にしておく必要がある。そこで、最初に空ベクトルを1つ作っておき、カウンタの値をベクトルに保存し、最後にplot()でグラフにする。穴埋めにしておくので、完成させて、実行してほしい。
     結果=c()    #空ベクトル
     青玉の頻度=0.5         #初期値は0.5
     for ( _ in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
       カウンタ = 0
       for ( _ in 1:10 ) {
         if (___(1) < 青玉の頻度 ) {
           カウンタ = カウンタ + 1
         }
       }
       青玉の頻度 = カウンタ/10  #10個の中の青玉の割合を算出
       結果=c(結果,青玉の頻度)       #青玉の頻度を結果というベクトルに追加
     }
     plot(結果)   #最後に結果をグラフにプロット
  5. グラフをもっとキレイにする
    これで遺伝的浮動のシミュレーションの結果をグラフにすることができたが、グラフは点々で表示されるだけだ。そこで、教科書に載っているような、きれいな折れ線グラフにしてみる。上のplot()のところを、下のように書き換えるだけで、折れ線グラフになる(type="l"というのが、折れ線グラフで表示するというオプションの指定)。
     plot(結果, type="l")   #type="l" が折れ線グラフにするオプション
    に置き換えるだけ。
    これで折れ線グラフにはなったが、Y軸が固定されていないのが気になってしまう。そこで、Y軸のメモリを0から1に固定しよう。
    y軸の目盛りを0から1に固定するには、plot()関数のylimというオプションを使う

     結果=c()    #空ベクトル
     青玉の頻度=0.5         #初期値は0.5
     for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
       カウンタ = 0
       for ( i in 1:10 ) {
         if (runif(1) < 青玉の頻度 ) {
           カウンタ = カウンタ + 1
         }
       }
       青玉の頻度 = カウンタ/10  #10個の中の青玉の割合を算出
       結果=c(結果,青玉の頻度)       #青玉の頻度を結果というベクトルに追加
     }
    plot(結果, type="l", ____=c(0,1))   #★★ylimを使ってy軸の目盛りを0から1に固定する

実行してみると、うまくいってるようだが、1回ずつ、グラフが変わってしまう。グラフを重ねて描くには、par(new=T) というコマンドをplot()の前に使う。もしいくつもグラフを重ねて描いて、グラフが汚くなってしまったら、

frame()
か
plot.new()

というコマンドを使ってグラフをクリアする。

 結果=c()    #空ベクトル
 青玉の頻度=0.5         #初期値は0.5
 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   カウンタ = 0
   for ( i in 1:10 ) {
     if (runif(1) < 青玉の頻度 ) {
       カウンタ = カウンタ + 1
     }
   }
   青玉の頻度 = カウンタ/10  #10個の中の青玉の割合を算出
   結果=c(結果,青玉の頻度)       #青玉の頻度を結果というベクトルに追加
 }
par(new=T)   #グラフを一度クリアする
plot(結果, type="l", ____=c(0,1))   #★★ylimを使ってy軸の目盛りを0から1に固定する


  1. さらに、ユーザー定義関数にして、グラフを重ねて書いてみよう
    なんども実行してみると、グラフがどんどん重なって書き換わるのが分かる。1回のキー操作が、1回の実験に相当している。キー操作で書き換えるのは、面倒なので、最初から実験回数を入力して、その回数だけ自動的に繰り返すような関数を作ってみよう。
    注意: 1回の実験が終わる度に、青玉の頻度を0.5に、結果のベクトルを空にリセットすることを忘れないこと
    drift=function(x){
       青玉の初期値=0.5         #初期値は0.5
     for(k in 1:x){              #関数に与えた値の数だけ実験を繰り返す
       青玉の頻度=青玉の初期値
       結果=c()    #空ベクトル
       for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
         カウンタ = 0
         for ( i in 1:10 ) {
           if (runif(1) < 青玉の頻度 ) {
             カウンタ = カウンタ + 1
           }
         }
         青玉の頻度 = カウンタ/10  #10個の中の青玉の割合を算出
         結果=c(結果,青玉の頻度)       #青玉の頻度を結果というベクトルに追加
       }
      par(new=T)   #グラフを一度クリアする
      plot(結果, type="l", ylim=c(0,1))   #★★ylimを使ってy軸の目盛りを0から1に固定する
     }
    }

発展:様々な変数を指定できる関数を作ってしまう

上のシミュレーションは集団の遺伝子の数が10の場合だけのシミュレーションだった。そこで、遺伝子の数(集団サイズ)、遺伝子頻度の初期値、繰り返し世代数、繰り返し実験数の全てを指定できるような関数を作ってみよう。複数の引数は、function()の中でカンマを使って指定できる

drift=function(実験回数, 初期頻度,集団サイズ,世代数){
   frame()     #グラフをクリアする
 for(k in 1:実験回数){              #関数に与えた実験回数だけ実験を繰り返す
   青玉の頻度=初期頻度
   結果=c()    #空ベクトル
   for ( j in 1:世代数 ) {   #さっきiという変数を使ったので、ここではjを使う
     カウンタ = 0
     for ( i in 1:集団サイズ ) {
       if (runif(1) < 青玉の頻度 ) {
         カウンタ = カウンタ + 1
       }
     }
     青玉の頻度 = カウンタ/集団サイズ  #集団サイズ個の中の青玉の割合を算出
     結果=c(結果,青玉の頻度)       #青玉の頻度を結果というベクトルに追加
   }
  par(new=T)   
  plot(結果, type="l", ylim=c(0,1))   #★★ylimを使ってy軸の目盛りを0から1に固定する
 }
}

発展・演習:

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

drift= function(num_repeats,num_generations,size_population, num_a_allele){
	 results=c()
        a=num_a_allele
	 for(【1】){
		  for(【2】){
		      count_a=0
		       for(k in 1:size_population){
			    if (【3】){
			     count_a=count_a+1
			    }
		       }
		       a=count_a
		       results=append(results, 【4】)
		  }
		  a=num_a_allele
	 }
	 rmatrix=matrix(results, nrow=num_generations, ncol=num_repeats)
	 return(matplot(rmatrix, type="l"))
}
  • 問1: プログラムが正常に動くよう、【1】〜【4】を補いなさい。
    • なお、プログラム中、それぞれの変数は次のような意味を持っています
      • num_repeats : シミュレーションを行う集団数(1つの集団についてある世代数だけ観察するという実験を、何回(何個の集団について)おこなうか)、すなわち観察集団数。何本の折れ線グラフを表示させるかはこの値で決まる。
      • num_generations : 1つの集団について何世代続けて観察するか。x軸の長さになります。
      • size_population : 集団のサイズ、つまり、1つの集団に含まれる遺伝子の総数
      • num_a_allele : aアリルの数(つまり、この値をsize_populationで割ったものが、0世代目でのaアリルの遺伝子頻度になる)
  • 問2: 上で作成したdrift()というユーザー定義関数を、次の2つの場合で実行しなさい。
    1. 実験繰返し10、1集団での観察世代数500、集団サイズ(遺伝子数)200, 対立遺伝子aの数の初期値100
    2. 実験繰返し10、1集団での観察世代数100、集団サイズ(遺伝子数)20, 対立遺伝子aの数の初期値10
    いずれの場合についても、できたグラフをWinShotで切り取って、課題提出ページに貼り付けなさい(注:1番目のシミュレーションの実行には、数分かかるかもしれません)(注:あまり大きい図(画面全面に広がる図とか)にはしないでください)
    ''上の取り消し線部分は、この演習問題を以前、課題として使っていたときの名残でした。消し忘れており、すみませんでした。今回、このシミュレーションの結果を提出する必要はありません。moodleページにある課題のみ、提出してください。
  • 問3: 上の2つのグラフを比較して分かることを述べなさい(集団遺伝学の問題)。

第11回授業の課題

  • 提出期限:7月3日水曜正午

Last-modified: 2015-05-13 (水) 16:45:42 (3269d)