このページは編集中です †

↑

Rを使ったプログラミング †

↑

【トピックス】帰無仮説・有意水準・有意差・仮説検定。。。が分からない †

前回授業でやった仮説検定のところがよくわからなかったという意見が多かったです。そこで、正確さは二の次にして、仮説検定で行っている作業をイメージで理解してもらえるようまとめてみました。それぞれどこの部分が、帰無仮説、有意水準、帰無仮説の棄却にあたるのか考えてみてください。また、なんとなく理解できたら、統計の教科書を読んでみることを強くお勧めします。

snmpenwl.gif
WS000002.JPG

1枚のコインを20回投げたら、オモテが5回しかでなかった。このコインはいかさまコインかなー？


もし、このコインがいかさまコインじゃない（正常なコイン）とすると...

20回中5回しかオモテが出ないのっていうのも、わりとフツーに起こることなんだろーね。


snmpensl.gif
Untitled-7.gif

わりとフツーに起こるこってどうやって数字で示す？確率で 95パーセントぐらいかな?...
正常なコインを20回投げた時のオモテの数ごとに、確率をグラフで表してみると...




snmpengi.gif

なんと！!95パーセント以上っていうのは、オモテの数が6枚から14枚までの場合だっ！

オモテの数が5枚の場合は、5%の範囲というフツーじゃ起こらないとした範囲に含まれちゃってる...


「20回中5回しかオモテが出ない」っていうのが、「正常なコイン」では滅多に起こらないっていうわけだから、

このコインが「正常なコイン」っていう前提が間違っていたと考える方がいいんだろうね。

うん、きっとこのコインは、いかさまコインなんだ


snmpenas.gif

でも、待てよ。フツーじゃ起こらないことの判断基準が5%って大きすぎない？
5%っていうと、、、100回中5回。。。 20回中に１回だ！
っていうことは、もしも「正常なコイン」という前提が正しくても、20回に1回は、間違った判断を下してしまう可能性があるってことだね。

↑

仮説検定：用語の整理： †

↑

前回やった２項検定の補足：どうしてp値が0.04139になるか？ †

> binom.test(15,20)
	Exact binomial test
data:  15 and 20 
number of successes = 15, number of trials = 20, p-value = 0.04139
alternative hypothesis: true probability of success is not equal to 0.5 
95 percent confidence interval:
 0.5089541 0.9134285 
sample estimates:
probability of success 
                  0.75

Untitled-8.gifUntitled-7.gif
file期待値.xls

↑

Rのおさらい †

　前回やったRを使った計算のおさらいを、サンプルデータをつかってやってみましょう。

↑

サンプルデータ：　仲岡先生の実習から †

下のデータは仲岡先生の実習でやった、調査コドラートの面積と、出現種数です

コドラート1    10*10    2種
コドラート2    15*15    4種
コドラート3    20*20    7種
コドラート4    25*25    7種
コドラート5    30*30    7種
コドラート6    35*35    9種
コドラート7    40*40    8種
コドラート8    45*45   12種
コドラート9   50*50   10種
コドラート10   55*55   10種
コドラート11   60*60   11種
コドラート12   65*65   11種
コドラート13   70*70   9種
コドラート14   80*80   10種
コドラート15  100*100  12種
（データ提供、してくれた方に感謝！）

　Rを使ってこのデータの、面積と種数の相関関係を知りたいとき、必要なのは、コドラートの１辺の長さと種数のベクトルです。

　今から何をするか、もうお分かりですよね？

予想通り、K2Editorを立ち上げて、要らないところを消し、

10   2
15   4
.    .

というデータに整形して、エクセルにコピペして、１列ごとのデータを scan() 関数を使ってオブジェクトに読み込みます。では、やってみましょう。（そろそろミミタコだろうから、説明は簡単に書いておきます）

↑

1. K2Editorでの操作 †

  1. 上の囲みのデータをコピペ
  2. 正規表現置換 
    検索文字列：　^.*\*　（説明：行頭から最初のアスタリスクまで）　　
置換文字列：＜無し＞（つまり削除するということ）

検索文字列：　種　　　　　　　　　　　
置換文字列：＜無し＞

