進化的アルゴリズムで関数を進化させてみる

生命の起源や進化論に関する本を読むと，我々人間は偶然の産物のようである．生まれたばかりの地球はマグマに覆われた火の玉であり，生物が発生するどころか住める環境ではない．クマムシにとってもあったかいどころではすまない星であった．やがて地球が冷めてくると水蒸気が水になって海ができるが，まだその海はただの水たまりであって生物はまだ発生していない．そんなある日になにかしら奇跡が起きてアミノ酸やらRNAやらができてそいつらがまた奇跡を起こして細胞を作り，細胞が遺伝子をごちゃごちゃいじっていたら偶然が偶然で重なっていつの間にか我々一人一人ができていたという感じらしい．

もちろん，生物が進化していったのは単に遺伝子に変異を重ねたからだけではなく，自然淘汰が重要な役割を果たす．遺伝子を改変した結果があまりに無茶苦茶であるとその生物はその命を維持することもできない．また，環境が急変した時には特定の表現型を持つ個体しか生き残ることができない．このように生存，生殖に優れた個体が生き残りやすいように自然淘汰という原理が進化を大いに加速させているのである．

だが，自然淘汰が優れた個体を残すとしても，そのシステムが功を奏するにはそもそも優れた個体がランダムな変異の中から発生していなくてはならない．そのような最適者はいつも生まれてくるのだろうか？その最適者は毎回変わってしまうのだろうか？地球と人間で話を置き換えると，もう一度地球を生まれさせた時，人間のような高度な文明を築き上げる生命体は現れるだろうか？また，そのもうひとつの地球に生まれた文明的生命体は人間なのだろうか？

最近出版された進化の謎を数学で解くという本には，この最適者の到来というダーウィン以来の未解決問題への一つの解が示されている．著者のアンドレアス・ワグナーのグループは，生物における生化学反応を5000個リストアップし，その反応の有無の組み合わせを表現型として持つ個の点からなる空間を考えた．その中で生存できる反応の組み合わせを持つ点を1つ選び，そこからランダムウォークで表現型空間を探ることにより，ある点からランダムに変異していってかつ生存できるという近傍をシミュレーションしたのだ．その結果として得られた近傍は表現型空間の中で広範囲に広がっており，様々な機能を持つ点と隣接していたらしい．つまり，様々な環境に対する最適者が誕生する余地があるということである．

この本を読んでいて思い出したのが，関数をランダムに変化させて進化させてみたいという考えである．プログラミングを学び始めた時から自分で自分を改良していくプログラムがあったら面白いなと思っていたのだが，その改良がランダムな改変であるプログラムというわけである．自分で自分を書き換えるプログラムというのはいわゆる自己書き換えコードというやつらしく，また進化的操作（交叉，突然変異，自然淘汰など）を組み込んだプログラミング手法は進化的プログラミングとしてかなり研究，応用も進んでいるようだ．なので，別に発想としては全く新しいことではないが，自分で思ったことをプログラムとして書いて動かすというのは最先端の研究でなくても楽しいものである．

というわけで，進化的アルゴリズムで関数を目的関数に進化させてみるというプログラムをRで書いてみた．Rを使ったのは私がある程度使い慣れているという他に，eval関数が用意されているというのが理由である．eval関数は文字列を式として評価する関数であり，例えば以下のコードを実行するとsin(x)が計算されて1という結果が得られる：

x <- pi/2
str <- "sin(x)"
eval(parse(text=str))

上のコードではxにすでに値が代入されているが，xを変数として扱って関数にすることもできる．次のコードでは，myfuncとして文字列strから作った関数を宣言している．そのため，myfunc(pi/2)は1を返し，またcurve関数にmyfunc(x)を入力するとサインカーブがプロットされる：

str <- "sin(x)"
myfunc <- function(x) eval(parse(text=str))
myfunc(pi/2)
curve(myfunc(x))

これを利用すれば，ランダムに生成した関数を目的の関数に進化させることができる．まず，関数に使う文字をリストアップする：

moji <- strsplit("xxxxxxxxxxxxxxxxxxxx+-*/^+-*/^+-*/^+-*/^<> ()%;<> ()%;<> ()%;1234567890","")[[1]]

これで変数mojiにはこれらの文字が１文字ずつのベクトルとなって代入される．数字が少なく，やたらxが多かったりするのは使用頻度の調節のためである．ここは個人のお好みになるだろう．

このmojiからランダムにサンプリングしてきて適当な式をざっと1000個ほど作る．式を作るというと複雑に聞こえるが，実際にはランダムな文字列を生成しているだけである．つくった文字列リストfuncsの最初の10個を見ると，やはり無茶苦茶な式ばかりである：

