![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
![[worried]](image/face/worried.png)
†![[smile]](image/face/smile.png)
†プログラミングの基本技の最後は条件分岐。if()を使って表現する。
i=4 #iに4を代入
if(i==3){ #iの値が3ならば
print("三") # 三 と表示する
} else { #iの値が3では無ければ
print("三以外") # 三以外と表示する
}
if命令の()の中の式を条件式といい、iの値を評価している。評価に使われるのは比較演算子というもので、
== 等しければ != 等しく無ければ >, >= 左辺が右辺より大きいなら、左辺が右辺以上ならば <, <= 左辺が右辺より小さいなら、左辺が右辺以下ならば
ということを意味している
''比較演算子''である''=='' に対して、代入の時に使った ''='' は''代入演算子''といいう。
上のif命令では、次のような条件分岐が行われている。
if(i==3){
print("三")
}
・もしiの値が３ならば、画面に三を表示。
・もしiの値が３で無いならば、何もしない
else { #iの値が3では無ければ
print("三以外") # 三以外と表示する
}
i=3 #iに3を代入
if(i==3){ print("三") } #iが3ならば三と表示上の例では、iの値が最初から決まっていたので、出てくる結果も１通りにしかならないため、プログラムとしては面白く無い。そこで、乱数を１つ発生させ、値が0.5以上なら「0.5以上」、0.5より小さいなら「0.5より小さい」と表示させるプログラムを考えてみよう。
乱数を1つ発生させるにはrunif()という関数を使う。試しに
runif(5)
と入力してみよう。0から1までの乱数が5個表示されたはずだ。乱数はシミュレーションではよく使うので、もう少し慣れ親しんでおこう。 例えば、runif()を使って乱数を10個発生させ、ヒストグラムを書いてみよう。同じことを、100個、1000個、10000個、100000個についてもやってみよう。
hist(runif(10))
バラツキがだんだんと無くなってきたのがわかるだろうか。
#穴埋め式練習問題
if (runif(1) >= _ ) { #条件式0.5以上ならば
print("0.5以上")
} _ {
print("0.5より小さい")
}
では、上の繰り返しと条件分岐を組み合わせて使ってみよう
練習問題：for命令を使って1から10までの数字を表示させなさい、5だけは"five"と英単語で表示させなさい
for(i in _:_){ # 10回繰り返し。iは1から10まで変化
if(i == _) { #iが5の場合
print("five") #fiveと表示
} _ { #それ以外の場合
print(_) #iの値を表示
}
} †「関数の定義」なんていうと難しそうだが、ようするに、先ほどまでに作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにしようというものだ。こういうときに関数を定義する関数、function()を使う。 例えば、円の面積を計算する関数を作るなら、
menseki=function(r){r*r*pai}
とすれば良い。これで自分独自の関数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)
goukei=0 #goukeiという変数に初期値0を代入
for(i in 1:10) {
goukei=_____+i #goukeiに前回までのgoukeiの値にiの値を足したもの代入
}
print(goukei)
選択肢： a) i b) goukei c)
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) <=
| 演算子 | 意味 | 使用例 |
| + | 足し算 | > 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 |
> 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を用いたシミュレーションに挑戦してみよう。
シミュレーションという言葉は「株式シミュレーション」とか「シミュレーションゲーム」とかで聞いたことがあるだろう。日本語にすると、「模擬実験」で、自分では実行が難しい繰り返し実験を、コンピュータに計算させることで行うものだ。コンピュータは単純な繰り返し計算は得意なので、砂粒をばらまいて、ある図形の中に入った数を数えるなんていう作業を、10万回でも、100万回でも繰り返すことができる。
左に１辺の長さが1の正方形と、それに内接する扇形がある。このとき、
正方形の面積は 1 扇形の面積は π/4
この図形の上に、砂粒をばらまいたところ、
正方形の中には33粒、そのうち扇形の中には25粒入っていた
砂粒が図形の上にランダムに散らばるとすると、砂粒をすごく沢山まいて数を数えれば、
正方形の中の砂粒の数(n) ： 扇形の中砂粒の数(m) ＝ 正方形の面積 ：扇形の面積 ＝ １ ： π/4
になると考えられる。つまり
n : m = 1 : π/4 だから、 π = 4m/n
もしあなたが、砂粒をたくさんランダムにまいて、正方形の中の砂粒と扇形の中の砂粒の数を数えることができれば、πの値がわかるということになる。
さて、砂粒を沢山まいて、それが円の中か外かを判定する実験を、コンピュータにやらせるにはどうすればよいだろうか。
ここで乱数を用いる。
今、0-1の範囲の乱数を２つ発生させたら、 (0.230 , 0.782) という２つの数字が得られたとする。この２つの数字が座標上の点を表していると考えると、上の図のようになる。
原点からこの点までの距離が1より小さいとき、この点は扇形の中に入っているということができる。この例では、
0.230^2+0.782^2 < 1
なので、ランダムに発生させた（図形の上に偶然落ちた）1つの点は、扇形の中だったということになる。
ではこのような点を大量に発生させるにはどうすればいいだろうか？
もうお分かりのように、「繰り返し」命令（for文）を使う。
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)
}
赤玉と青玉が同じ数だけたくさん入った筒から、20個の玉を取り出して、その中に含まれる赤玉の割合に従って筒に含まれる赤玉の個数を変え、また20個取り出す。。。こういう操作を何度も繰り返すというモデルで、何度も実験すると、袋の中の玉の割合がどのように変化するかをシミュレーションで示してみよう。
「赤玉と青玉の同じ数入っている筒から無作為に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文をどうプログラミングすれば良いかが、今回のシミュレーションの最重要ポイントだろう。この部分は乱数を使って表現できる。
つまり、0-1の範囲の乱数を１つ発生させ、それが0.5よりも小さい値だったら、赤玉が取り出されたとみなせばいい。
乱数を１つ発生させるには runif(1) を使う。
カウンタ = 0
for ( i in 1:20 ) {
if (runif(1) < 0.5 ) {
カウンタ = カウンタ + 1
}
}
print(カウンタ) #この命令で'カウンタ'という変数の内容を表示
これで20個取り出したうちの赤玉の個数を数えるという部分は完成です
上の図で、最初の矢印の部分の操作を行い、筒から取り出した赤玉の数を数えるシミュレーションのプログラミングは終わりました。次に、取り出した赤玉の割合になるように筒の中の赤玉の数を変更して、同様の操作を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(カウンタ) #この命令で'カウンタ'という変数の内容を表示
}
あとは結果をプリントする文を入れれば完成です。何回目に固定したかが分かるように、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(結果)
}
一昨年、同様の授業を行ったところ、受講生から次のような質問がありました。
「こういうシミュレーションのプログラミングって、これから生物学を学ぶのに、どのくらい必要なんですか？」
これまでプログラミングなどしたことの無い人にとって、とても素朴な疑問だろうと思います。そこで、次のように答えておきました。
画一的な答えは難しいです。シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。 ただ、私は、プログラミングは、生物学における発展的な情報処理をする上での基本技だと考えています。 千葉大の研究室に限って言えば、プログラミング技術が無いと卒業研究が終わらないというところは 現時点では（たぶん）無いでしょう。但し、プログラミングができると、卒業研究や修士の研究でより 進んだことができたり、理解が深まったり、研究が楽になるようなところはいくつもあります。 将来、数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、 敢えて、不可欠な技術だと言っておきます。
./11
*第11回授業・基本課題 **氏名： **課題への回答 -今日の授業の進み方は？（はやい、丁度いい、おそい） --回答： -今日の授業の難しさはどう感じましたか（簡単 丁度いい 難しい）： --回答： -難しいと答えた人は、特にどの点が難しかったですか？： --回答： -今日の授業は（分かった 半分ぐらいは分かった 分からなかった）： --回答： -分からないと答えた人は、特にどの点が分からなかったですか？： --回答： -今日の講義で分からなかった用語があったら挙げてください： --回答： -授業に関する要望・質問があったらなんでもどうぞ： --回答：
下のプログラムは遺伝的浮動をシミュレートする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"))
}
num_repeats, num_generations, size_population, num_a_allele
1. 実験の繰り返し数10回、1集団での観察世代数500世代、集団サイズ（遺伝子数）200, 初期状態での対立遺伝子ａの数100 2. 実験の繰り返し数10回、1集団での観察世代数100世代、集団サイズ（遺伝子数）20, 初期状態での対立遺伝子ａの数10 いずれの場合についても、できたグラフをWinShotで切り取って、課題提出ページに貼り付けなさい（注：1番目のシミュ レーションの実行には、数分かかるかもしれません）（注：あまり大きい図（画面全面に広がる図とか）にはしないでく ださい）