検索文字列：  +　（説明：左に書かれているのは半角のスペースと半角のプラス）　
置換文字列： \t　（説明：タブ）
  3. 終わったら全体を選択してコピー。　エクセルを起動
↑

2. エクセルでの操作 †

  1. K2Editorで整形したデータをペースト
  2. ２列に分かれてデータが表示される 
    10  2
15  4
.   .
  3. 1列のデータをコピーして、Rに移動
↑

3. Rでの操作 †

  1. scan()関数を使って、コドラートの一辺の長さを qに、種数をyに代入 
    > q=scan()
[1]　　　　←　ここでペースト。yについても同様
> y=scan()
  2. （確認のため、それぞれ表示しておきます） 
    > q
[1]  10  15  20  25  30  35  40  45  50  55  60  65  70  80 100
> y
[1]  2  4  7  7  7  9  8 12 10 10 11 11  9 10 12
  3. 欲しいのは面積なので、xにqの2乗を入れる 
    x=q^2
  4. それぞれの常用対数を使ってグラフにプロット（対数グラフ用紙必要無し！ [smile]) 
    > plot(log10(x), log10(y))
  5. 回帰分析 
    > result=lm(log10(y) ~ log10(x))
> result
Coefficients:
(Intercept)     log10(x)  
   -0.1375       0.3250  
ここに表示された結果は、log(y)=0.325*log10(x)-0.1375 という回帰直線の式
  6. 回帰直線の描画 
    > abline(result)
  7. 解析のサマリー表示 
    > summary(result)
Call:
lm(formula = log10(y) ~ log10(x))
Residuals:
     Min        1Q    Median        3Q       Max 
-0.211429 -0.054045  0.006334  0.053616  0.142160 
Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.13750    0.15254  -0.901    0.384    
log10(x)     0.32498    0.04713   6.896 1.09e-05 ***　　←有意水準0.1%で有意
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 　←これ、有意水準
Residual standard error: 0.09927 on 13 degrees of freedom
Multiple R-Squared: 0.7853,	Adjusted R-squared: 0.7688 
F-statistic: 47.55 on 1 and 13 DF,  p-value: 1.092e-05
  8. 相関の検定 
    > cor.test(log10(y),log10(q2))
Pearson's product-moment correlation
data:  log10(x) and log10(y) 
t = 6.8957, df = 13, p-value = 1.092e-05
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.6847590 0.9618161 
sample estimates:
     cor 
0.886173 　　←　相関係数

　この授業では、それぞれの分析の意味については説明しません。仲岡先生に聞いて下さい。（グラフ用紙で提出しなくてもいいと仰っていたような、いなかったような。。。）

＊大事なこと＊

「こんな解析、自分でやるのは絶対無理！」と思うかもしれませんが、
 そんなことは無いです！

　上のような解析の方法、私だってRでできるとは知りませんでした。 では今回、私はどうやったのでしょう？「回帰分析行いなさい」なんて言われたら、どうすればいいのでしょうか？もしかしてRでできるかなと思ったら、

Adobe Readerを使って、マニュアルを検索します
↑

仮説検定：補足 †

前回アンケートで、仮説検定について「帰無仮説や有意水準などが分からない。」、「Rの操作は分かったけど、その意味が分からない」という意見があります。でも、私は、今の段階では、

全員がRを使って仮説検定の操作を行うことができただけで十分！

と考えています。生物学では収集したデータをもとに何らかの判断を下す場合、ほとんどの場合仮説検定が必要になります。これは別に系統学や生態学に限ったことではありません。最新号のNatureに載っている染色体不安定性を示すマウスの腫瘍に関する論文でもt検定が使われています。

↑

Rを使ったプログラミング演習 †

　さて、今回で３週目になるRという計算ソフト。コマンドを入力するのはちょと難しいけど、エクセルみたいにマウスを使って面倒な操作をしなくても、一瞬にしてグラフが描けてしまう優れものですね。これなら、良いレポートを仕上げるために、何度も解析をやりなおして、グラフを描き直すっていう作業も面倒ではなさそう。。。 　初めて使ったソフトなのに、皆さん、すでにRの操作には慣れてきているし（全員、課題ができていた）、興味も感じてくれたよう。そこで、今回は要望の多かったプログラミングをRを使って勉強します。

↑

