複雑システム系演習１（鈴木泰博 担当分）

第１回

ガイダンス

鈴木（泰）が担当する後半部では，数理モデル化とシミュレーション，の演習をおこないます．この演習ではプログラミング言語「R」をつかいます．４回分のコース内容は以下の通りです．

1. R言語・入門
2. 数理モデル化入門（ルール力学系によるモデル化）
3. 確率過程によるモデル化
4. R言語によるセルラーオートマトン

評価　以下の２つ課題から総合的に評価します．

i) 毎回行った”履歴”をそのままWORDなどに貼り付けて提出する．

ii) 最終課題：自分で数理モデルを作成し，シミュレーションを行い，結果をまとめる．

最終課題の評価基準

以下の評価基準で評価を行う．C以上が60点以上（つまり鈴木担当分での合格点）である．

コースのすすめかた

コース資料に課題が出題されているので，それを解いてく．課題を行った履歴をすべてコピーして，WORDなどに貼り付け，NUCTから提出する．

最終課題はWORDなどをもちいてレポートを作成する．締め切りは第４回の次週の木曜日の午後５時とする．

R言語によるシミュレーション入門　鈴木泰博

本ドキュメントの参考文献

渡辺利夫　誰にもできるらくらくR言語、ナカニシヤ出版

樋口千洋　Rによるバイオインフォマティクスデータ解析（第２版）共立出版

高階知巳　プログラミングR　オーム社

以下では"＞"（たとえば＞1+2）はRのコンソールのプロンプトをあらわします。なので、**[＞は入力する必要はありません]。また、以下ではコンソールではなく、Rエディタを用いて作業することを学びます。

Rのインストール．

