- 追加された行はこの色です。
- 削除された行はこの色です。
*！編集中! Rを使ったプログラミング演習2： シミュレーション [#uf5bc6d0]
*Rを使ったプログラミング演習 1 ：プログラミングの基礎[#zc5b0e21]
#contents
はじめに~
前々回から使いはじめたRというソフトウェア。コマンドを入力するのはちょと難しいけど、エクセルみたいにマウスを使って面倒な操作をしなくても、一瞬にして統計解析をやってくれるし、グラフも簡単に描けてしまう優れもの。これなら、レポートや卒論のために、何度も解析をやりなおして、グラフを描き直すという作業も面倒ではないだろう。
初めて使ったソフトなのに、受講生の皆さんは、すでにRの操作には慣れてきているし（全員、課題ができていた）、興味も感じているよう。そこで、今回はRを使ったプログラミングに挑戦する。
Rでプログラミングの授業を始めたころ、受講生から次のような質問があった。
「シミュレーションのプログラミングって、これから生物学を学ぶのに、どのくらい必要なんですか？」
**第11回授業の獲得目標：&worried; [#z281e914]
-''1. プログラミングとは何かを理解する''
-''2. プログラミングにおける%%%繰り返し%%%と%%%代入%%%と%%%条件分岐%%%を修得する''
-''3. Rで基本的なプログラムを作成し、実行できるようになる''
これまでプログラミングなどしたことの無い人にとって、とても素朴な疑問だろうと思う。私からの答えとしては、
シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。
プログラミングは、生物学における発展的な情報処理として、重要なものだと思います。
将来、数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、
きっと必要になるでしょう。
**プログラミングって何？⌣ [#if4bd58f]
プログラミングっていったい何をすることだろうか？IT用語辞典で調べてみると、、、
**第11回授業の獲得目標：&worried; [#j12b03a4]
-''1. プログラミングの基本：繰り返し・代入・条件分岐に習熟する''
-''2. ユーザー定義関数を理解する''
-''3. シミュレーションの考え方を習得する''
-''4. シミュレーションのプログラミングを体験する''
-''5. 遺伝的浮動のシミュレーションに挑戦する！''
載っていない
唯一近いところでは、「プログラミング言語」というのがある。IT用語辞典の説明文からプログラミングに関わるところを抜き出して、「プログラミング」の説明文を作ってみると、
**繰り返し・代入・条件分岐のおさらい ⌣ [#ma15d60f]
次のそれぞれのプログラムの＿部分に何を入れるべきかを選択肢から選びなさい。また、Rを使って実行してみよう。
***繰り返し [#uaa3b826]
-演習： 1から10までの数字を表示させるプログラムを作りなさい
for(i in 1:__){ #{と}の間を10回繰り返す。i の値は毎回1ずつ増える
print(i)
}
選択肢： a) 1 b) 5 c) 10
***代入 [#wb6c9fb4]
-演習： 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
-''プログラミング'': 人間に理解できるように英語などを元に作られたプログラミング言語を用いて、コンピュータに実行させる命令を作ること。できた命令はソースコード、プログラム、スクリプトなどの名前で呼ばれる。命令文をコンピュータに実行させるためには、コンピュータが理解可能な機械語の羅列(オブジェクトコード)に翻訳する必要がある。プログラミング言語によっては、コンパイルという翻訳作業を別にしなければならないコンパイラ型言語と、リアルタイムに翻訳作業をしてコンピュータに命令を渡すインタプリタ型言語がある。
***条件分岐 [#p7f27ed4]
-演習： 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) <=
つまり、プログラミングとは、簡単に言うと、「%%%言葉を使ってコンピュータが理解できる命令を作ること%%%」。
今日使うRは統計解析ソフトでもあるけれど、インタプリタ型のプログラミング言語でもある。
***プログラミングTips [#ka02fd43]
・１行に１つの命令を書いたら、必ず改行
・{ } のそれぞれの後で改行すると見やすい（このページの例を参照）
**おさらい： プログラミングは手順書き！ 子供にお買い物を言葉で教えるようなもの [#g3fd6e2c]
前の授業で、&size(16){「''プログラミングは数学というよりも、コンピュータにことばで手順を伝えることだ''」};という説明をした。~
例えば、おつかいに行く子供に、次のような買い物の手順を説明してみよう。
お財布の中に500円入っているから、これを持って買い物に行ってね。
この通りをまっすぐに進むとお店が３軒あるから、順番に買い物をするんだよ。
最初のお店で、100円のゴーヤを１本もらって、100円を渡すんだよ。
次のお店で100円のお豆腐を一丁もらって、100円を渡すんだよ。
次のお店で100g 100円の豚肉を、200gもらって、200円を渡すんだよ。
残ったお金はいくらだった？あとでお菓子をかってもいいわよ。
このお買い物の手順を、先週学んだ「繰り返し・代入・条件分岐」 を使って書いてみると、例えば、下のようなプログラムになる。
お財布の中 = 500
for (i in 1:3) { #３回買い物を繰り返す
if (i != 3) { #3件目以外は100円ずつ使う
お財布の中 = お財布の中 - 100
} else {
お財布の中 = お財布の中 - 200
}
}
print(お財布の中)
プログラミングが「手順」を伝えるものであることが、分かって貰えただろうか？
**コンピュータにどんな命令ができるのか？ or どんな命令がしたいか？ [#g9c1a760]
「Wordを立ち上げて文章を作成し、印刷する」とか、「Excelで家計簿を管理する」とかは、市販のソフトウェアを使った作業。いずれも、非常に高度で複雑な命令が、使用者からコンピュータに向かって発せられている。例えば、~
・マウスを動かしてマウスカーソルを動かし、
・メニューをクリックして印刷を選択して、
・ファイルを印刷する
なんていうのは、かなり複雑な「命令」の集まり。でも、多くの場合、人がコンピュータで行っている作業の多くは、「命令」を意識しなくても、「カーソルを動かす」とか「クリックする」など、利用者が理解しやすく・使いやすいように設計されている。。
**ユーザー定義関数⌣ [#tca11c6a]
「ユーザー定義関数」なんていうと難しそうだが、ようするに、先ほどまでに作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにしようというものだ。こういうときに関数を定義する関数、''function()''を使う。
例えば、円の面積を計算する関数を作るなら、
menseki=function(r){r*r*pi}
とすれば良い。これで自分独自の関数menseki()ができた。
では、実行してみよう。
menseki(10) #好きな数字を入れて円の面積を計算
でも、本当に市販のソフトウェアを使うだけでいいのだろうか？良くないという点が、すぐに２つ思く：
1. 自分の目的にあったソフトウェアがいつも存在するというわけではない
2. ソフトが存在する場合、目的の数だけ、ソフトの使い方を覚えなければならない
''演習 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を表示
} #関数定義の終了
そこで、この授業で次の３つを目的としてプログラミングを勉強する。
+自分で作った独自の命令をコンピュータに与えて、目的を達成する
--プログラミングを行って命令を実行させることで、初めて、''コンピュータに作業をさせている''ことが実感できる
--ワープロや表計算とかは、どちらかというと、道具を使って作業をさせられている感が強い
+プログラミングを行うことで、コンピュータができることの理解が深まる
--あなたのコンピュータは実は、すごい性能をもっている。それを理解しておかないと使い倒せない
+プログラミングの過程で、コンピュータが行う処理の意味を論理立てて考えることを学べる
--コンピュータはすごく石頭、物事を本当に論理立てて説明してやらないと、命令を聞いてくれない
&br;
実を言うと、生物学の研究者でも、プログラミングなどを経験したことの無い人は大勢いる。今や、DNAシーケンスの決定や整列、系統樹作成にだって、専用のソフトウェアがある。。。しかし、そういうソフトウェアも、多くの場合は研究者自身が作成したもの（学生が作ったものもたくさんある）。この授業では、頭の柔らかい１年生のうちにプログラミングについての動機付けを行うことで、このクラスの中から将来は新しいソフトウェアを開発できる人が出ることも期待している。。
**授業で使ったRの基本関数 [#h60d69ca]
//#include(授業/20/情報処理/R関数一覧,notitle)
→[[授業/H20/情報処理/R関数一覧]] これまでに学んだRのいろんな関数を一覧できる。
**初めて(?)のプログラミング："Hello World!" [#m2b829c6]
「Hello World!」だなんて、なんか変なタイトル。誰が使い始めたのかは知らないが、多くのプログラミング言語の教科書で、最初に出てくるのがこのプログラム（プログラムを勉強したことのある人なら、誰でも知っている）。では、早速やってみよう。Rを立ち上げて、次のように入力。
**Rを用いたシミュレーション⌣ [#pb5e91d7]
print("Hello World!")
これまでに学んだ繰り返し・代入・条件分岐の３つの命令それぞれは単純だが、組み合わせれば、かなり複雑なことでも表現できる。そこで、今回は、これらの命令を組み合わせて、''Rを用いたシミュレーション''に挑戦してみよう。~
シミュレーションという言葉は「株式シミュレーション」とか「シミュレーションゲーム」とかで聞いたことがあるだろう。日本語にすると、「模擬実験」で、自分では実行が難しい繰り返し実験を、コンピュータに計算させることで行うものだ。コンピュータは単純な繰り返し計算は得意なので、砂粒をばらまいて、ある図形の中に入った数を数えるなんていう作業を、10万回でも、100万回でも繰り返すことができる。下の例では、そういうシミュレーションの方法を用いて、円周率（π）の値を求める。このように、シミュレーションをうまく使えば、たとえ数学的に解くのが難しい問題についても、近似解を求めることができる。また、ある現象について作成したモデルを用い、パラメーターの変化がどのような状態の変化をもたらすかを予測することができるため、進化生物学や生態学でもよく使われる。~
~
今回は、乱数を用いたシミュレーションを行う。乱数とは、ランダム（無作為）に得られる数値のことで、次にどんな数字が現れるのか予測できないような数字の集合のことを言う。Ｒでは、
runif()
という関数を使って、好きな数だけ乱数を得ることができる。~
例えば、
runif(5)
と入力すれば、乱数が５個得られる。せっかくなので、得られた乱数をグラフにすることで、視覚的に乱数の性質を理解してみよう。
plot(runif(1000)) #枠内に、1000個の点がランダムにプロットされる
hist(runif(1000) #発生させた1000個の乱数のヒストグラム； 1000個ぐらいだと、けっこうばらつく
#数を増やすと、バラツキが小さくなるのが分かるだろうか？（これを大数の法則という）
うまくご挨拶できただろうか？
***円周率(π）を「繰り返し」・「乱数」・「条件分岐」で求める [#u365271e]
&ref(授業/H20/情報処理/11/#10_1.jpg,60%);~
左に１辺の長さが1の正方形と、それに内接する扇形がある。このとき、
正方形の面積は 1 扇形の面積は π/4
この図形の上に、砂粒をばらまいてみる。
砂粒が図形の上にランダムに散らばるとすると、砂粒をすごく沢山まいて数を数えれば、
正方形の中の砂粒の数(n) ： 扇形の中砂粒の数(m) ＝
正方形の面積 ：扇形の面積 ＝
１ ： π/4
になると考えられる。つまり
n : m = 1 : π/4 だから、 π = 4m/n
もしあなたが、砂粒をたくさんランダムにまいて、正方形の中の砂粒と扇形の中の砂粒の数を数えることができれば、πの値がわかるということになる。自分で実際に実験することを想像すると、恐ろしく大変だが、そういう面倒な実験はコンピュータにやらせてしまえばいい！
ここでやったのは、print()という関数＜注：''命令''とほとんど同じ意味＞を使って"Hello World!"と画面に表示させただけだが、これも
画面にHello Worldと表示させる
ということを目的とした立派なプログラムだ。
~「えっ？それなら、Rの最初の授業でやったオブジェクトの内容を画面に出すのもプログラミング？」と聞かれるかもしれませんが、その通り。次の囲みの中を全てコピーして、Rのコンソールにペーストしてみよう。
x=c(1,2,3,4,5,6,7,8,9,10)
sum(x)
上でやったプログラミングは、
画面に1から10までの整数の合計を表示させる
というもの。先週やったRを使った操作も、何らかの処理の結果を画面に表示させるプログラミングだったわけだ。なんとなく、プログラミングが身近になっただろうか？
***シミュレーションの考え方 [#y94eb6a1]
さて、砂粒を沢山まいて、それが円の中か外かを判定する実験を、コンピュータにやらせるにはどうすればよいだろうか。
まず、点（砂粒）を大量に発生させる必要がある。さて、どうすればいいだろうか？~
もうお分かりのように、「繰り返し」命令（for文）を使う。
#とりあえず、1,000つぶの砂をまくことをfor文で表現する
for (i in 1:____){
print(i) #＜砂を１粒図形の上にまく＞
}
**複数行の命令文 [#r702dbae]
上でやった２つのプログラミングは、見比べてみると違いがある。"Hello World!"の方は１行だけの命令だが、sumの方は2行になっている。プログラムは通常、いろんな処理を組み合わせて作るので、複数行になることが多い。~
では３行からなるプログラムを作ってみよう。
q=c(10,15,20,25,30,35,40,45,50,55,60,65,70,80,100)
y=c(2,3,4,5,6,7,8,9,10,11,10,11,9,10,12)
plot(log="xy",log10(q^2), log10(y))
このプログラムは、
２つの数値ベクトルについて計算を行って、両対数グラフを描く
というものだ。使ったデータは、まるでどこかの大学の生物学科１年生がコドラート調査でとったようなもので、面積が広くなると種数がどう増えるかを示したものだ。Rを使えば対数グラフ用紙を使って手作業で面倒なプロットをする必要もなくなる。
次に、砂を１粒図形の上にまいて、それが扇形の中か外かを判定することを考える。~
**Rで使う命令 [#q91e1e2f]
授業では一つ一つ解説するが、復習課題をやるときとか、自分で新たな解析に挑戦するときは、どういう命令があるかが書かれた説明書が必要になる。前回も紹介したが、
-http://cse.naro.affrc.go.jp/takezawa/r-tips.pdf ~
には船尾さんという方が書かれた、日本語のRの解説がある。知りたい項目はAdobe Readerの検索機能でサーチしてみよう。
また、Rを使っている最中は
> ?plot
> example(plot)
という風に、?の後に関数名を打てば、説明や使用例が英語で表示される。説明文の一番最後にはExample（使用例）が載っているので、けっこう参考になる。
砂を１粒まくことは、無作為に（ランダムに）図形の上に点を打つことと同じことだ。「ランダム」というと、我々はすでに乱数を発生させる方法をしっている。図形の上に１つの点を打つということは、２つの乱数を発生させることで表現できるだろう。
&ref(授業/H20/情報処理/11/#10_2.jpg,60%);~
例えば、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]
**プログラミングの基礎: 繰り返しと代入と条件分岐⌣ [#p8b86305]
上の例で、紹介したプログラムは、個々の命令を順番に並べて、一度にペーストしただけだった。こんなことなら、Rに命令を1行ずつ与えるのと、ぜんぜん変わらない。そこで、もっと%%%プログラムらしい命令%%%をRに与えてみよう。それは、''繰り返し''と''代入''と''条件分岐''というもので、プログラミングの基本中の基本。
砂粒（点）が扇形の中かどうかを判断するにはどうすれば良いか？これは条件分岐で表現できる。つまり、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増やす
}
***繰り返し⌣ [#b97e46d7]
人が不得意でコンピュータが得意なことの一つは、単純な繰り返し作業を際限なく（文句も言わずに）行うことだろう。つまり、プログラミングの一つの目的は、面倒な繰り返し計算をコンピュータにやらせることだと言っても言い過ぎではない。では、%%%何回繰り返すか%%%という命令をどうやってコンピュータに与えるかというと、 ''for'' という命令を使う。
以上のように、これまでに学習した、繰り返し、条件分岐、代入で、図形の上に砂粒を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)
-''繰り返しには for() を使う''
for(i in 1:10) { print(i) }
この命令は、
{ } で囲まれた部分の命令を10回繰り返しなさい
というものだ。
{ } の中は何をする命令か分かるだろうか？
さらにユーザー定義関数：''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)
では、上の囲みの中をコピーして、Rのコンソールにペーストしよう。~
1から10までの整数が表示された。ここでもう一度命令文をよく見てみると、
for と (i in 1:10) と { print(i) }
という３つのパートからなっている。
for は これから繰り返し命令が始まるよということを示す
(i in 1:10) は 繰り返しの回数が10回であることを示す
また、iは1回の繰り返しのたびに1ずつ、10まで増加する
もし、50回繰り返したかったら (i in 1:50) と書く
iという文字の代わりに、jでも、kaisuでも、構わない
{ print(i) } は 繰り返しの度に{ }の中を実行せよという命令で、
print(i)で i という変数の値を画面に表示させなさいという意味
{ } のところは複数行に分かれていてもよく
for(i in 1:10) {
print(i)
}
> 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
と書いても結果は同じ。
***せっかくだから、点をうつところを図で表現してみよう [#lec8b399]
''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万回の実行結果。
なお、1:10 は
1 2 3 4 5 6 7 8 9 10
という数字の集まりを示す。次の囲みの中を１行ずつRに実行させてみよう。
1:10
x=c(1:10)
x
**参考： 円周率の計算を乱数を使わずに行うプログラム [#deed5be0]
-これまでの説明が簡単すぎて時間の余った人は、次のような考え方でπを求めてみよう。
-上のプログラムとよく似ているが、乱数を使わず、1x1の正方形をマス目に分割して、扇形の中に入っている交点の数を数えるもの。穴埋め問題にしておくので、プログラムを完成させてほしい。得られる値は、100のとき3.1、1000のとき3.137524、10000のとき3.14119。100000のときは計算に時間がかかりすぎて途中でストップ。
練習問題： 2から18までの偶数を全て表示させるプログラムをfor命令を使って作りなさい
--プログラムのことをスクリプトとも言いう。Rのファイルメニューから「新しいスクリプト」を選ぶと、「スクリプトエディター」というエディターが開く。スクリプトエディターでスクリプトを書いてから、マウスで選択して右クリックすると、選択中のスクリプト（命令とかコードとも呼ばれる）が実行できる。
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)
解答例：（穴埋め形式）
for(i in _:_){ # _ のところには数字が入る
print(i*_)
}
**遺伝的浮動のシミュレーションに挑戦する！⌣ [#ne685455]
次に、生物学科の学生なら、入試のときに必ず勉強している''遺伝的浮動''のシミュレーションをやってみよう。~
遺伝的浮動とは、
大きさの決まった集団で起きる、配偶子のランダムサンプリングによる遺伝子頻度の変動
のことだ。例えば、集団中の２つの対立遺伝子を、青玉と赤玉で表すと、袋（配偶子の集団： 配偶子の数はものすごく多く、青と赤の比率はそれぞれ50%）から取り出される青玉・赤玉の割合（集団中の遺伝子頻度）には、下の図に示したような変動が生じる（下の例で、青玉の遺伝子頻度は、0.8, 0.7, 1.0, と変動している）。もし、袋から取り出されたのが青玉ばっかり（あるいは、赤玉ばっかり）だと、それ以降、何度繰り返しても変化は生じない（下の例だと、「集団は青玉に固定した」などという）。~
&ref(授業/H21/情報処理/11/#11_1.gif,50%);
***代入⌣ [#z81c679c]
上の for(i in 1:10) { print(i) } というプログラムでは、iの値は繰り返しの度に違う値になる。
1回目 iの値は1
2回目 iの値は2
3回目 iの値は3
.... ........
10回目 iの値は10
では、一つ前の回のiの値を使いたい場合はどうすれば良いだろうか？例えば、1から10までの合計をfor命令を使って計算することを考えると、
//***集団遺伝学的問題：「固定までにようする時間」の求め方 [#ue34022b]
//上の例では、10個の遺伝子からなる成体の集団が、青玉：赤玉＝0.5:0.5という遺伝子頻度から始めて、青玉のみに固定するまでに数世代かかっている。それでは、ここで上と同じ状況を考えて、問題を１つ。
// 青玉5, 赤玉5からなる成体の集団が、青玉あるいは赤玉のいずれかに
// 固定するまで要するの世代数は、平均すると何世代かを答えなさい
//
//
***遺伝的浮動のシミュレーションをプログラミングして、赤玉と青玉の変動の度合いをグラフ化する [#gb07b277]
それでは、赤玉と青玉が同じ数だけたくさん入った袋から、10個の玉を取り出して、その中に含まれる青玉の割合に従って袋に含まれる青玉配偶子の個数を変え、また10個取り出す。。。こういう操作を何度も繰り返すと、取り出した玉のうちの青玉と赤玉の割合の変動がわかる。~
自分で何度も実験を繰り返してもいいけれど、とても面倒！
そこで、この実験をコンピュータにやらせ、取り出した玉の割合がどのように変化するかをグラフにしよう！
-でも、シミュレーションを始める前に、
質問： この実験を何度も繰り返しておこなうとき、取り出した玉が赤か青かどちらか１色になってしまうのは、
何回目だと思いますか？
（袋の中の青玉と赤玉の数は同数からはじめ、 袋から取り出すのは10個）
現在のiの値に、１つ前のiの値の合計を加えて表示させたい
1回目 iの値は1 １つ前までの合計は0 1 を表示
2回目 iの値は2 １つ前までの合計は1 3 を表示
3回目 iの値は3 １つ前までの合計は3 6 を表示
.... ........
10回目 iの値は10 １つ前までの合計は45 55 を表示
***シミュレーションの考え方 [#g214160e]
もう皆さんは''繰り返し、代入、条件分岐''というプログラミングの基本技を修得しているので、今問題にしている実験を、この３つのワザを使ってどのように表現するかが、プログラミングの重要なポイントだ。
こういう場合に、''代入''という命令を使う。代入は、前回すでにオブジェクトへの数値やベクトルの代入で学んだように、
=
を使う。我々が通常行う計算では、'='は「3+5=8」というように、「左辺の計算結果が右辺に等しくなる」というる意味で用いられるが、プログラミングにおいては、
左辺の変数に右辺の計算結果を代入する
という意味で用いられる。上の合計の計算プログラムは次のようになる。
goukei=0 #goukeiという変数に初期値0を代入
for(i in 1:10) { #以下を10回繰り返す
goukei=goukei+i #goukeiに前回までのgoukeiの値にiの値を足したもの代入
print(goukei) #goukeiの内容を画面に表示
}
では、上の囲みの中をRに実行させてみよう。
+'' まず袋から玉を10個、ランダムに取り出すところを考える''~
~&ref(授業/H21/情報処理/11/#11_2.gif,50%);
「赤玉と青玉の同じ数入っている袋から無作為に10個の玉を取り出したら、青玉が7個入っていた」
-----------------------------------------------------
上でやった操作には、「10個取り出して青玉かどうか判断する」というところに「繰り返し」が含まれている。
「繰り返し」がはっきりと分かるように説明すると、
1. 袋から玉を一つ取り出す
2. その玉が青玉だったら数える（１つ目だったら「一つ」、２つめは「２つ」..)。赤玉だったら数えない
3. 1 と 2 の操作を10回繰り返す
この部分を大雑把でいいので、プログラミングしてみる。繰り返しはfor文で表現できる
カウンタの値は最初は0にしておく
for ( 10回繰り返す ) {
if (取り出したのが青玉だったら) {
カウンタの値を１増やす
}
上の命令の中には２箇所、代入が使われている。最初の行では、goukeiという変数に0を代入するという単純なもので、goukeiの初期値を決めている。~
初期値というのも耳慣れない言葉かもしれない。プログラムに登場する変数は、最初の値が何かをはっきりコンピュータに教えてやらないと、メモリに残っている数値がそのまま使われてしまう。試しに、上のプログラムから１行目を省いたプログラムを何度か実行させてみよう。初期値が入っていないので、合計値がどんどん加算されてゆくのが分かる。
次の、''goukei=goukei+i''という命令文では、%%%右辺の「goukei+i」を先に計算%%%して左辺の「goukei」に代入している。~
''「右辺と左辺に同じ変数が出てくる」''ところが分かりにくいかもしれないが、理解して欲しい。言葉で説明すると、
今の goukei の値に 1 を加えたものを、goukei に代入する（これまでのgoukeiは上書きされる）
という意味だ。
練習問題： 1から10までの数字を１つずつ順々にかけ合わせたた結果を全て表示させるプログラムを作りなさい
（1x1=1, 1x2=2, 2x3=6, 6x4=24, 24x5=...... と計算してゆくということ）
以下、穴埋め形式でヒントを書いておく。_の部分（１文字とは限らない）を補ってプログラムを完成しなさい。
seki=_ #sekiという変数に初期値1を代入（かけ算だから）
for(i in 1:10) { #以下を10回繰り返す
_=_*_ #sekiに前回までのsekiの値にiの値を掛けたもの代入
print(seki) #sekiの内容を画面に表示
}
うまくできただろうか？
縦に長くプリントされるのは格好悪いので、%%%毎回の計算結果をオブジェクトにベクトルとして代入してから表示させてみよう%%%。
//ベクトルに要素を追加するには''append''という命令を使う。
seki=1 #sekiという変数に初期値1を代入
kekka=c() #kekkaというオブジェクトに空ベクトルを初期値として代入
for(i in 1:10) { #以下を10回繰り返す
seki=seki*i #sekiに前回までのsekiの値にiの値を掛けたもの代入
kekka=c(kekka,seki) #kekkaというベクトルにsekiの値を要素として追加
}
print(kekka) #kekkaの内容を画面に表示
ここでは、計算結果を順々に保存しておく入れ物を、kekkaという名前で作成し、初期値に空ベクトル（データは無いが、ベクトル）を作っている。ベクトルの作成は、前回の授業でやったように
c()
という関数を使う。
また、ベクトルに要素を追加するにも、同じ関数 c() を使っている。
では、上の囲みの中をRに実行させてみよう。
-２つの代入文の比較~
上のプログラムの''kekka=c(kekka,seki)''と１つ前の足し算のプログラムにあった''goukei=goukei+i''という２つの式はよく似た働きをしていることを理解して欲しい。
kekka=c(kekka,seki) kekkaにsekiを要素として追加したものをkekkaに代入
goukei=goukei+i goukeiにiを加えたものをgoukeiに代入
いずれの命令でも、現在の変数の値が''上書きされて変更されている''ことに注目。
-せっかくRを使ってプログラミングをしているので、結果をグラフに表示させてみよう。上のprint(kekka)の部分を、plot(kekka)（あるいは、plot(log(kekka)))に変えるだけで、グラフとして出力することができる。値も表示させたい場合は、両方の命令を書いておけば良い。
練習問題（穴埋め式）：
seki=1 #sekiという変数に初期値0を代入
kekka=c() #kekkaというオブジェクトに空ベクトルを初期値として代入
for(i in 1:10) { #以下を10回繰り返す
seki=seki*i #sekiに前回までのsekiの値にiの値を掛けたもの代入
kekka=c(kekka,seki) #kekkaというベクトルにsekiの値を要素として追加
}
print(kekka) #kekkaの内容を画面に表示
_(log10(kekka)) #kekkaの対数をグラフにプロット
&size(16){ずいぶんとプログラムらしくなってきた};
-for文は２重にして使うこともできる。繰り返しのループを二重にすることで、２つの変数の総当たり計算が可能になる。
#＊＊次のプログラムの実行結果を予測してみよう。＊＊
for(i in 1:9) { #iの値を1から9まで1ずつ変化させる
for(j in 1:9){ #jの値を1から9まで1ずつ変化させる
print(i*j) #iとjを掛けたものを表示する
}
カウンタのところは代入文で表現できる。for文も、Rで実行できる形式に直す
カウンタ = 0
for ( i in 1:10 ) {
if (取り出したのが青玉だったら) {
カウンタ = カウンタ + 1
}
}
+''次に、「もし、取り出したのが青玉だったら」という部分を、if文を使ってプログラミングする''~
この点が、シミュレーションの最重要ポイントだ。~
ここでは、前回の予習課題で使った''乱数''である''runif()''という関数を使う。
~&ref(授業/H21/情報処理/11/#11_3.gif,50%);~
つまり、
0-1の範囲の乱数を１つ発生させ、それが0.5よりも小さい値だったら、
青玉が取り出されたとみなす
乱数を１つ発生させるには 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, ....)
という形式にしておく必要がある。そこで、最初に空ベクトルを１つ作っておき、カウンタの値をベクトルに保存し、最後に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に固定する
予想できたら、実行しよう。
実行してみると、うまくいってるようだが、１回ずつ、グラフが変わってしまう。グラフを重ねて描くには、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(結果,青玉の頻度) #青玉の頻度を結果というベクトルに追加
-予想通りでしたか？小学生のころ呪文のように何度も唱えものが表示されたはずだ。上のままでは縦に長く表示されて格好が悪いので、matrixという関数をつかって、ベクトルを行列に変換してで出力してみよう。
''matrix()''という関数はベクトルを行列に変換する働きがある。例えば、
x=c(1, 2, 3, 4, 5, 6, 7, 8, 9)
というベクトルを、3 x 3の行列に変換したいときは、
matrix(x, nrow=3, ncol=3) #nrow= は行の数、 ncol= は列の数を指定
matrix(x, 3, 3) #nrow=, ncol=を省略して、こういう書き方もできる
#演習：九九の表をmatrixという関数を使って表示さなさい。_を埋めなさい
result=c() #resultに空ベクトルを初期値として代入
for(i in 1:9) { #iの値を1から9まで1ずつ変化させる
for(j in 1:9){ #jの値を1から9まで1ずつ変化させる
result=c(_, _*_) #iとjを掛けたものをresultの要素として追加
}
par(new=T) #グラフを一度クリアする
plot(結果, type="l", ____=c(0,1)) #★★ylimを使ってy軸の目盛りを0から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に固定する
}
}
print(result)
print(matrix(result, nrow=i, ncol=j)) #iとjはそれぞれ9で終わっているので、9x9のマスになる
***発展：様々な変数を指定できる関数を作ってしまう [#a815b0ae]
上のシミュレーションは集団の遺伝子の数が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に固定する
}
***条件分岐⌣ [#c35ba266]
プログラミングの基本技の最後は条件分岐。if()を使って表現する。
i=4 #iに4を代入
if(i==3){ #iの値が3ならば
print("三") # 三 と表示する
} else { #iの値が3では無ければ
print("三以外") # 三以外と表示する
}
***発展・演習： [#l7cc8b03]
下のプログラムは遺伝的浮動をシミュレートする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"))
if命令の()の中の式を''条件式''といい、iの値を評価している。評価に使われるのは''比較演算子''というもので、
== 等しければ
!= 等しく無ければ
>, >= 左辺が右辺より大きいなら、左辺が右辺以上ならば
<, <= 左辺が右辺より小さいなら、左辺が右辺以下ならば
ということを意味している
--r-tips.pdfで比較演算子を検索してみよう！
''比較演算子''である''=='' に対して、代入の時に使った ''='' は''代入演算子''といいう。
上のif命令では、次のような条件分岐が行われている。
if(i==3){
print("三") # ・もしiの値が３ならば、画面に三を表示。
} else { #もしiの値が3では無ければ
print("三以外") # 三以外と表示する
}
-問1: プログラムが正常に動くよう、【1】〜【4】を補いなさい。
--なお、プログラム中、それぞれの変数は次のような意味を持っています
---num_repeats : シミュレーションを行う集団数（１つの集団についてある世代数だけ観察するという実験を、何回（何個の集団について）おこなうか）、すなわち観察集団数。何本の折れ線グラフを表示させるかはこの値で決まる。
---num_generations : 1つの集団について何世代続けて観察するか。x軸の長さになります。
---size_population : 集団のサイズ、つまり、１つの集団に含まれる遺伝子の総数
---num_a_allele : aアリルの数（つまり、この値をsize_populationで割ったものが、0世代目でのaアリルの遺伝子頻度になる）
-問2: 上で作成したdrift()というユーザー定義関数を、次の２つの場合で実行しなさい。
1. 実験繰返し10、1集団での観察世代数500、集団サイズ（遺伝子数）200, 対立遺伝子ａの数の初期値100
2. 実験繰返し10、1集団での観察世代数100、集団サイズ（遺伝子数）20, 対立遺伝子ａの数の初期値10
いずれの場合についても、できたグラフをWinShotで切り取って、課題提出ページに貼り付けなさい（注：1番目のシミュレーションの実行には、数分かかるかもしれません）（注：あまり大きい図（画面全面に広がる図とか）にはしないでください）
-問3: 上の２つのグラフを比較して分かることを述べなさい（集団遺伝学の問題）。
-else以下の処理が必要無いときは省略できる。また、1行に書くことも可能。
i=3 #iに3を代入
if(i==3){ print("三") } #iが3ならば三と表示
**繰り返しと条件分岐 [#qa2abd24]
では、上の繰り返しと条件分岐を組み合わせて使ってみよう
練習問題：for命令を使って1から10までの数字を表示させなさい、5だけは"five"と英単語で表示させなさい
--ヒントを穴埋め式で示しておく。
for(i in _:_){ # 10回繰り返し。iは1から10まで変化
if(i == _) { #iが5の場合
print("five") #fiveと表示
} _ { #それ以外の場合
print(_) #iの値を表示
}
}
**簡単なシミュレーション： 繰り返し・条件分岐・代入でランダムウォークを表現 [#ve13052d]
これまでやってきたプログラミングだと、まだ、「コンピュータに何かをさせている」ということが、あまり実感できないかもしれない。そこで、第13回から本格的にとりくむシミュレーションのさわりを、時間の許す限り解説しておこう。とりあげるのは「ランダムウォーク」で、次に移動する位置が確率的に無作為に決まるというものだ。~
例えば、葉っぱの真ん中にポトリと落ちた虫が、全く無作為に葉っぱを食べ進むとき、食べ後はどのようなパターンを描くだろうか？~
自分がすでに通った場所にも無作為に戻ると仮定して、ランダムウォークのシミュレーションを作ると、次のようになる。~
ここで想定しているのは、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個の整数乱数を発生させる
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の正方形の空間をはみ出すには、何回ぐらいかかるかを予想してからやってみると面白い。~
また、ここにあげたのは非常にシンプルな条件でのランダムウォークのシミュレーションだが、虫が実際に葉っぱをたべる様子をシミュレートする場合、葉の縁から食べ始めるとか、一度たべたところには戻らないとか、自分が落ちてしまわないように食べるなどの条件を付け加える。シミュレーションのプログラミングは、単純なモデルからはじめ、徐々に複雑にして行けば良い。~
**第11回授業の課題 [#q9067699]
-提出期限：''7月6日水曜正午''
--moodleページ：http://bean.bio.chiba-u.jp/moodle から提出して下さい。
**第12回授業の課題 [#zec6bcd4]
-提出期限：''7月13日水曜正午''
***課題１． 意見調査 [#u13a8f4f]
-2点満点で評価されます。http://bean.bio.chiba-u.jp/moodle から必ず提出して下さい。
***課題2：コインを投げる実験のシミュレーション [#o13d2c1c]
http://bean.bio.chiba-u.jp/moodle から提出して下さい。
***課題3:（発展課題:加点有り）自作シミュレーション [#q390e6d0]
http://bean.bio.chiba-u.jp/moodle から提出して下さい（ プログラミングをもっとやってみたいという人のみ提出）。
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")