プログラミングって何？ †

　ところで、プログラミングって何をすることでしょうか？IT用語辞典で調べてみると、、、

　　載っていません

　唯一近いところでは、「プログラミング言語」というのがあります。説明文からプログラミングに関わるところを抜き出して、「プログラミング」の説明文を作ってみると、

　今日の授業でコンピュータにあれこれ命令を与えるために使うRは、計算ソフトですが、インタプリタ型のプログラミング言語ということもできます。

　つまり、プログラミングとは、簡単に言うと、「言葉を使ってコンピュータが理解できる命令を作ること」　です。

↑

コンピュータにどんな命令ができるのか？　or　どんな命令がしたいか？ †

　「Wordを立ち上げて文章を作成し、印刷する」とか、「Excelで家計簿を管理する」とかは、市販のソフトウェアを使った作業ですが、いずれも、非常に高度で複雑な命令が、使用者からコンピュータに向かって発せられています。マウスを動かしたらポインタが動くなんていうのは、けっこう複雑な命令ですよ。でも、ほとんどの場合、一般の人がコンピュータでやりたいと思うほとんどの作業は、わざわざ「命令」なんていうものを意識する必要は無いはずです。利用者が理解しやすく・使いやすいように、専用のソフトウェアが開発され、販売されているでしょう。

　生物学の研究者でも、プログラミングなどを経験したことの無い人は大勢います。だって、今や、DNAシーケンスの決定や、整列や、系統樹作成にだって、専用のソフトウェアがあるのですから。。。それぞれのソフトウェアの使い方を覚えれば、それですむわけです。私たちも新しい解析方法が開発されるたびに、そのソフトウェアの使い方を覚えるのに苦労しています。

　でも、それだけで、本当にいいのでしょうか？良くないという点が、まず、２つ思いつきます。

1. 自分の目的にあったソフトウェアがいつも存在するというわけではない
2. ソフトが存在する場合、目的の数だけ、ソフトの使い方を覚えなければならない

　ということです。特に１番目は致命的ですね。ソフトウェアが開発されてなかったら、目的が達成できないのですから。例えば、集団遺伝学で学ぶ、遺伝的浮動の解析ソフトなんていうのは、パソコンショップに行っても売ってないです（でも、じつは、このくらいメジャーなものなら、ネットで探せばフリーのものがあります）。

　そこで、この授業でプログラミングを勉強する理由は次の３点です。

  1. 市販のソフトではカバーできないような目的を、コンピュータに命令を与えることで達成する
    • あなたのコンピュータは実は、すごい性能をもっています。その性能はできる限りつかわないと損
  2. プログラミングを行うことで、目的とする生物学的現象の理解が深まる
    • コンピュータは石頭なので、物事を本当に論理立てて説明してやらないと、命令を聞いてくれません
      　　プログラミングの過程で、生物学的現象の意味を論理立てて考える必要に迫られます
  3. プログラミングを行って、命令を実行させることで、初めて、コンピュータに作業をさせていることが実感できる
    • ワープロや表計算とかは、どちらかというと、道具を使って作業をさせられている感が強いですよね
↑

初めて(?)のプログラミング："Hello World!" †

「Hello World!」だなんて、なんか変なタイトルだなーと思いますよね。誰が使い始めたのかは知りませんが、様々なプログラミング言語の教科書で、最初に出てくるのがこのプログラミングです（プログラムを勉強したことのある人なら、たいてい誰でも知っています）。では、早速やってみましょう。Rを立ち上げて、次のように入力してください。

print("Hello World!")

うまくご挨拶できましたか？

ここでやったのは、print()という関数を使って"Hello World!"と画面に表示させただけですが、これも

画面にHello Worldと表示させる

ということを目的としたプログラミングです。

えっ？それなら、Rの最初でやったオブジェクトの内容を画面に出すのもプログラミング？

と聞かれるかもしれませんが、その通りです。次の囲みの中を全てコピーして、Rのコンソールにペーストしてみましょう。

x=c(1,2,3,4,5,6,7,8,9,10)
sum(x)

上でやったプログラミングは、

画面に1から10までの整数の合計を表示させる

