ここで"ベクトル"とは、数学の教科書に出てくるベクトルをそのままさしているのではない（そのため"ベクトル"としている）。JavaやCなどの、いわゆるALGOL型のプログラミング−パラダイム（プログラムの仕方についての哲学／方針のようなもの）では、データや変数間の計算は「１つずつ行う」と考える。たとえば、「ベクトル(1,2,3,4,5)の要素を２倍する」、場合であれば、
f<-function(X)
{
listx<-NULL
for(i in 1:length(X))...length関数はベクトルXの長さを返す関数
listx<-append(listx, 2*X[i])...append関数は結果を追記する関数
listx
}
となる。この場合には、ベクトルXから１つずつデータが読み出されて２倍されてから、その結果をlistxに書き込んでいく。X=(1,2,3)とすると、最初はlistxは初期化されているので、
NULL
(2)
(2, 4)
(2, 4, 6)
と処理がなされていく。
これにたいして"ベクトル"での演算とは、入力されたベクトルX＝(1, 2,3)を「バラバラにせずにそのままデータの塊り•集まり（集合）として扱う」と考える。そして、
「データの集合」---処理→「データの集合」−処理→「データの集合」−処理→...
と「データ処理した結果を入力として処理を行っていく」と考える。ここでの"処理"を"関数"とおきかえると、
f(データの集合)→出力...処理されたデータの集合
となるので、この"出力"を入力としてg(出力)と新たな計算を行うのであるが、出力とはf(データの集合であるから)これは、g(f(データの集合)となる。よって、以上のように「処理したデータの集合を入力として処理を繰り返す」ことは、先述した関数型プログラミングの仕方である「合成関数による計算」と同じことになる。
マッピング
こうした処理は"マッピング"とよばれ、マッピングを行うための関数が"apply系"の関数である。apply系の関数には以下のようなものがある。
sapply(xs, f)： xsの各要素にfを適用して、その結果をベクトルでかえす
lapply(xs, f)：xsの各要素にfを適用して、その結果をリストでかえす
rapply(xs,f)：xsの各要素に対して再帰的に関数fを適用して、その結果をリストでかえす
#これら以外にも、mapply, apply, tapply, eapplyなどがある（詳細は参考文献[高階, p184など]を参考にされたい）。
このapply関数をつかうと、f_1(2, 4, 5, 3) → (f_1(2), f_1(4), f_1(5),f_1(3))は、
sapply(c(2,4,5,3), ho)
[1] 4 7 8 6
となる。
フィルタ
"ベクトル"によるプログラミングでは、マッピングの他に「フィルタ」と「還元」の手法が用いられる。還元とは２つの引数をとる関数、例えば、x + yのような関数、をベクトルに対して適用するための処理である。だが、R言語には還元の処理を行う関数は用意されていない。その理由は、R言語では、基本的に引数としてベクトルをとるように設計されているので、引数が２つの関数を使うような場合が少ないからである。そこで、ここではフィルタについて紹介することにする。
フィルタとは、入力された集合の各要素を調べて、あらかじめ指定された条件をみたすものだけを取り出して、出力する処理を行う。そのためにはベクトルや行列の処理のところでの"条件付き参照"を用いる。もう一度、復習してみると、ベクトルa=(1, 2, 3, 4)から2以上の要素を取り出すための条件つき参照は
＞a[a>2]
とすればよかった。
以上の復習をふまえて、ベクトルから"偶数のみ"を取り出す処理を行うことを考えよう。そのために偶数か否かを判定する関数is.evenを定義し、条件つき参照をもちいてフィルタをつくってみよう。
演習
ヒント: "偶数"とはxを２で割った余りが0と考えられる。除算の余りを求める演算子は %% (xを2で割った余りならば、x%%2), また、論理比較の演算子で aとbが等しければTRUEをかえす演算子は ==である。
ヒント: 条件つき参照とis.evenをもちいる
数理モデルを構築する方法は２つに大別される。例えば、生命システム（システム生物学、人工生命、生命の起源など）や、社会系などの複雑系の数理モデルを構築する場合に、
の２つが考えられる。
相互作用によるモデル化については、マルチエージェントシステム(Multi Agent System, MAS)がよく用いられてきた。MASとは相互作用の主体、分子や動物など、の特徴や相互作用の仕方をエージェントとして抽象化して、これらエージェント間での相互作用によりモデル化を行う方法である。
このMASをさらに単純化してエージェントを"分子"とみなすこともできるだろう。ただこの分子にはMASのような内部構造はなく、ただ分子間の相互作用のみで記述されるとすると、化学反応をアナロジーとしてもちいるとモデル化を行うことができる。このように、相互作用のモデル化のために化学反応をアナロジーとしてもちいた数理モデルは"抽象化学系 (Abstract Chemistry, AC)"とよばれる。
抽象化学系
反応槽を記号"a", "b"などの集合とする。普通の集合の場合(素朴集合, simple set)には{a,a}は{a}と表記するが、化学系のモデル化においては{a,a,b,b,b}のように要素の重複を許す集合、マルチ集合(multiset)を用いる。そして、反応規則を「マルチ集合の書き換え」としてモデル化する。たとえば、a,b→cはa分子1つとb分子1つからなるマルチ集合{a,c}を{c}に書き換えるとする。また a,a,b→c,cのように"a"が2分子による反応は{a,a,b}から{c,c}へのマルチ集合の書き換えとして記述される。
よって抽象化学系とは、
○反応槽：マルチ集合
○反応規則：マルチ集合の対
としてモデル化される。
例えば
反応槽：{a,b,b}
反応規則：a,b→c, b,c→a
である場合にこの２つの反応規則のうちa,b→cは適用できるが、b,c→aは適用できない（ことはわかるであろう）。この例から反応規則の適用条件を考えてみると、もし、反応規則の左辺が反応槽に真に包含される場合、右辺と書き換える。そうでない場合は、反応規則は適用されない。
まずは１回の反応をモデル化してみよう。そのためには、
を行えばよい。これを記号処理で行ってもよいが、反応槽や反応規則は記号の数により、そのならび(abやba)は関係がない。そこで、これらをベクトルとしてモデル化すると扱いやすくなる。例えば、
反応槽 {a,a,b}は
(a,b,c)=(2a,b,0)=(2,1,0)
とすることができる。
反応規則については、
a,b → c
であれば、
(a,b,c)=(a,b,0)=(1,1,0)
(a,c,c)=(0,0,c)=(0,0,1)
とベクトルの対として表現することができる。
演習 １回の反応についてモデル化せよ
反応規則
ほとんどの場合、反応規則は１つではなく複数である。その場合の扱いについてはさmざまな方法が考えられるが、ここでは"行列"として考えて扱うとしよう。上述したように、反応規則とは"ベクトルの対"である。
例えば、a,b->c,cの規則であれば、"a,b"は(1,1,0),"c,c"は(0,0,2)となるから、この反応規則は＜(1,1,0), (0,0,2)＞となる。これを1,1,0,0,0,2と１つの行ベクトルとあらわすことにする。
演習
① a->b, b,c->a,aをベクトルの対として表現せよ
ヒント：行ベクトルを積み重ねるとは。。
(1,1,0,0,0,2)
(0,0,1... )
(1,0,..... )
から行列
1,1,0,0,0,2
0,0,1,.....
1,0,... ..
を求めることである。この行列を次の演習のために行列Aとする。
演習
①ルールの行列Aから、ルール１の右辺と左辺を求める、スクリプトを示せ。
ヒント：A[行, 列]などの行列やベクトルのフィルタを示せばよい。
② ①の演習でもとめたスクリプトをもちいて、n番目のルールの左辺と右辺を返す関数をつくれ
ヒント:
get.nth.rule<-function(R, th){
（ここに演習①でつくったスクリプト）
}
とすればよい。ここでRとはルールの行列A, thはn番目のルールが対応する。たとえば、以上でつくったAの場合
＞get.nth.rule(A, 1)
とすると
1 1 0 0 0 2
を返すようにせよ。
③ ①、②をもとにマルチ集合msをn番目のルールでかきかえるプログラムを作成せよ。
ヒント：
rewriting.1step<-function(R, th, ms){
lhs<- NULL #書き換え規則の左辺のためのオブジェクト
rhs<- NULL #書き換え規則の右辺のためのオブジェクト
tmp<- NULL #get.nth.ruleの結果代入のためのオブジェクト
tmp<-get.nth.rule(R, th)
lhs<-(ここにルールの左辺を求めるフィルタ)
rhs<-(ここにルールの右辺を求めるフィルタ)
ms<-(ここに演習でつくった書き換えを行うスクリプト)
}
これで、ARMSの基盤となる部分のプログラムは完成したことになる。複数回ルールを適用する場合に問題になるのは「どの規則を適用するのか？」である。だが、それにすすむ前に、相互作用でモデル化すること、について慣れておこう。
相互作用ルールは以下の３つに大別される、
増加：a -> a,a
減少：a,a->a
変化：a -> b
そして、この３つをくみあわせる、たとえば「増加と変化」をくみあわせると、
a-> a,b
のようになる。これはaが変化してbになっている。aは反応の前後で変化しないが、bは反応の後で増加している。
演習：数理モデルをつくる
ARMSから抽象化学系へ
MASのような離散系で複数のルールがある場合に"どのようにルールを適用するか？"は重要なポイントである。ルールの適用の仕方により系の時間発展は大きく変わってしまう。また、ルールを"逐次的"に適用するのか？ルールを"並列に"適用するのか？によっても挙動は大きく変化してしまう。ARMSは反応速度論に基づいて反応規則を定めることにする。
反応速度論とは化学反応による濃度変化の動学であるが、直感的には"反応基質の濃度が高ければ反応しやすい（反応速度が速くなる）"し、"濃度が低ければ"反応がしにくくなる（反応速度が遅くなる）。「なぜ、濃度が高くなると反応がしやすくなるのだろう？」"濃度が高い"とは溶液中にその基質が多く溶け込んでいることを示している。よってそれだけ溶液中で"出会いやすく"なる。だが、それらの基質が"味噌"分子のように水に比べて重い分子はやがて反応槽の底にたまってしまうだろう。反応速度論を考える場合には、溶液はよく撹拌されているとし"味噌分子"のような重い分子が反応槽の底にたまってしまう場合は考えない（だが、もちろんそのような場合を考えることは可能である。それにMASのような局所的な相互作用で系が時間発展していく場合にはエージェントの空間配置も重要になるだろう）。
さて、溶液中で分子同士が衝突する場合、！☆！☆ブバーーーーン！☆！☆と思いっきり衝突する場合もあるだろうし、その逆にちょっとかする程度に衝突する場合も考えられる。ここでの「どの程度、思いっきり衝突するのか？」は反応をかんがえる場合に重要である。なぜなら、化学反応を生じさせるには"そのためのエネルギー"が必要だからである。"そのためのエネルギー"とは化学反応の生じやすさを与えるものである。もし、このエネルギーが小さい場合には"ちょっと強めに衝突すれば"反応が生じてしまう。その一方で、このエネルギーが大きい場合にはちょっと強めに衝突したぐらいでは反応は生じず、かなり"思いっきりぶつからないと"反応が生じない。この"反応を生じさせるためのエネルギー"は反応の速度を左右する（もちろんそれだけで定まるのではなく、基質の濃度も大きく影響するわけだが...）。
では以上を数理モデル化してみよう。まず以下では基質の濃度...いま私たちは離散系で考えているから"分子の数"とみなしてもよいだろう、例えば基質Aの濃度（または分子の数）は[A]のように[,]を用いてあらわす。たとえば A,B-(k2)->B,Bという反応をかんがえてみよう。この反応ではA分子とB分子が衝突することにより反応が生じる、その総当たりの衝突の組み合わせは、
[A] × [B] = 10通り
となる。たとえば[A]=5, [B]=2の場合にはトータル4×2通りの衝突の仕方がある。また、以上でk2とは反応速度定数を表す。この定数はこれらの衝突のうち反応を生じさせる反応の確率を示している。たとえば、k2=0.1ならば10回衝突するうちの1回は反応が生じることをしめす。よってこの場合に反応が生じる回数は期待値として与えられ、それは「全衝突の組み合わせ」×k2 = 10×0.1=1回、となる。
いまはルールが１つしかないので、どのルールを適用するか？を考える必要はないが、以上にたいしてさらに、
A --k1-> A,A
B --k3->
の２つのルールを追加する。
以上よりルールの適用回数の期待値からA,B->A,Aのルールを適用する確率をもとめることができる：
P_{r1} = k2[A][B] / (k1[A]+k2[A][B]+k3[C]) = 1 / (0.5 + 1 +0.2) = 0.59
これと同様にして他のルールの適用の確率も求めることができるであろう。
さて、では以上をR言語で記述していくことにしよう。すでに前回までにARMSの基盤となる部分のプログラムは完成しているので、このプログラムを修正していくことによりプログラムを作成していこう。
＃注意！＃ このように既にあるプログラムを改変したり編集したりして作業していくことがよくあるが「必ず、改変する前のプログラムは別に保存しておくこと！」これを怠ると、改変していってプログラムが全く動かなくなってしまった場合に"どうしょうもなくなる"...なので、変更を加えていく場合にはかならずオリジナルのプログラムは残しておいて、いつでもそのプログラムに戻れるようにしておくことが重要なのである。毎年、卒業研究や修士論文の追い込みの時期になって「プログラムが全く動かなくなってしまいました」といってくるひとがいる。この"天災"は締め切りにちかければちかいほどにおきやすいのである。「天災は忘れたことにやってくる」。
さてさて、いま手元にあるプログラムは「ルール集合のnth番目を指定するとマルチ集合を書き換えてくれる」ものであるはずだ。なので、方針としては、まずnth番目を求めるのにあたって、上述したような仕方で確率を計算してくれるスクリプトを作成することである。つまり、ルールの"番目"書いてある「イカサマさいころ」を作成してこれを"振って"...nth番目を得る。なぜ"いかさま"か？といえば、このさいころの目の出方はルールの適用回数の期待値により、その期待値は基質の濃度によって変化していくからである。
このさいころの計算を行うにはルールの左辺のみがあればよい。A,B->B,Bにおいて右辺のB,Bはさいころの計算を行うためには関係がない。右辺が関係するのは反応が実際に生じる場合であり、とりあえず必要になるのはルールの左辺である。
A,B->B,B
からA分子とB分子の総当たりの衝突の組み合わせ[A]×[B]をもとめた計算であるが、もしこれがA,A,B->となっていたらどうなるであろう？この場合には[A][A][B]=[A]^2[B]となる。このように同種分子同士の衝突はベキ乗で計算を行うことになる（A,Aのような同種分子の衝突について、その衝突の組み合わせは正しくは[A][A-1]である。だがここでは煩雑になるので特に[A][A]としている。分子数が多い場合にはそれほど気にしなくてよいであろう）。
さて、上述したルールをベクトルで表記すると、
＜(1,1),(0,2)＞
となる。また、マルチ集合をms={A,A,A,A,A,B,B}=(5,2)とすると、ここで行っている計算とは、マルチ集合の対応する各々の要素のベキ乗を掛け合わせる、ことになる。
そこで、まずベキ乗のスクリプトを作成する。
pow<-function(x,a){
if (a==0)
1
else
x^a}
演習
① このスクリプトを実行させて、動作を確認せよ。
② mapplyを用いてpowをベクトルに作用させることにより、ルールの左辺を計算することができる。
＞prod(mapply(pow,ms,lhs))
を実行させて動作をたしかめよ。msはマルチ集合, lhsはルールの左辺である。
以上より、ルールを確率的に得るスクリプトを作成すると、例えば、
pow<-function(x,a){
if (a==0)
1
else
x^a}
prob_ls<-function(R,ms,rconst){
nume<-NULL
denom<- 0
len<-dim(R)[1]
if ((dim(R)[2] %% 2) != 0)
print ("Error in the rule list")
else
{wd<-dim(R)[2] / 2
for (i in 1:len){
lhs<-R[i,][c(1:wd)]
h_i <- prod(mapply(pow,ms,lhs)) * rconst[i]
nume <- c(nume,h_i)}
nume / sum(nume)}}
btw<-function(rnum, p_ls){
len<-length(p_ls)+1
i<-1
while(i < len){
if (rnum <= p_ls[i])
{kek<-i
i <- len}
else
i<- i+1}
kek}
roulet<-function(prob_ls){
len<-length(prob_ls)
p_ls<-NULL
tmp<-0
rand<-0
for(i in 1:len){
tmp<-tmp+prob_ls[i]
p_ls<-c(p_ls, tmp)}
rand<-runif(1)
btw(rand, p_ls)
}
となる。
演習:このスクリプトを読んで理解しなさい。そして、
r_1<-c(c(1,0),c(2,0))
r_2<-c(c(1,1),c(0,2))
r_3<-c(c(0,1),c(0,0))
R_1<-rbind(r_1,r_2,r_3)
などとして適当なルールを定義して、
＞roulet(prob_ls (R_1, c(10,10), c(0.1,0.1,0.1))
を実行させて、ルーレットによりルールが確率的に選択されていることを確認しなさい。また、マルチ集合や反応速度定数などを変えて、ルールの選択が変化することを確認しなさい。
以上で、確率的なARMS系を構築することができた。この計算モデルはマルチ集合の要素数が充分に多い場合には、その振る舞いは微分方程式と一致する。
演習：
①確率的ARMSで1ステップの書き換えを行うスクリプトは、例えば、
arms.1step<-function(R,ms,rconst){
th <- roulet(prob_ls(R,ms,rconst))
ms<-rewriting(R,th,ms)
ms}
となる。このスクリプトを用いて複数回書き換えを行うプログラムを作成せよ。
② ARMSで要素数の変動をグラフとして表示するプログラムを作成せよ。
（その１）
A,X---k1->X,X
B,X-k2->Y
X,X,Y-k3->X,X,X
X-k4->
はBelousov-Zhabotinskii反応とよばれる、非線形化学振動系のモデルである。このモデルではk2,k3のパラメータをうまく変えていくと、安定解（収束）が振動生ずる(Hopf分岐)することが知られている。モデル化した確率的ARMSを用いてその振る舞いを調査しなさい。
（その２）
P --k1-> P,P
P,H -k2-> H,H,V
V,H,C-k3->C,C
C --k4->
は香りにより植物が天敵の誘因を行うことで、害虫を間接的に駆除するある化学生態系をモデル化したものである。Pは植物、Hは害虫、Vは香り、Cは天敵をそれぞれ表している。