*Rを使ったプログラミング演習2: シミュレーション [#l740a591] #contents **はじめに [#d6278ce5] Rでシミュレーションの授業を始めたころ、受講生から次のような質問があった。 「シミュレーションのプログラミングって、これから生物学を学ぶのに、どのくらい必要なんですか?」 これまでプログラミングなどしたことの無い人にとって、とても素朴な疑問だろうと思う。私からの答えとしては、 シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。 プログラミングは、生物学における発展的な情報処理として、とても重要なものだと思います。 将来、数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、 必須技術といっても良いでしょう。 ということで、今回は、Rを使ったシミュレーションに挑戦してみよう。 **第11回授業の獲得目標:&worried; [#rf3df5e5] -''1. プログラミングの基本:繰り返し・代入・条件分岐に習熟する'' -''2. ユーザー定義関数を理解する'' -''3. シミュレーションの考え方を習得する'' -''4. シミュレーションのプログラミングを体験する'' -''5. 遺伝的浮動のシミュレーションに挑戦する!'' **まずはおさらい1: プログラミングの基本: 繰り返し・代入・条件分岐 ⌣ [#wc8dcad4] 次のそれぞれのプログラムの_部分に何を入れるべきかを選択肢から選びなさい。また、Rを使って実行してみよう。 ***繰り返し [#mcc6a699] -演習: 1から10までの数字を表示させるプログラムを作りなさい for(i in 1:__){ #{と}の間を10回繰り返す。i の値は毎回1ずつ増える print(i) } 選択肢: a) 1 b) 5 c) 10 簡単にできた人は、結果を横一列に表示させてみよう。 ***代入 [#n018f210] -演習: 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 ***条件分岐 [#x42b9425] -演習: 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 [#s50684f1] ・1行に1つの命令を書いたら、必ず改行 ・{ } のそれぞれの後で改行すると見やすい(このページの例を参照) **おさらい2: プログラミングは手順書き! [#md800d1e] 前の授業で、&size(16){「''プログラミングは数学というよりも、コンピュータにことばで手順を伝えることだ''」};という説明をした。~ 例えば、おつかいに行く子供に、次のような買い物の手順を説明してみよう。 お財布の中に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(お財布の中) このプログラムの中に含まれている「繰り返し・代入・条件分岐」 の働きを、十分に理解しておこう。 **おさらい3:ユーザー定義関数⌣ [#wb2372b3] 「ユーザー定義関数」というのは、関数を定義する関数''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) } ransu(50) などと入力すれば、0.5より小さい乱数がいくつか表示される **おさらい4. 授業で使ったRの基本関数 [#m920a376] //#include(授業/20/情報処理/R関数一覧,notitle) →[[授業/H20/情報処理/R関数一覧]] これまでに学んだRのいろんな関数を一覧できる。 **Rを用いたシミュレーション⌣ [#fe9ffeb6] これまでに学んだ繰り返し・代入・条件分岐の3つの命令それぞれは単純だが、組み合わせれば、かなり複雑で、研究に役立ちそうなことでも表現できる。そこで、今回は、これらの命令を組み合わせて、''Rを用いたシミュレーション''に挑戦してみよう。~ シミュレーションという言葉は「株式シミュレーション」とか「シミュレーションゲーム」とかで聞いたことがあるだろう。日本語にすると、「模擬実験」で、自分では実行が難しい繰り返し実験を、コンピュータに計算させることで行うものだ。コンピュータは単純な繰り返し計算は得意なので、砂粒をばらまいて、ある図形の中に入った数を数えるなんていう作業を、10万回でも、100万回でも繰り返すことができる。下の例では、実際に、ある図形の上に砂粒をまく実験をコンピュータに繰り返し行わせるシミュレーションを行って、円周率(π)の値を求める。~ この例ように、シミュレーションをうまく使えば、たとえ数学的に解くのが難しい問題についても、近似解を求めることができる。また、ある生命現象について、自分で作成したモデルを用い、パラメーターの変化がどのような状態の変化をもたらすかを予測することができるため、進化生物学や生態学でもよく使われる。2つめのシミュレーションでは、遺伝的浮動のシミュレーションをやってみよう~ ***円周率(π)を「繰り返し」・「乱数」・「条件分岐」で求める [#s0ee7a6e] &ref(授業/H20/情報処理/11/#10_1.jpg,60%);~ 左に1辺の長さが1の正方形と、それに内接する扇形がある。このとき、 正方形の面積は 1 扇形の面積は π/4 この図形の上に、砂粒をばらまいてみる。 砂粒が図形の上にランダムに散らばるとすると、砂粒をすごく沢山まいて数を数えれば、 正方形の中の砂粒の数(n) : 扇形の中砂粒の数(m) = 正方形の面積 :扇形の面積 = 1 : π/4 になると考えられる。つまり n : m = 1 : π/4 だから、 π = 4m/n もしあなたが、砂粒をたくさんランダムにまいて、正方形の中の砂粒と扇形の中の砂粒の数を数えることができれば、πの値がわかるということになる。自分で実際に実験することを想像すると、恐ろしく大変だが、そういう面倒な実験はコンピュータにやらせてしまえばいい! ***シミュレーションの考え方 [#xf7b0c77] -''砂粒をたくさんまく部分''~ さて、砂粒を沢山まいて、それが円の中か外かを判定する実験を、コンピュータにやらせるにはどうすればよいだろうか。 まず、点(砂粒)を大量に発生させる必要がある。さて、どうすればいいだろうか?~ もうお分かりのように、「繰り返し」命令(for文)を使う。 #とりあえず、1,000つぶの砂をまくことをfor文で表現する for (i in 1:____){ print(i) #<砂を1粒図形の上にまく> } -''砂粒をランダムにまく部分''~ 次に、砂を1粒図形の上にまいて、それが扇形の中か外かを判定することを考える。~ 砂を1粒まくことは、無作為に(ランダムに)図形の上に点を打つことと同じことだ。先週学んだ、''runif()'' という乱数を発生させる関数を使えば、%%%図形の上に1つの点を打つことは、2つの乱数を発生させることで表現できるだろう%%%。 &ref(授業/H20/情報処理/11/#10_2.jpg,60%);~ 例えば、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 //http://www.bunkyo.ac.jp/~nemoto/lecture/simulation/97/random/ppframe.htm ***せっかくだから、点をうつところを図で表現してみよう [#o33bd569] ''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) } &ref(授業/H20/情報処理/11/simpai100000.jpg,60%);~ 10万回の実行結果。 **発展演習: 円周率の計算を乱数を使わずに行うプログラム [#ucb45d16] -これまでの説明が分かった人は、次のような考え方でπを求めてみよう。 -上のプログラムとよく似ているが、乱数を使わず、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) } **遺伝的浮動のシミュレーションに挑戦する!⌣ [#a2a76046] 次に、生物学科の学生なら、大学入試にそなえて必ず勉強した''遺伝的浮動''のシミュレーションをやってみよう。~ 遺伝的浮動とは、 大きさの決まった集団で起きる、配偶子のランダムサンプリングによる遺伝子頻度の変動 のことだ。例えば、集団中の2つの対立遺伝子を、青玉と赤玉で表すと、袋(配偶子の集団: 配偶子の数はものすごく多く、青と赤の比率はそれぞれ50%)から取り出される青玉・赤玉の割合(集団中の遺伝子頻度)には、下の図に示したような変動が生じる(下の例で、青玉の遺伝子頻度は、0.8, 0.7, 1.0, と変動している)。もし、袋から取り出されたのが青玉ばっかり(あるいは、赤玉ばっかり)だと、それ以降、何度繰り返しても変化は生じない(下の例だと、「集団は青玉に固定した」などという)。~ &ref(授業/H21/情報処理/11/#11_1.gif,50%); ***出席確認演習+グループワーク: 0か1に固定するまでの世代数 [#x99dae1d] まず、上のような袋と玉をつかった実験を手作業でやった場合、取り出した玉が赤か青かどちらか1色になってしまうのは、何回目だろうか?下の条件それぞれで、グループで相談してみよう。また、相談後に自分が考えてた固定までの回数を、授業Moodleページの出席確認演習に記入しよう。 -最初、袋の中の青玉と赤玉の数は同数で、袋から取り出すのは10個 -最初、袋の中の青玉と赤玉の数は同数で、袋から取り出すのは100個 -最初、袋の中の青玉と赤玉の数は同数で、袋から取り出すのは1000個 ***遺伝的浮動のシミュレーションの学び方 [#zb1eb43c] では、上で説明した実験を、これまでに習った''繰り返し、代入、条件分岐''というプログラミングの基本技を使って表現してみよう。~ ここで、実験をどうやって表現するするかが、プログラミングの重要なポイントだ。この解説では、まず最初に完成版ではあるが、部分的に抜け落ちのあるプログラムをお見せする。このプログラムを見て、シミュレーションの働きを考え、完成させて、実際にシミュレーションを行ってみよう。 ***演習・グループワーク: 遺伝的浮動のシミュレーション: 完成版(ただし、一部、抜け落ち) [#o8ffec18] 下のプログラムは遺伝的浮動をシミュレートするRのプログラムだが、一部、抜け落ちている。また、このプログラムでは、頭の中での置き換えを簡単にするために、オブジェクト(変数)名が、日本語で表示されている。 -''観察集団数'': このシミュレーションで観察する独立した集団の数 -''観察世代数'': 1集団について観察する世代数 -''玉の総数'': 取り出した玉の数。集団内の遺伝子数に相当) -''青玉の初期数'': 0世代目の玉の総数中の青玉の数. 集団内の対立遺伝子aの数に相当) これらの値を与えて実行すれば、ランダムに選ばれた遺伝子が、次世代の集団の遺伝子頻度にどのように影響するかを、シミュレートする観察集団数の数だけ、色つきの折れ線グラフで表示することができる。~ -問1. 下のプログラムの下線部には、観察集団数、観察世代数, 玉の総数, 青玉初期数 のどれかが入る。どれが、どこに入るかを、グループで相談してみよう。 -問2. 上の出席確認演習でやった、3つの条件のシミュレーションを、10集団について、100世代、観察してみよう。 --玉の数10個の場合: drift(10,100,10,5) --玉の数100個の場合: drift(10,100,100,50) --玉の数1000個の場合: <自分で考えること> -問3. [[2008年度前期試験に出た問題>http://hiw.oo.kawai-juku.ac.jp/nyushi/honshi/08/ch1.html]]のシミュレーションをやってみよう。 集団:1万人(遺伝子の集団サイズ20000)、Hbsの遺伝子頻度、約0.1 世代数 100 地球上からマラリアが撲滅され、鎌状赤血球貧血症で死ぬことも無くなった場合 drift = function(観察集団数,観察世代数,玉の総数,青玉初期数){ 結果 = c() #結果を入れる空ベクトルを準備 for( i in 1 : ________ ){ 青玉の数 = ________ #青玉の数は観察集団が代わる度に初期値に戻す for( j in 1 : ________ ){ 結果 = append(結果, 青玉の数 / 玉の総数) #青玉の頻度を「結果」ベクトルに入れる #最初の繰り返しでは青玉の初期頻度が入る 青玉計数 = 0 #袋の中の青玉の数を青玉計数に入れる。初期値は0 for( k in 1 : ________ ){ if (runif(1) < 青玉の数 / 玉の総数 ){ #「青玉の数/玉の総数」は青玉の頻度 青玉計数 = 青玉計数 + 1 #青玉計数の値を1つ増やす(青玉を数えるということ) } } 青玉の数 = 青玉計数 #次の世代の袋の中の「青玉の数」に、今数えた「青玉計数」を入れる } } 結果の表 = matrix(結果, nrow = 観察世代数, ncol = 観察集団数) print(結果の表) return(matplot(結果の表, type="l", ylim=c(0,1))) } ***演習問題の答えと解説 [#m466ca08] drift = function(観察集団数,観察世代数,玉の総数,青玉初期数){ 結果 = c() #結果を入れる空ベクトルを準備 for( i in 1 : 観察集団数 ){ 青玉の数 = 青玉初期数 #青玉の数は観察集団が代わる度に元に戻す for( j in 1 : 観察世代数 ){ 結果 = append(結果, 青玉の数 / 玉の総数) #青玉の頻度を「結果」ベクトルに入れる #最初の繰り返しでは青玉の初期頻度が入る 青玉計数 = 0 #袋の中の青玉の数を青玉計数に入れる。初期値は0 for( k in 1 : 玉の総数 ){ if (runif(1) < 青玉の数 / 玉の総数 ){ #「青玉の数/玉の総数」は青玉の頻度 青玉計数 = 青玉計数 + 1 #青玉計数の値を1つ増やす(青玉を数えるということ) } } 青玉の数 = 青玉計数 #次の世代の袋の中の「青玉の数」に、今数えた「青玉計数」を入れる } } 結果の表 = matrix(結果, nrow = 観察世代数, ncol = 観察集団数) print(結果の表) return(matplot(結果の表, type="l", ylim=c(0,1))) } +'' まず袋から玉を10個、ランダムに取り出すところを考える''~ ~&ref(授業/H21/情報処理/11/#11_2.gif,50%); 「赤玉と青玉の同じ数入っている袋から無作為に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 } } +''次に、「もし、取り出したのが青玉だったら」という部分を、if文を使ってプログラミングする''~ この点が、シミュレーションの最重要ポイントだ。~ ここでは、前回の予習課題で使った''乱数''である''runif()''という関数を使う。 ~&ref(授業/H21/情報処理/11/#11_3.gif,50%);~ つまり、 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(カウンタ) #この命令で'カウンタ'という変数の内容を表示 &size(16){これで10個取り出したうちの青玉の個数を数えるという部分は完成!};~ Rのスクリプトエディタにコピーして、実行してみよう。何度も実行すると、値は変化するのが分かる。 +''上の実験を100回繰り返して、遺伝子頻度の変動を見る''~ ~&ref(授業/H21/情報処理/11/#11_2.gif,50%);~ 上の図で黄色く塗った部分:「袋から取り出した青玉の数を数えるシミュレーションのプログラミング」は終わった。次に、取り出した青玉の割合になるように袋の中の配偶子の割合を変更して、同様の操作を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に固定することがわかる。~ &size(16){''これで、遺伝的浮動のシミュレーションは概ね完成した!!''};~ あとは結果をプリントする文を入れれば完成。何回目に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(結果) } +''では、最後に、結果をグラフにして表示してみよう''~ グラフにするには''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(結果) #最後に結果をグラフにプロット +''グラフをもっとキレイにする''~ これで遺伝的浮動のシミュレーションの結果をグラフにすることができたが、グラフは点々で表示されるだけだ。そこで、教科書に載っているような、きれいな折れ線グラフにしてみる。上の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回の実験が終わる度に、青玉の頻度を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に固定する } } ***発展:様々な変数を指定できる関数を作ってしまう [#r3d25287] 上のシミュレーションは集団の遺伝子の数が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に固定する } } ***発展・演習: [#p16201fc] -問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番目のシミュレーションの実行には、数分かかるかもしれません)(注:あまり大きい図(画面全面に広がる図とか)にはしないでください)%%~ (&color(red){''上の取り消し線部分は、この演習問題を以前、課題として使っていたときの名残でした。消し忘れており、すみませんでした。今回、このシミュレーションの結果を提出する必要はありません。[[moodleページ>http://bean.bio.chiba-u.jp/moodle24]]にある課題のみ、提出してください。}; -問3: 上の2つのグラフを比較して分かることを述べなさい(集団遺伝学の問題)。 **第11回授業の課題 [#tccc92c9] -提出期限:''7月3日水曜正午''