というプログラミングです。皆さんがこれまでにやってきたRを使った操作も、何らかの処理の結果を画面に表示させるプログラミングだったわけです。なんとなく、プログラミングが身近になったでしょ？

↑

複数行の命令文 †

上でやった２つのプログラミングですが、違いがありますね。"Hello World!"の方は１行だけの命令だし、合計値の方は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(log10(q^2), log10(y))

このプログラムは、

２つの数値ベクトルについて計算を行って、対数グラフを描く

というものでした。さっき、上でやったのと同じことを、複数行にまとめて書いて、まとめてRに命令したわけです。

↑

Rで使う命令 †

　授業では一つ一つ解説しますが、復習課題をやるときとか、自分で新たな解析に挑戦するときは、どういう命令があるかが書かれた説明書が必要になるでしょう。

↑

プログラミングの基礎: 繰り返しと代入と条件分岐 †

上で複数の行で様々な命令を一度にＲに与えられるということを解説しました。でもやったことは、それぞれ個々の計算を順番に並べて、一度にペーストしただけですよね。そこで、次にもっとプログラムらしい命令をRに与えてみましょう。それは、繰り返しと代入と条件分岐というもので、プログラミングの基本中の基本です。

↑

繰り返し †

人が不得意でコンピュータが得意なことの一つは、単純な繰り返し作業を際限なく（文句も言わずに）行うことです。要するに、プログラミングの一つの目的は、そういう面倒な繰り返し計算をコンピュータにやらせることです。では、何回繰り返すかというのをどうやってコンピュータに指示するかというと、forという命令を使います。

↑

代入 †

上の for(i in 1:10) { print(i) } というプログラムでは、iの値は繰り返しの度に違う値になります。

1回目  iの値は1
2回目　iの値は2
3回目　iの値は3
....   ........
10回目 iの値は10

では、一つ前の繰り返しの時のiの値を使いたい場合はどうすれば良いでしょうか？例えば、1から10までの合計をfor命令を使ってプログラミングすることを考えると、

今の繰り返し回数のiの値に１つ前回数までの合計を加えて表示させたい
1回目  iの値は1   １つ前までの合計は0  1 を表示
2回目　iの値は2   １つ前までの合計は1  3 を表示
3回目　iの値は3   １つ前までの合計は3  6 を表示
....   ........
10回目 iの値は10  １つ前までの合計は45 55 を表示

こういう場合に代入という命令を使います。その命令は

=

です。通常の計算では'='というのは「3+5=」というように、計算を行わせる意味で用いられますが、プログラミングにおいては、

左辺の変数に右辺の計算結果を代入する

という意味で用いられます。上の計算のプログラムは次のようになります。

goukei=0           #goukeiという変数に初期値0を代入
for(i in 1:10) {   #以下を10回繰り返す
  goukei=goukei+i  #goukeiに前回までのgoukeiの値にiの値を足したもの代入
  print(goukei)    #goukeiの内容を画面に表示
}

goukei=goukei+iという命令文では、右辺の「goukei+i」を先に計算して左辺の「goukei」に代入しているところが、分かりにくいかもしれません。

練習問題：　1から10まで累積してかけ算をした結果を全て表示させるプログラムを作りなさい
↑

条件分岐 †

プログラミングの基本技の最後は１つは条件分岐です。if()を使って表現します

i=3                     #iに3という数値を代入
if(i==3){print("三")}  #iの値が3ならば画面に三と表示

この例の２つめのif命令で、iの値を評価しています。評価に使われるのは比較演算子というもので、

==　    等しければ
!=      等しく無ければ
>, >=   左辺が右辺よりも大きければ、左辺が右辺以上ならば
<, <=   左辺が右辺よりも小さければ、左辺が右辺以下ならば

ということを意味しています

↑

繰り返しと条件分岐 †

　では、上の繰り返しと条件分岐を組み合わせて見ましょう

練習問題：for命令を使って1から10までの数字を表示させなさい、5だけは"five"と英単語で表示させなさい
↑

Rを使ったシミュレーション：「遺伝的浮動」を実感してみよう！ †

drift.gif

上で説明してきたことで、プログラミングの基本は終わりです。もう皆さんは、繰り返しと条件分岐を使って、かなり複雑な処理ができるはずです。そこで、いよいよ実践に移りましょう。