> funcs[1:10]
 [1] " ^;;^x x<x12^( *+x+x/ %x^x4( x*x2;+>)%2; %(-x*;(("                                         
 [2] " x+x;x; x;*<+*x<x4/+(<x<xx;x>4x"                                                           
 [3] "<x)x *(x^^x>x%x%x**x xx%<(^^^x%<^;/;1;*>*//^-x+x^x^x%-*%*xx(x--xx-x^%xxx^-x^x"             
 [4] "x)- <+2xxx*x+;x< ^x^^/xx2(%xx^>x1x/%>/^%<%(+x^ >xxxxx;x%;xx xxx>-xx(x%)^)xx3^xx<x%+;()/)*x"
 [5] "x;>-;<x4<x)*x<<x+-x^x )x*/x^x*-x --2^ x*^;(x+-;- x++x+"                                    
 [6] "(; >;(/x<*/*x(x*^)-%4 *xx*>x*xx/ )<*xx%x- -2x(%x/x/*x;xx*x3)<x-/-*;xx)-<x);-"              
 [7] "xx+*x-(;"                                                                                  
 [8] "-(x)x2x/x%%x>)+;x2*<x++x)^) x--x*<xxxx1+(x^ ^/<-^x3x))x/x-x<x4((2x/)/2>;x*-x+;/< x"        
 [9] "xxx 1(xxx);+(/xxx-+/x>4/*;<x1-/+*+^x*/x-/*<^x+x;)x*(^%)><x>)-+x+ )4+;^;)x*3x>()x>x+++/"    
[10] "4/^+%x-x--xx*^/^x;xx;%%)^<2)x-*-x)/<xx//x"       

ここから構文エラーのないもののみを抜き出す．

suppressWarnings(te <- try(testf(target.kizami)))
te <- try(testf(target.kizami))

testfをfuncsの式から作った関数とし，そこに引数target.kizamiを代入してみる．ここで，tryをかましておいてエラー判定を行う．これを1000個もやるとエラーや警告が大量に表示されることが予想されるので，suppressWarningsやoptions(show.error.messages=FALSE)を使うことをおすすめする．1000個作って意味がある文は高々10個くらいである：

> length(funcs)
[1] 7
> funcs
[1] "--1" "3"   "-3"  "x "  "x"   "x"   "x*4"

あとはこの文字列をいじりつつ何かしら順序をつけて淘汰するというサイクルを繰り返せば良い．今回は関数を目標関数として，との]のノルム（つまり作った関数をとしてを計算する）を評価関数に設定した．構文エラーのない文字列は，の距離で順位をつけて半分だけを生き残らせ，生き残った文字列に文字列の挿入，削除，変異をランダムに加えて新たに1000個の文字列を作る．こうして1000世代作るシミュレーションを12回行った．

結果はなんと，12回とも同じ関数が出現した．その証拠に，12回のシミュレーション後に残った世代のうち，距離の上位10位を拾ったデータが次である（番号付けが5から始まっているのは，1~4番は条件を変えているため）：

> load("smL1-5.RData")
> Fx[order(Fx$L1dist)[1:10],]
         expr len   L1dist
1   8*x*-x+20   9 95.13315
4   8*x*-x+20   9 95.13315
5  8*-x* x+20  10 95.13315
7   8*x*-x+20   9 95.13315
8   8*x*-x+20   9 95.13315
9  8*x*-x +20  10 95.13315
10  8*x*-x+20   9 95.13315
12  8*x*-x+20   9 95.13315
14 8*x*-x++20  10 95.13315
16  8*x*-x+20   9 95.13315
> load("smL1-6.RData")
> Fx[order(Fx$L1dist)[1:10],]
        expr len   L1dist
3   20-8*x*x   8 95.13315
4   20-8*x*x   8 95.13315
9   20-8*x*x   8 95.13315
11  20-8*x*x   8 95.13315
14  20-8*x*x   8 95.13315
15  20-8*x*x   8 95.13315
16 20-8* x*x   9 95.13315
19 20-8 *x*x   9 95.13315
20 +20-8*x*x   9 95.13315
21  20-8*x*x   8 95.13315
> load("smL1-7.RData")
> Fx[order(Fx$L1dist)[1:10],]
         expr len   L1dist
4    20-x*8*x   8 95.13315
5   20-x*8* x  10 95.13315
6    20-x*8*x   8 95.13315
7    20-x*8*x   8 95.13315
8    20-x*8*x   8 95.13315
10  20-x* 8*x   9 95.13315
11   20-x*8*x   8 95.13315
13  20-x*+8*x   9 95.13315
14 20 -x* 8*x  10 95.13315
15   20-x*8*x   8 95.13315
> load("smL1-8.RData")
> Fx[order(Fx$L1dist)[1:10],]
          expr len   L1dist
