![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
はじめに
Rでプログラミングの授業を始めたころ、受講生から次のような質問があった。
「シミュレーションのプログラミングって、これから生物学を学ぶのに、どのくらい必要なんですか？」
これまでプログラミングなどしたことの無い人にとって、とても素朴な疑問だろうと思う。私からの答えとしては、
シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。 プログラミングは、生物学における発展的な情報処理として、重要なものだと思います。 将来、数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、 きっと必要になるでしょう。
![[worried]](image/face/worried.png)
†![[smile]](image/face/smile.png)
†次のそれぞれのプログラムの＿部分に何を入れるべきかを選択肢から選びなさい。また、Rを使って実行してみよう。
for(i in 1:__){ #{と}の間を10回繰り返す。i の値は毎回1ずつ増える
print(i)
}
選択肢： a) 1 b) 5 c) 10goukei=0 #goukeiという変数に初期値0を代入
for(i in 1:10) {
goukei=_____+i #goukeiに前回までのgoukeiの値にiの値を足したもの代入
}
print(goukei)
選択肢： a) i b) goukei c) printgoukei=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) <=・１行に１つの命令を書いたら、必ず改行
・{ } のそれぞれの後で改行すると見やすい（このページの例を参照）
前の授業で、「プログラミングは数学というよりも、コンピュータにことばで手順を伝えることだ」という説明をした。
例えば、おつかいに行く子供に、次のような買い物の手順を説明してみよう。
お財布の中に500円入っているから、これを持って買い物に行ってね。 この通りをまっすぐに進むとお店が３軒あるから、順番に買い物をするんだよ。 最初のお店で、100円のゴーヤを１本もらって、100円を渡すんだよ。 次のお店で100円のお豆腐を一丁もらって、100円を渡すんだよ。 次のお店で100g 100円の豚肉を、200gもらって、200円を渡すんだよ。 残ったお金はいくらだった？あとでお菓子をかってもいいわよ。
このお買い物の手順を、先週学んだ「繰り返し・代入・条件分岐」 を使って書いてみると、例えば、下のようなプログラムになる。
お財布の中 = 500
for (i in 1:3) { #３回買い物を繰り返す
if (i != 3) { #3件目以外は100円ずつ使う
お財布の中 = お財布の中 - 100
} else {
お財布の中 = お財布の中 - 200
}
}
print(お財布の中)
プログラミングが「手順」を伝えるものであることが、分かって貰えただろうか？
†「ユーザー定義関数」なんていうと難しそうだが、ようするに、先ほどまでに作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにしようというものだ。こういうときに関数を定義する関数、function()を使う。 例えば、円の面積を計算する関数を作るなら、
menseki=function(r){r*r*pi}
とすれば良い。これで自分独自の関数menseki()ができた。 では、実行してみよう。
menseki(10) #好きな数字を入れて円の面積を計算
演習 1： 1から入力した数値までの全てを横一列に表示させるプログラムを作りdisplayという名前の関数として定義しなさい
不足している部分を補充して、Rで実行すること。
display=____(a){ #関数定義の始まり
kekka=c() #kekkaに空ベクトルを代入して初期化
for(i in 1:a){ #a回（iの値を1からaまで変化させる間）繰り返し
kekka=c(kekka,i) #kekkaというベクトルにiを要素として代入
}
print(kekka) #kekkaの内容を表示
} #関数定義の終わり
↑を使って何回も実行してみると、結果がいろいろ変わるのがわかる。例えば10を入れて実行するには、
display(10)
演習2: 上の関数の定義方法に従って、入力した数までの合計値を計算するsumupという名前の関数を作成する。下の囲みの中の_の部分（１文字に対応するとは限らない）を埋めて、プログラムを完成しなさい。
sumup = ______ #関数sumpuを定義
kotae = _ #kotaeを初期化
_ (i in ___){ #1からaまで繰り返し
kotae = _____ #kotaeにiの値を足したものとkotaeに代入
} #繰り返し終了
_____ #kotaeを表示
} #関数定義の終了
→授業/H20/情報処理/R関数一覧 これまでに学んだRのいろんな関数を一覧できる。
†これまでに学んだ繰り返し・代入・条件分岐の３つの命令それぞれは単純だが、組み合わせれば、かなり複雑なことでも表現できる。そこで、今回は、これらの命令を組み合わせて、Rを用いたシミュレーションに挑戦してみよう。
シミュレーションという言葉は「株式シミュレーション」とか「シミュレーションゲーム」とかで聞いたことがあるだろう。日本語にすると、「模擬実験」で、自分では実行が難しい繰り返し実験を、コンピュータに計算させることで行うものだ。コンピュータは単純な繰り返し計算は得意なので、砂粒をばらまいて、ある図形の中に入った数を数えるなんていう作業を、10万回でも、100万回でも繰り返すことができる。下の例では、そういうシミュレーションの方法を用いて、円周率（π）の値を求める。このように、シミュレーションをうまく使えば、たとえ数学的に解くのが難しい問題についても、近似解を求めることができる。また、ある現象について作成したモデルを用い、パラメーターの変化がどのような状態の変化をもたらすかを予測することができるため、進化生物学や生態学でもよく使われる。
今回は、乱数を用いたシミュレーションを行う。乱数とは、ランダム（無作為）に得られる数値のことで、次にどんな数字が現れるのか予測できないような数字の集合のことを言う。Ｒでは、
runif()
という関数を使って、好きな数だけ乱数を得ることができる。
例えば、
runif(5)
と入力すれば、乱数が５個得られる。せっかくなので、得られた乱数をグラフにすることで、視覚的に乱数の性質を理解してみよう。
plot(runif(1000)) #枠内に、1000個の点がランダムにプロットされる hist(runif(1000) #発生させた1000個の乱数のヒストグラム； 1000個ぐらいだと、けっこうばらつく #数を増やすと、バラツキが小さくなるのが分かるだろうか？（これを大数の法則という）
左に１辺の長さが1の正方形と、それに内接する扇形がある。このとき、
正方形の面積は 1 扇形の面積は π/4
この図形の上に、砂粒をばらまいてみる。 砂粒が図形の上にランダムに散らばるとすると、砂粒をすごく沢山まいて数を数えれば、
正方形の中の砂粒の数(n) ： 扇形の中砂粒の数(m) ＝ 正方形の面積 ：扇形の面積 ＝ １ ： π/4
になると考えられる。つまり
n : m = 1 : π/4 だから、 π = 4m/n
もしあなたが、砂粒をたくさんランダムにまいて、正方形の中の砂粒と扇形の中の砂粒の数を数えることができれば、πの値がわかるということになる。自分で実際に実験することを想像すると、恐ろしく大変だが、そういう面倒な実験はコンピュータにやらせてしまえばいい！
さて、砂粒を沢山まいて、それが円の中か外かを判定する実験を、コンピュータにやらせるにはどうすればよいだろうか。
まず、点（砂粒）を大量に発生させる必要がある。さて、どうすればいいだろうか？
もうお分かりのように、「繰り返し」命令（for文）を使う。
#とりあえず、1,000つぶの砂をまくことをfor文で表現する
for (i in 1:____){
print(i) #＜砂を１粒図形の上にまく＞
}
次に、砂を１粒図形の上にまいて、それが扇形の中か外かを判定することを考える。
砂を１粒まくことは、無作為に（ランダムに）図形の上に点を打つことと同じことだ。「ランダム」というと、我々はすでに乱数を発生させる方法をしっている。図形の上に１つの点を打つということは、２つの乱数を発生させることで表現できるだろう。
例えば、0-1の範囲の乱数を２つ発生させたら、 (0.230 , 0.782) という２つの数字が得られたとする。この２つの数字が座標上の点を表していると考えると、上の図のようになる。
原点からこの点までの距離が1より小さいとき、この点は扇形の中に入っているということができる。
Rでは平方根を表すとき、 sqrt()という関数を使うので、この例では、
sqrt(0.230^2+0.782^2) < 1
となり、ランダムに発生させた（図形の上に偶然落ちた）1つの点は、扇形の中だったということになる。
実際に２つの乱数を発生させて、それを、pointというオブジェクトに代入するには、次のようにする。
point=runif(2) #この命令で、0-1の範囲の乱数を２つ発生させ、pointというベクトルに入れる
このとき、pointに入ったそれぞれの値を指定するには、
point[1] point[2]
のようにカギ括弧［］を使う。
例： []の使い方を、試してみよう
> x=runif(2)
> x
> x[1]
> x[2]
砂粒（点）が扇形の中かどうかを判断するにはどうすれば良いか？これは条件分岐で表現できる。つまり、runif(2)で発生させた２つの数値のそれぞれを2乗して足しあわせたものの平方根が1より小さいか否かを判断する。
#扇形の中かどうかを判断する条件文
distance=sqrt(point[1]^2+point[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)
}
10万回の実行結果。
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)
}
†次に、生物学科の学生なら、入試のときに必ず勉強している遺伝的浮動のシミュレーションをやってみよう。
遺伝的浮動とは、
大きさの決まった集団で起きる、配偶子のランダムサンプリングによる遺伝子頻度の変動
のことだ。例えば、集団中の２つの対立遺伝子を、青玉と赤玉で表すと、袋（配偶子の集団： 配偶子の数はものすごく多く、青と赤の比率はそれぞれ50%）から取り出される青玉・赤玉の割合（集団中の遺伝子頻度）には、下の図に示したような変動が生じる（下の例で、青玉の遺伝子頻度は、0.8, 0.7, 1.0, と変動している）。もし、袋から取り出されたのが青玉ばっかり（あるいは、赤玉ばっかり）だと、それ以降、何度繰り返しても変化は生じない（下の例だと、「集団は青玉に固定した」などという）。
それでは、赤玉と青玉が同じ数だけたくさん入った袋から、10個の玉を取り出して、その中に含まれる青玉の割合に従って袋に含まれる青玉配偶子の個数を変え、また10個取り出す。。。こういう操作を何度も繰り返すと、取り出した玉のうちの青玉と赤玉の割合の変動がわかる。
自分で何度も実験を繰り返してもいいけれど、とても面倒！ そこで、この実験をコンピュータにやらせ、取り出した玉の割合がどのように変化するかをグラフにしよう！
質問： この実験を何度も繰り返しておこなうとき、取り出した玉が赤か青かどちらか１色になってしまうのは、 何回目だと思いますか？ （袋の中の青玉と赤玉の数は同数からはじめ、 袋から取り出すのは10個）
もう皆さんは繰り返し、代入、条件分岐というプログラミングの基本技を修得しているので、今問題にしている実験を、この３つのワザを使ってどのように表現するかが、プログラミングの重要なポイントだ。
「赤玉と青玉の同じ数入っている袋から無作為に10個の玉を取り出したら、青玉が7個入っていた」 ----------------------------------------------------- 上でやった操作には、「10個取り出して青玉かどうか判断する」というところに「繰り返し」が含まれている。 「繰り返し」がはっきりと分かるように説明すると、 1. 袋から玉を一つ取り出す 2. その玉が青玉だったら数える（１つ目だったら「一つ」、２つめは「２つ」..)。赤玉だったら数えない 3. 1 と 2 の操作を10回繰り返すこの部分を大雑把でいいので、プログラミングしてみる。繰り返しはfor文で表現できる
カウンタの値は最初は0にしておく
for ( 10回繰り返す ) {
if (取り出したのが青玉だったら) {
カウンタの値を１増やす
}
}
カウンタのところは代入文で表現できる。for文も、Rで実行できる形式に直す
カウンタ = 0
for ( i in 1:10 ) {
if (取り出したのが青玉だったら) {
カウンタ = カウンタ + 1
}
}
つまり、
0-1の範囲の乱数を１つ発生させ、それが0.5よりも小さい値だったら、 青玉が取り出されたとみなす乱数を１つ発生させるには runif(1) を使って、上のプログラム「取り出したのが青玉だったら」の部分を書き換える
カウンタ = 0
for ( i in 1:10 ) {
if (runif(1) < 0.5 ) {
カウンタ = カウンタ + 1
}
}
print(カウンタ) #この命令で'カウンタ'という変数の内容を表示
これで10個取り出したうちの青玉の個数を数えるという部分は完成！
上の図で黄色く塗った部分：「袋から取り出した青玉の数を数えるシミュレーションのプログラミング」は終わった。次に、取り出した青玉の割合になるように袋の中の配偶子の割合を変更して、同様の操作を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に固定することがわかる。 青玉の割合=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(結果)
}(0.65, 0.7, 0.75, 0.85, ....)という形式にしておく必要がある。そこで、最初に空ベクトルを１つ作っておき、カウンタの値をベクトルに保存し、最後に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") #type="l" が折れ線グラフにするオプションに置き換えるだけ。
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に固定する実行してみると、うまくいってるようだが、１回ずつ、グラフが変わってしまう。グラフを重ねて描くには、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に固定する
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)、１集団あたり観察する世代数(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(【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. 実験繰返し10、1集団での観察世代数500、集団サイズ（遺伝子数）200, 対立遺伝子ａの数の初期値100 2. 実験繰返し10、1集団での観察世代数100、集団サイズ（遺伝子数）20, 対立遺伝子ａの数の初期値10いずれの場合についても、できたグラフをWinShotで切り取って、課題提出ページに貼り付けなさい（注：1番目のシミュレーションの実行には、数分かかるかもしれません）（注：あまり大きい図（画面全面に広がる図とか）にはしないでください）
http://bean.bio.chiba-u.jp/moodle から提出して下さい。
http://bean.bio.chiba-u.jp/moodle から提出して下さい（ プログラミングをもっとやってみたいという人のみ提出）。