右の図は見たことはありますか？これは、遺伝的浮動のシミュレーションの結果で、集団サイズN=100, 対立遺伝子が２つのとき、片方の対立遺伝子の頻度の変化を表したものです。。最初の世代の対立遺伝子の数は50:50。このグラフを使って、遺伝的浮動とはどういうものか、説明できる人はいますか？また、このグラフはどうやって作ったものですか？

　こんな図を作るのはとうてい無理と思っちゃうかもしれません。でも、Rを使ってプログラミングすると、わりと簡単に描けてしまいます。

　そこで、ここから先は実践編として、「遺伝的浮動」のシミュレーションをやることにしました。このテーマを選んだ理由は、「生物学科の皆さんになじみが深い」、「Rをつかってやってみて見た目が面白い」、「思ったよりはプログラミングが簡単」というこです。さらに、私にとっては初めて自分で挑戦したシミュレーションで、できあがったグラフを見たときにはかなり嬉しかったものです。　

　なお、シミュレーション（simulation, 「シュミレーションでは無いですよ！」）は、数学モデルなどを用いて現実に似た状況をコンピュータ上で再現することです。模擬実験とも呼ばれています。（大量の計算を必要とするので、Rのようなインタプリタ型のソフトでは、大規模なシミュレーションは難しいです）

↑

0. シミュレーションの前段階：モデルの想定 †

　これからシミュレーションを行うのですが、つぎのような状況を想定してみましょう。

mona2.gif
mona1.gif

まず、対象とする架空の生物（仮に、「モナー」と呼びます（注：モナーについて詳しく知りたい人は、ウィキペディアで「モナー」を検索してみてください。））は2倍体の生物で、雌雄の区別が無く、繁殖期になると集団内の個体が一斉に配偶子を数千億個ずつ放出して死亡します。配偶子は生息域の空間内をランダムに漂った後、24時間後に、最も近くにある配偶子と接合し、新しい繁殖個体になります。集団の生息域は空間的に限られているため、親世代と同数の子供個体しか、成熟個体になれません。どの子供が成熟個体になるかは、ランダムに決まるものとします。つまり、このモデルだと、世代間で集団サイズに変動は無いということです。モナーにはトンガリ耳と丸耳の２タイプが1対立遺伝子(A, a)で決まることが知られています。トンガリ耳の遺伝子型は(AA, Aa)、丸耳の遺伝子型は(aa)です。

なお、上のモナーの説明は、私が今回のシミュレーションのために、適当に設定したものです。だって、「任意交配をする有限サイズの生物集団から遺伝子をランダム抽出して、同サイズの次世代集団を作った」なんて言い方をするよりは、生物を扱ってる感じがするでしょ？

↑

1. 初めてのシミュレーション：10個体からなる集団における、遺伝子頻度変化のシミュレーション †

　モナー10個体からなる、１つの集団を考えてみましょう。いまこの集団で、対立遺伝子Aの頻度が0.5, 対立遺伝子aの頻度が0.5のとき（10個体だから遺伝子の総数は、20個です。そのうち半分の10個がAで、10個がaです）。集団中の配偶子の接合はランダムに起こると仮定して、100世代の間に、集団中の対立遺伝子の頻度はどのように変化するでしょうか？

slide5.gif slide6.gif

今から実験をするのですから、前もって予想しておいた方が面白そうです
この集団の遺伝子頻度が1か0に固定するまでに、何世代ぐらいかかると思いますか？
0-25  26-50  51-75  76-100  100<

それでは、次のプログラムを走らせて、20個の遺伝子中対立遺伝子aが10個ある任意交配集団の、100世代の間のaの遺伝子頻度の変化をみてみましょう。

numa=10                     #対立遺伝子aの集団中の個数
for(i in 1:100){            #この行から対応する}までの間を 100世代分繰り返す
 counta=0                        #次世代にランダムに残るaの値を入れる変数を初期化
   for(j in 1:20){               #数千億からなる配偶子プールから1個の配偶子を取り出す試行を20回行う
     if ( runif(1) < numa/20  ){    #乱数を一つ発生させ、それが対立遺伝子aの頻度(a/20)よりも小さければ
       counta=counta+1              #aの個数を入れるcountaという変数の値を１つ増やす
     }
   }
 out=c(i, counta/20)             #何世代目かを表す数字(i)と対立遺伝子aの頻度(a/20)をoutというオブジェクトに入れる
 print(out)                      #outの内容を画面に表示する
 numa=counta                     # aの値を新しい世代の対立遺伝子aの個数(counta)で置き換える
}
↑