筑波大学のミラーサイト等からダウンロードできる (http://cran.md.tsukuba.ac.jp)。

インストールの方法

1. 筑波大学のサイトにアクセスする

2. Download R for Windowsをクリック

3. baseをクリック

4. 最新版をダウンロード

5. 「このファイルを実行する」か「ファイルをコンピュータに保存する」かの選択を尋ねられる→保存する　を選択し、ディレクトリをきめて「実行」をクリック。インストール中に使用する言語としてJapaneseを選択して「OK」をクリック

6. 「R for windows

   ...」セットアップウィザードの開始のウィンドウが開くので「次へ」をクリック

7. 「GNU F_2NERAL PUBLIC

   LICENSE」を尋ねてくるので、内容を確認して「次へ」

8. 「インストール先の指定」を尋ねられるので「C:\Program Files\R\R-（ここにはバージョン名がはいる）」を選択し「次へ」

9. あとはすべてに「次へ」を選択

10. 「完了」が表示されたら「完了」をクリック

11. デスクトップ上にショートカットアイコンが表示されていることを確認

計算•アルゴリズム•シミュレーション

本実習では現象の計算モデル化とシミュレーションの仕方について学びます。モデル化とシミュレーションは（複雑系）科学一般で用いられています。でも、やにわに「はい、ではモデル化しましょう。」と云われると。。。戸惑います。なので、"他の誰かがやったモデル"を理解して、それを少し修正して。。。となりがち。それはそれで、研究のひとつだし、よい勉強の仕方でもあります。でも、"そもそもモデル化やシミュレーションとはなにか"について、イメージがあったほうがクリエイティブになりますし、理解の助けにもなるとおもいます。

モデル化とシミュレーション

似ているような、全く違うようなイメージをもつかとおもいますが、どうなんでしょうね。。似ているのかもしれないし、なんか違うかもしれないし。モデル化ってナンダロウ...教科書なり論文を拡げると微分方程式や確率過程、その他さまざまなモデルがあって数式やフローチャートがあってナンダカヨクワララナイ。そして、本なり論文なりを閉じる。

　モデル化とは「モノコトの見方です」。誰にとっての？それはもう..."あなたにとってのモノコトの見方"です。えー！でもそんな、むつかしい数学とかプログラミングとかできないし...いーんです！ 例えば...履修登録をするとき、もしあと12単位とれないと留年したり卒業ができない...なんてちょっと切羽詰まった状況になるとすると、情報を集めに集め「すんごい厳しい単位がとれるかとれないか微妙な科目」とか「（いわゆる）楽勝科目」などを選別し、12単位をとるためには...と「計算」します（この場合の観測系はなんでしょうか？）。ここでの「厳しい科目」、「楽勝科目」こそが"モデル化"です。そして、厳しい科目の場合はとれるかわからないから勘定に入れず、楽勝科目をカウントして「絶対に12単位がとれる（であろう）」履修登録を行う（かもしれません。そんな"モデル化"をせずに興味と関心、内容を踏まえて履修するというのも一つの"モデル化"です）。ここでのモデル化は、自分とは能力も関心も異なる他の誰かのモデルは余り役に立ちません。大学院に進学し理論物理や数学を専攻しようとしている人の"モデル"と、文化やアートに興味もあり将来は美術系に進みたいと考えている人の"モデル"は...違いますよね。

Rの基本操作

起動と終了

アイコンをクリックするとR Console画面が表示される。するとプロンプト＞が表示され、入力待ちの状態になる。このプロンプトの右にコマンドを入力してEnterを押すとコマンドが実行され結果が表示される。

＞ 1 + 2

＞ [1] 3

が表示される。ここで[1]は行め出力であることを示している。出力数が多い場合には、ウィンドウに収まりきらずにデータが"折り返されて"しまう。

演習

①＞letters

とするとアルファベット26文字が表示されることを確認する。

②ウィンドウのサイズを広げて（できるだけ横幅いっぱいに広げる）＞lettersとする。そして、[ ] の表示がどのように変わるかを試してください。

コマンドの補完機能pl + TAB(２回)

作業をしていると、コマンド名がうろ覚えなことがある。そんなときに、この補完機能が便利である。

演習: plからはじまるpl... とのコマンド名をど忘れしてしまったので、この補完機能を用いてplからはじまるコマンドをリストアップしてみる

①＞pl

ここでTABを２回おす。

すると以下のようにリストアップされることを確認する。

以上の補完機能はディレクトリにおいても有効である。たとえば、list.filesコマンド（ディレクトリ内のファイルを表示させる関数）をつかうときに、表示させるディレクトリ名がはっきりしない場合、たとえば現在のディレクトリ"./"にあるファイルを表示させたい場合には、

演習

①＞list.files("./ TABを２回押す

②./にあるファイルが表示される

入出力

Rは大文字と小文字を区別する。

演習 以下のコマンドを入力し実行結果をくらべなさい

① ＞letters

② ＞LETTERS

③ ＞Letters

④ ＞1+2

⑤ ＞1 + 2

⑥ ＞ 1+ 2 - 3

⑦ ＞　f

⑧ ＞ 3 -

1

⑧について。関数の引数がたくさんあったり、式が長くなってしまったりした場合にウィンドウにはいりきれない場合があります。そのときには、**[適当に改行してしまってかまいません（読みやすいように）]、式などの途中で改行した場合には[プロンプトが ＞から+に]{.underline}かわります。式の最後までを入力すると、Rはそれを認識してコマンドを実行します。

⑧の場合には　3 -　と3ではじまる減算が途中で改行されているので、プロンプトが+にかわっています。なので１を入力すると３−１となってコマンドが実行され[1]2となります。

演習

1. ＞list.files("./ としてEnter
2. + 1+2などのコマンドを入力してみる

*list.files("ディレクトリ名")はディレクトリにあるリストをファイルするコマンドです。

①をするとプロンプトが+に変わってからコマンドの入力をうけつけてくれなくなってしまいます。ここから"脱出"して、プロンプトを＞に戻して、再びコマンド入力を受け付けられるようにしてください。そして＞に戻すことができた理由も述べなさい。

困ったときの強制終了

コンソールにあるSTOPボタンをおすとプログラムが強制終了される。

終了する場合には

＞q( )

とする「作業スペースを保存しますか？」と尋ねられる。Rではそれまでの作業で定義したオブジェクト（Rでは変数のことをオブジェクトとよぶ）を保存することができ、作業スペースを保存しておくと次に作業をはじめるときにオブジェクトを再定義する必要がない。

　後に作業を続けるならばスペースを保存し、そうでない場合には保存しないほうがよい（新たに作業をはじめるときに古いオブジェクトの定義が残っていないように）。

加減乗除〜演算子

加減乗除の演算は, +, -, *, / の演算子により行う。演算のための括弧は (, )の使用のみができる。{,}や[,]はR言語では別の目的につかわれる。

演習　以下の演算を行ってください。

1. 2 + 3, ②2 -- 3, ③2 * 3, ④2/3, ⑤2+3*(4 + 2), ⑥2 + 3/(4 + 2),

⑦2 + 3 / ((4 + 5) * (1 + 6)) ⑧2 %/% 1 ⑨2%% 1

(%/%は除算の整数部、%%は除算の余り)

このように加減乗除を行ったり、除算の整数部や余りを求めたりするために用いているものは"演算子"とよばれる。他のよくつかわれる演算子として列挙の演算子「:」がある。「:」は、"列挙開始の値 : 列挙終了の値"として用いられる。

演習 以下の演算を行ってください。

i) ① 1:1, ②1:10, ③1:10-4, ④1:10+4, ⑤1:(10-4), ⑥1:(10+4)

ii) ③と⑤、④と⑥を比較して、これらの演算の違いについて述べなさい。

オブジェクトと代入

Rでは変数のことを"オブジェクト"とよぶ。オブジェクト名は、半角英数字、「.」、「_」のみ使用できる（例: ys, y.suzuki, x_1など）。また"大文字と小文字は区別される"よって、ABCとabcは異なるオブジェクト名となる。

オブジェクトへの値の代入は -> により行う。->は -- と>を組み合わせたものである。-> は左辺から右辺への代入を、<-は右辺から左辺への代入をあらわしている。

オブジェクトaに1を代入するには、

a<-1もしくは1->a

となる。オブジェクトに代入されている値を確かめるには、オブジェクト名をそのまま入力してEnterキーをおせばよい。

＞ a

[1]1

演習

① aに1をbに2を代入し、a + bの結果をcに代入してcの値を表示させなさい。

2. cの値からaを減算してdに代入して、dの値を表示させなさい

3. e <- lettersとアルファベットのベクトルをeに代入して、eの値を表示させなさい

   lettersコマンドとは全アルファベットを要素とするベクトルをつくるコマンド

システム関数

R言語では、データ解析用にさまざまな関数がシステム関数として用意されている。

ベキ乗 ^

＞2^3

＞2^0.5

Rでは指数として少数をとることができるよって2^0.5のように平方根を計算することができる。

平方根sqrt

＞sqrt(2)

Rでは平方根の計算のためのシステム関数もあり、これを用いても平方根を計算することができる。

$e^{x}$はexp(x)$\ \log_{e}x$はlog(x)、対数関数で底が10の場合はlog10(x)のほか三角関数や逆三角関数も用意されている。また、四捨五入用の関数としてroundが用意されている。この関数はround (x, a)と定義されxには数値が、aには四捨五入したい少数の位がはいる。例えば、1.4145を少数第４位で四捨五入して少数第３位までもとめたい場合にはround(1.4145, 3)とすればよい。

演習 次のシステム関数を用いた計算を行いなさい。

1. 2^3 ② sqrt(2) ③ log(2) ④ exp(2) ⑤ sin(0.5) ⑥ cos(0.5) ⑦ tan(0.5) ⑧round(sqrt (2), 2)

その他の便利な関数

データ解析用の関数ではないが、定義したオブジェクトをリストしてくれる関数ls( )は、新たにオブジェクトを定義する際などに便利である。

たとえば以下の実行例で、character (0) とはオブジェクトがなにも定義されていないという意味である。次にオブジェクトに3を代入し（x<- 3）、あらためてオブジェクトのリストを確認してみると、xがオブジェクトとして定義されていることが示されている。

＞ls () ...オブジェクトのリストの表示

character(0) ... 「定義されているオブジェクトがない」

＞ x <- 3 ...オブジェクトxに3を代入

＞ ls () ...オブジェクトのリストを表示

[ 1 ] "x" ...xがオブジェクトのリストにはいっていることがわかる

すでに使用したオブジェクトを消去したい場合、たとえばオブジェクトxを消去した場合には

＞rm (x)

とすればよい。

＞ls()

として確認してみれば、character(0)となっている。

演習

1. オブジェクトのリストを表示させなさい。
2. 新たなオブジェクトに3を代入しなさい。
3. オブジェクトのリストを表示させ、新たなオブジェクトが追加されているか確認しなさい。
4. オブジェクトxを消去しなさい。
5. xがオブジェクトのリストから消去されているか確認しなさい。

実行したコマンドの保存と呼び出し

（その１）

R言語で実行したコマンドを保存して、記録しておきたい場合がある（よくある）。行数が多い場合や、手続きが複雑な場合には、プログラムのファイルとして保存し、それを呼び出して処理を行うことが必要となる。その場合には、

1. コンソールの画面の保存したい部分を選択してコピーしてメモ帳などのテキストエディタにペーストして保存する。
2. ファイル名を「ファイル名.R」として保存する。拡張子に".R"をつけておくと、コマンドライン（コンソール）から呼び出すことができる。

保存したRファイルを呼び出すには、

＞source("/Users/prof_ys/Desktop/test.R")

のようにsourceコマンドで呼び出すことができる。もしくは、メニューからファイル→ソースを読み込む　でファイルを指定して呼び出すこともできる。

演習

1. メモ帳などエディタを立ち上げる。
2. y <- 2 + 3と書き込み、これをtest.Rとして名前を付けて保存する。
3. R consoleのメニュー→ソースを読み込む　でtest.Rファイルを読み込む。
4. ＞y としてyが5になっているか確かめる。

（その２）Rエディタを用いる

ファイル→新しいスクリプトを選択すると、Rエディタが起動する。

エディタ内で、例えば2 + 3と入力した場合、実行するには

1. コマンドを選択（複数行でもよい）
2. Ctrl -- rで実行され、実行結果はRコンソールに表示される。

練習　例にあげたRエディタの起動と実行をやってみてください。もし、こちらの方法が気に入ったら、以降はRエディタを用いて作業しても結構です（私はRエディタを使ってプログラミングすることが多いです）。

データ構造

ベクトルの演算

ベクトルとは２つの以上の数値をもとに構成される。例えば、ベクトルをx_1、その要素を(3, 5, 6, 8)とすると、Rではこれをc(3, 4, 5, 6)と、関数cを用いて定義する。

＞X_1 <- c(3, 5, 6, 8)

ベクトルが、２つの要素からなるのに対して、１つの要素からなる場合はスカラーとよばれる。ベクトルの一部を表示したい場合には, [ ]を用いる。X_1の２番目の要素を表示したい場合には

＞X_1[2]

[1]5

とする。X_1の"２番目と4番目"のように複数の要素を表示させたい場合、まず表示させたい位置をベクトルとして表現する（この場合であれば　c(2,4)）、そして、このベクトルを用いてベクトルの要素を表示させる：

＞X_1[c(2, 4)]

[1]5 8

演習

1. 1, 2, 3, 4, 5を要素とするベクトルX_2をつくれ。
2. X_2の１番目と４番目の要素を表示せよ。
3. X_2の1番目と６番目の要素を表示せよ。

ベクトルの要素の書き換え

ベクトルへあたらな値を代入することにより、要素の書き換えを行うことができる。

たとえば、X_1の2番目の要素を5ではなく9にしたい場合には、

＞X_1[2]<-9

とすればよい。すると、

＞X_1

[1]3 9 6 8

となっていることが確認できる。

ベクトルの要素の削除

ベクトルX_1について、[-1]とマイナスをつけて位置を表示すると、その位置の要素が削除される。

演習

1. 5, 4, 3, 2, 1を要素とするベクトルX_3をつくれ
2. X_3の1番目の要素を６に書き換えよ
3. X_3の２番目と４番目の要素をそれぞれ10と100にせよ（ヒント: 複数の要素を表示させる場合にベクトルを用いたことを参考にせよ）。
4. X_3から1番目と3番目の要素を削除したベクトルをX_4に代入せよ（ヒント③での計算の仕方）。
5. X_4を表示せよ。

ベクトル同士の演算

ベクトル同士での演算はスカラーでの演算とどうように加減乗除の演算子（+ -* /）により行うことができる。ベクトルの内積は%*%を用いる。

演習

①1 2 3 4 5を要素とするベクトルX_1をつくる

②5 4 3 2 1を要素とするベクトルX_2をつくる

③X_1とX_2の加減乗除を行う

④X_1とX_2の内積を計算せよ

ベクトルの比較

ベクトルの比較は以下の演算子によって行うことができる。また、これらの演算子はスカラー間の場合にも使用することができる。

== 等しい

<, >　大小

<=, =< 等しいか大小

&　対応するそれぞれの要素同士で論理積した結果のベクトル

&&　最初の要素同士で論理積

|　対応するそれぞれの要素同士で論理積した結果のベクトル

||　最初の要素同士で論理和

! 否定

%in% ２つのベクトルのTRUE / FALSEが一致している位置をかえす

which TRUEの位置をかえす

any 1つでも条件をみたしていればTRUE、そうでなければFALSEをかえす

演習

1. 1 2 3を要素とするベクトルX_1をつくれ
2. 3 2 1を要素とするベクトルX_2をつくれ
3. X_1 == X_2
4. X_1 > X_2, X_1 < X_2
5. X_1 >= X_2, X_1 <= X_2
6. any(X_1 < X_2)
7. T F T Fを要素とするベクトルX_1をつくれ
8. F F F Fを要素とするベクトルX_2をつくれ
9. X_1 & X_2
10. X_1 && X_2
11. X_1 | X_2
12. X_1 || X_2
13. !FALSE, !c(T,F,T)
14. X_1 %in% X_2
15. which(X_1), which(X_2)

ベクトルの包含関係

長さが同じ２つのベクトル, X_1とX_2について、X_1とX_2のそれぞれの要素を比較して、どの要素も「等しい、もしくは小さい」場合には[X_1はX_2に真に包含されると定義する。一方で、どの要素も「等しい、もしくは大きい」場合にはX_2はX_1に真に包含されると定義する。また、すべての要素が等しい場合にはX_1とX_2は等しい（X_1= X_2）と定義する。

演習

①真の包含関係にある2つのベクトル、一部のみが包含される(真の包含関係にはない)２つのベクトルの例を挙げなさい。

②真の包含関係にあればTRUEを、そうでなければFALSEを返す計算を、以上の演算子をもちいて示しなさい（動作の例も示すこと）。