6  20+-+x*8* x  11 95.13315
7     20-x*8*x   8 95.13315
8   20-x *8*+x  10 95.13315
12    20-x*8*x   8 95.13315
13   +20-x*8*x   9 95.13315
14    20-x*8*x   8 95.13315
15    20-x*8*x   8 95.13315
23    20-x*8*x   8 95.13315
24    20-x*8*x   8 95.13315
26    20-x*8*x   8 95.13315
> load("smL1-9.RData")
> Fx[order(Fx$L1dist)[1:10],]
        expr len   L1dist
1   20-x*x*8   8 95.13315
2  20 -x*x*8   9 95.13315
3   20-x*x*8   8 95.13315
5  20-+x*x*8   9 95.13315
6   20-x*x*8   8 95.13315
7   20-x*x*8   8 95.13315
12  20-x*x*8   8 95.13315
14 20-+8*x*x   9 95.13315
16  20-x*x*8   8 95.13315
19  20-x*x*8   8 95.13315
> load("smL1-a.RData")
> Fx[order(Fx$L1dist)[1:10],]
        expr len   L1dist
1  20-8*x* x   9 95.13315
2   20-8*x*x   8 95.13315
3   20-8*x*x   8 95.13315
4  20- 8*x*x   9 95.13315
5   20-8*x*x   8 95.13315
7   20-8*x*x   8 95.13315
9   20-8*x*x   8 95.13315
10  20-8*x*x   8 95.13315
11  20-8*x*x   8 95.13315
13  20-8*x*x   8 95.13315
> load("smL1-b.RData")
> Fx[order(Fx$L1dist)[1:10],]
          expr len   L1dist
2     20-8*x*x   8 95.13315
3     20-8*x*x   8 95.13315
5  20 -8* +x*x  11 95.13315
6     20-8*x*x   8 95.13315
7    20+-8*x*x   9 95.13315
8     20-8*x*x   8 95.13315
10    20-8*x*x   8 95.13315
11    20-8*x*x   8 95.13315
13   20-8*x*x    9 95.13315
14    20-8*x*x   8 95.13315
> load("smL1-c.RData")
> Fx[order(Fx$L1dist)[1:10],]
         expr len   L1dist
2    20-x*8*x   8 95.13315
3    20-x*8*x   8 95.13315
4    20-x*8*x   8 95.13315
5  20-+x*8 *x  10 95.13315
7    20-x*8*x   9 95.13315
8   20 -x*8*x   9 95.13315
10  20-x*+8*x   9 95.13315
12  20-x*8*x    9 95.13315
14  20-x*8*+x   9 95.13315
17  20-x*8 *x   9 95.13315
> load("smL1-d.RData")
> Fx[order(Fx$L1dist)[1:10],]
         expr len   L1dist
3    20-x*x*8   8 95.13315
5   20-x*x*8    9 95.13315
6    20-x*x*8   8 95.13315
7   20-x*x* 8   9 95.13315
9  20+-x* x*8  10 95.13315
10  20-x*x *8   9 95.13315
12   20-x*x*8   8 95.13315
14   20-x*x*8   8 95.13315
16   20-x*x*8   8 95.13315
17   20-x*x*8   8 95.13315
> load("smL1-e.RData")
> Fx[order(Fx$L1dist)[1:10],]
        expr len   L1dist
1   20-x*8*x   8 95.13315
2   20-x*8*x   8 95.13315
3   20-x*8*x   8 95.13315
4   20-x*8*x   8 95.13315
5   20-x*8*x   8 95.13315
8  20-x*+8*x   9 95.13315
9   20-x*8*x   8 95.13315
10  20-x*8*x   8 95.13315
11 20-x*+8*x   9 95.13315
12  20-x*8*x   8 95.13315
> load("smL1-f.RData")
> Fx[order(Fx$L1dist)[1:10],]
        expr len   L1dist
2   20-x*8*x   8 95.13315
3  20- x*8*x   9 95.13315
5   20-x*8*x   8 95.13315
6   20-x*8*x   8 95.13315
7   20-x*8*x   8 95.13315
8  20-x*8*+x   9 95.13315
10  20-x*8*x   8 95.13315
12 20-x*8*+x   9 95.13315
14  20-x*8*x   8 95.13315
18 20+-x*8*x   9 95.13315
> load("smL1-0.RData")
> Fx[order(Fx$L1dist)[1:10],]
         expr len   L1dist
1   x*-8*x+20   9 95.13315
2   x*-8*x+20   9 95.13315
3   x*-8*x+20   9 95.13315
4   x*-8*x+20   9 95.13315
8   x*-8*x+20   9 95.13315
10  x*-8*x+20   9 95.13315
13 x*-8*x+20   10 95.13315
15  x*-8*x+20   9 95.13315
18  x*-8*x+20   9 95.13315
19  x*-8*x+20   9 95.13315