遺伝的浮動のプログラムをステップ・バイ・ステップで理解する †

  1. まず、前回説明した「モナーの生活環」をみて、シミュレーションの条件を自分で考える。最初のシミュレーションなので、簡単に、 
    20個の遺伝子からなる集団
0世代目の対立遺伝子aの数が10（つまり対立遺伝子aの頻度は0.5）
この集団を100世代観察する
    という条件にしておく
  2. 今、上に書いた数字に自分で考えた変数名をつける。名前を付けておくと、あとでコンピュータに命令すときに便利だから。数値が決まっているものは代入しておく。名前は、Rの内部で使われているもの以外なら、何でも良い。自分に分かりやすいようにつける。 
    ＊自分で考えた変数名と定数
numa=10   #0世代目の対立遺伝子aの数。2世代目以降は、シミュレーション実験で抽出されたaの数が代入される
counta=0    #配偶子プールから、次世代の集団に取り出された対立遺伝子aの数。１世代が終わったあとで、
             数字が入るので、とりあえず今は0を入れておく
20  #集団の大きさ（つまり、集団に含まれる遺伝子の数）
100  #繰り返しの世代数
    slide10-1.gif slide10-2.gif
  3. 変数名が決まったら、プログラムの構造を考える。紙に絵をかいて、どういうプログラムなら、「モナー」の生活環を模倣できるか考えたら、言葉にしてみる。 
    ＊自分で考えたプログラムの構造。コマンドを使って書かなくても、
  とりあえず、構造が理解できるように書いてあればよい。
    slide10-3.gif
    1. 何世代分か、以下の操作を繰り返す
    2. 　１つの世代の中では、配偶子プールから配偶子を１つ選ぶという操作を20回（遺伝子の数）繰り返す
    3. 　　　取り出した配偶子が対立遺伝子aだったら、変数countaの値を一つ増やす。そうでなかったら何もしない
    4. 　集団内で、20回配偶子を取り出す繰り返しが終わったら、対立遺伝子aの数を numaに代入する
    5. 次世代の集団における対立遺伝子aの頻度(num/20)をプリントする
  4. プログラムの構造ができたのだから、次ぎに、一つ一つのステップを、コンピュータが理解できる文で表現してみる。
    1. まず、100世代この集団を観察するのだから、繰り返し文で100回繰り返す構造が必要になる。 
      ＊100回繰り返す繰り返し文の構造
for(i in 0:100) {
      ＜※ここに繰り返される内容が入る＞
}
    2. １世代で行われる操作（上の※の部分）は、１つの集団から配偶子を親の集団と同数（つまり、20個）取り出すという繰り返し操作。なので、上の100回繰り返すfor文の中には、さらに、集団中の遺伝子の数だけ繰り返すfor文が入る 
      ※集団中の遺伝子数分、つまり、20回繰り返す文
   for(j in 0:20) {
   　　＜※※ここに繰り返される内容が入る＞
   }
    3. 上のfor文の中では、配偶子プールからpopsize個取り出した配偶子の中に、対立遺伝子aがいくつ含まれているかを数えること。そのとき、 
      　※※の操作：
  「配偶子を一つ取り出すとき、遺伝子頻度の確率で、対立遺伝子aをひく。もし、対立遺伝子aを引いたら、
    aの数を入れておく変数(counta)の値を１つ増やす」
      　 ＝　　slide10-4.gif　　　slide10-5.gif 
        上の『対立遺伝子aを引いたら』（このことは対立遺伝子頻度aの確率で生じる）というのは、
   『乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら』（これも対立遺伝子頻度aの確率で生じる）と
　　同じことと考える。そこで、
　「乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら、
    aの数を入れておく変数(counta)の値を１つ増やす」
      • ここが理解されなかったら、Rを使ってrunif(20)を表示させる。0.5より小さい値の数を数えて、次世代の集団の対立遺伝子aの頻度とし、もう一度Rを使って乱数発生。この操作を3回ぐらいやる
    4. 上の操作を、20回繰り返したら、countaという変数には、取り出された対立遺伝子aの数（countaという名前にしてある）が入っているはず。この値は、次世代の集団における対立遺伝子aの数なので、次世代（つまり、外側のfor文の次の繰り返し）に渡さなければならない。そこで、numaという変数にcountaの値を代入して、次の繰り返しに渡す。countaという変数は、次の繰り返しでまた、0からaを数えなければならないから、ここで0を代入して、初期化しておく 
       ・内側のfor文の繰り返しが１回終わったら
  numa = counta
  counta = 0
    5. 次の世代でも上と同じことを繰り返すが、集団内の対立遺伝子頻度は、前の世代の対立遺伝子頻度で決まる。
↑

第10回授業の課題 †

↑

課題１．アンケート調査 †

  1. http://bean.bio.chiba-u.jp/joho19/ に、「自分のID」/10　という新しいページを作成し、下の囲みの中にあるアンケートをコピー・ペーストして、「回答：」の後に答えを書き込むこと。
↑

課題１．意見調査 †

　下の囲みの中にあるアンケートをコピー・ペーストして、「回答：」の後に答えを書き込むこと。

*第10回授業アンケート
**氏名：
**課題への回答
-今日（6月22日）の授業の進み方は？（はやい、丁度いい、おそい）
--回答：
-今日の授業の難しさはどう感じましたか（簡単すぎ　簡単　丁度いい　難しい　難しすぎ）：
--回答：
-難しいと答えた人は、特にどの点が難しかったですか？：
--回答：
-今日の授業は（よく分かった　分かった　分からなかった）：
--回答：
-分からないと答えた人は、特にどの点が分からなかったですか？：
--回答：
-今日の講義で理解できなかった用語があったら挙げてください：
--回答：
-Ｒで何ができるのか、だんだん分かってきたと思います。今後自分の勉強にＲを使いたいですか？
--回答：
-Ｒでも、Ｒ以外の言語でも、何かやってみたいプログラミングはありますか？
--回答：
-今回やったプログラミングはおもしろかったですか？
--回答：
↑

復習課題：Rによるプログラミング：繰り返しと条件分岐 †

次の問題文の指定にあるプログラムを作成しなさい。作成したプログラムは課題提出ページに貼り付けなさい。 行頭には必ず半角の空白を入れること。

  1. 1から20までの総和を計算するプログラム（?のこと）（注：繰り返し文(for())を使ってプログラムを書いてください。sum()関数は使わないこと）
  2. 　乱数を20個発生させて、0.5未満のものだけを画面に表示させるプログラム
↑

予習課題：Rを使った関数の定義 †

「関数の定義」なんていうと難しそうですが、ようするに、作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにすることです。例えば、

1から入力した数値までの全てを横一列に表示させるプログラムを作り
displayという名前の関数として定義する

どうやるかというと、

display=function(a){       #関数定義の始まり
 kekka=c()                 #kekkaに空ベクトルを代入して初期化
 for(i in 1:a){              #a回（iの値を1からaまで変化させる間）繰り返し
   kekka=append(kekka,i)     #kekkaというベクトルにiを要素として代入
 }
 print(kekka)              #kekkaの内容を表示
}                          #関数定義の終わり

上の方で作ったプログラムの上と下を関数をfunction(){}で囲んで、オブジェクトに代入するだけですから、簡単です。この関数を実行するには、最初aの値を()の中に入力して、

display(10)

とします。 ↑を使って何回も実行してみると、結果がいろいろ変わるのがわかりますよね。

そこで問題です。上の関数の定義方法に従って、入力した数までの合計値を計算するsumupという名前の関数を作成します。下の囲みの中の_の部分（１文字に対応するとは限らない）を埋めて、プログラムを完成して提出しなさい。

 sumup = ______       #関数sumpuを定義
   kotae = _          #kotaeを初期化
   _ (i in ___){      #1からaまで繰り返し
     kotae = _____    #kotaeにiの値を足したものとkotaeに代入
   }                  #繰り返し終了
  _____               #kotaeを表示
  }                   #関数定義の終了