目標関数であったと，を]でプロットすると，ぴったりではないがまあいい線を行っていることがわかる：

[8,10]は諦めたのか？

また，世代ごとの最小距離をプロットすると，12回ともだいたい100世代行かないくらいで急激に進化し，その後進化は落ち着くという具合になっている：

以上を踏まえて最初の疑問「最適者はいつも生まれてくるのだろうか？」そして「その最適者は毎回変わってしまうのだろうか？」を考えてみる．まず，「最適者はいつも生まれてくるのだろうか？」から考察しよう．「最適者」が何を意味するかは考えどころであるが，今回現れたは目標関数からまあそんなに外れてはいない．だが，「最適」なのはやはり目標関数そのものである．その意味で言うと今回の計算実験では「最適者」は現れなかった．次に，「その最適者は毎回変わってしまうのだろうか？」であるが，「最適者」をとするならば明らかにこの疑問に対する答えは「毎回同じ」である．

これを地球と人間に置き換えてみると，こうなる：「地球は何回生まれてもやはり人間が現れるが，人間は最適な存在ではない．」まあ，あんまりわかったようなことを言うとそれを専門にしている方々から怒られそうなのでそろそろ記事を終わることにしよう．

最後におまけで，今回実行したソースコードを貼り付けておく．趣味で作ったので作りはだいぶ雑でコメントもない．

options(show.error.messages=FALSE)
target.f <- function(x) x*(x-1)*(x-3)*(x-10)
target.intv <- c(0,10)
target.N <- 10000
target.kizami <- seq(target.intv[1],target.intv[2],length=target.N)
target.val <- target.f(target.kizami)

moji <- strsplit("xxxxxxxxxxxxxxxxxxxx+-*/^+-*/^+-*/^+-*/^<> ()%;<> ()%;<> ()%;1234567890","")[[1]]
moji.l <- length(moji)

N <- 1000
lmax.init <- 50
funcs <- as.character(numeric(N))
funcs.limit <- 100
funcs.L1dist <- numeric(N)
l.init <- sample(1:lmax.init,N,TRUE)

for(i in 1:N){
	funcs[i] <- paste(moji[sample(1:moji.l,l.init[i],TRUE)], sep="", collapse="")
	testf <- function(x) { eval(parse(text=funcs[i])) }
	suppressWarnings(te <- try(testf(target.kizami)))
	if(class(te) == "numeric"){
		funcs.L1dist[i] <- sum(abs(te-target.val))/target.N
	}else{
		funcs.L1dist[i] <- Inf
	}
}

F0 <- data.frame(expr=funcs[funcs.L1dist<Inf], len=nchar(funcs[funcs.L1dist<Inf]), L1dist=funcs.L1dist[funcs.L1dist<Inf])

F0 <- F0[!is.na(F0$L1dist),]

Fx <- F0

nextF <- function(beforeF){
F0 <- beforeF

F0.ord <- order(F0$L1dist)
oya.n <- round(length(F0$expr)/2)
funcs.oya <- as.character(F0$expr[F0.ord[1:oya.n]])

for(i in 1:N){
	hoge <- strsplit(funcs.oya[(i%%oya.n)+1],"")[[1]]
	funcs[i] <- paste("d", paste(hoge,rep("d",length(hoge)),sep=""),sep="",collapse="")
	for(j in 1:nchar(funcs[i])){
		if(runif(1)<0.05){
			substr(funcs[i],j,j) <- moji[sample(1:moji.l,1)]
		}else if(runif(1)<0.05){
			substr(funcs[i],j,j) <- "d"
		}
	}
	funcs[i] <- gsub("d","",funcs[i])
	testf <- function(x) { eval(parse(text=funcs[i])) }
	suppressWarnings(te <- try(testf(target.kizami)))
	if(class(te) == "numeric"){
		funcs.L1dist[i] <- sum(abs(te-target.val))/target.N
	}else{
		funcs.L1dist[i] <- Inf
	}
	if(nchar(funcs[i])>funcs.limit){
		funcs.L1dist[i] <- Inf
	}
}

F1 <- data.frame(expr=funcs[funcs.L1dist<Inf], len=nchar(funcs[funcs.L1dist<Inf]), L1dist=funcs.L1dist[funcs.L1dist<Inf])
return(F1[!is.na(F1$L1dist),])
}

Ngen <- 10000
dist.gen <- numeric(Ngen)

dist.gen[1] <- min(Fx$L1dist)

for(i in 2:Ngen){
	Fx <- nextF(Fx)
	dist.gen[i] <- min(Fx$L1dist)
}
save.image("smL1-5.RData")