![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
RでGA(遺伝的アルゴリズム)
(工事中)
今のところ Ca' Foscali University(Venezia)のga package、Sydney工科大学のgafit packageおよび rgenoud そして、seaoが存在する。
library(ga)
Decode (binary文字列(gene)を数値にする)
binary2integer binary2real
Encode (数値を binary文字列(gene)にする)
integer2binary real2binary
扱う対象の数値範囲などを考えながら binary文字列長あるいは Decode, Encode方式を考える。
例えば、関数に与える数値が整数でもよいなら integer系の encode/decodeでよいが、 実数にしたいなら 変数の区間(from, to)や binary長(size)を定義した上で real系の encode/decodeを 使う必要がある。
関数の最小/大値を探査するためには genetic関数を使って geneticオブジェクトを 生成しないといけない。以降、このオブジェクトを使って、 新世代の生成などを実施することになる。
例えば binary文字列長が5でよいのであれば
exa <- genetic(encode=integer2binary, par.encode=list(size=5), decode=binary2integer, length.string=5, fitness=fit)
とする。
対象とする関数は
fit <- function(x) x^2
などとして*1定義する。
> exa Alphabet 0 1 Wildcard * Length String 5 Size Generation 100 Size Population 20 Encode integer2binary Parameter Encode $size [1] 5 Decode binary2integer Fitness fit Type of Selection roulette Prob. Crossover 0.8 Prob. Mutation 0.01
交配確率(Prob. Crossover = 0.8)や突然変異確率(Prob. Mutation = 0.01)、 一世代中の個体数(Size Population = 20)や世代数(Size Generation = 100) などは指定しないとDefault値を充てられるので、気に入らなければ適宜 指定すること。
ga packageは学習用だけあって、内部でどのようなことをなしているのかを こちらに見えるようにすることも可能である。
が、今はとりあえず手っ取り早く解いてみて、このパッケージの適用性を 示してみる*2。
garesult <- ga.simple(exa) plot(garesult)
結果は plot.ga関数を使用して (上では plot(garesult)となっているが やっていることは同じこと)プロットできる。
このパッケージに附属している関数を使えば 遺伝的アルゴリズム内部で何をやっているのかを 見ることができる。いくつかを紹介しておく。
以下は 値21の示すbinary文字列と 値13の示すそれを交差させたものである。 交差する位置ははじめは1-2の位置となっているが 2回目の試行では3-4となっている。この位置はランダムに決定される。
> x <- encode(21,exa) > x 1~0~1~0~1 > y <- encode(13,exa) > y 1~0~1~1~0 > crossover(x,y,exa) Parents: 1|0~1~0~1 1|0~1~1~0 Offspring: 1|0~1~1~0 1|0~1~0~1 > crossover(x,y,exa) Parents: 1~0~1|0~1 1~0~1|1~0 Offspring: 1~0~1|1~0 1~0~1|0~1
ちなみに、複数の位置で交差させるのは現時点では無理かと思われる。
例えば一世代 5個体とする。世代中の個体の組 x は mutation関数により変異させられる。
> x 1 1 0~0~0~0~1 2 2 1~1~1~1~0 3 3 0~0~1~0~1 4 4 0~1~1~1~1 5 5 0~1~1~1~0 > mutation(x,exa) String before mutation: 0~1~0~0~0~0~1~0~1~1~0~1~1~1~1~0~1~0~1~1~1~0~1~1~0 String after mutation: 0~0~0~0~1~0~1~0~0~1~0~0~1~1~1~0~1~0~1~1~1~0~1~1~0
xをある世代の個体集合、fitを対象とする関数、exaをgeneticオブジェクトとして 与えると、世代の中から比較的優位な個体二つをランダムに選択して親が決定される。
selection(x,fit,exa) $x 0~1~1~1~0 $y 1~0~0~0~0 $parents.position [1] 17 8
彼らの子どもは前述の crossoverを使って作られる。
まずは AdamとEve(たち)がいないとはじまらないので generate.population関数を使って 第一世代を作る。
gen00 <- generate.population(exa)
以降はこの第一世代を用いて generation関数を使って 新たな世代を作っていくことが可能となる。
gen01 <- generation(gen00, exa)
世代に含まれる個体を評価するのはfitness.population関数を使う
fitness.population(gen01, exa)
世代の平均値を知りたい場合は以下を用いる。
fitness.mean(gen01, exa)
なお、先に示した ga.simpleの結果のプロットのうち左図は世代の fitnessのうち最大値を実線、平均を点線で示している。
系譜を取得することが可能である。 優良な個体の血統を知ることができる。 例えば、先ほど実施した simple.gaの結果のうち、3世代から(to) 5世代まで(from) の系譜を調べたいのなら以下のようにする。
この時 fromの方が toより大きくないといけないことに注意(遡る感じで)
> y <- genealogy(,offspring=1,from=5,to=3,exa) > plot(y)
この場合の plotには、plot.genealogyが使われている(が意識する必要はない)
探査対象関数の最小値を到達可能な範囲で捜査する。 探査対象の関数(target)および、 探索開始点(start)、thermal(突然変異発生確率)、maxiter(世代数上限)、 samples(一世代辺りの個体数)、step(step size?)を与える。
こちらは ga packageとは違って収束過程、個体の性質などは 内部で処理しており R上で見えるわけではない。
library (gafit)
まずは関数の定義から。 まず、リファレンスにあるように cos(θ) + sin(θ)から。
e <- expression( cos( theta ) + sin( theta )) gafit( e, list( theta=3 ), step=0.1 ) $theta [1] 3.922567 attr(,"score") [1] -1.414200
この場合 thetaが 3.926991のとき -1.414214で最小値を取る。
x <- seq(3,5,by=0.001) f <- function(x) (cos(x) + sin(x)) plot(x, f(x), type="l") s <- sort(f(x),index=T) x[s$ix[1]]
などで検証可能である。
なお、パラメータが二つになっても適用可能 (ga packageは基本的に一つ)
ここでは有名なRosenbrockのbanana関数に適用することにする。 この関数は (x, y) = (1, 1)で最小値 0となる。 最小値近傍で傾斜が複雑になっているため、最適化ルーチンのベンチマークに 使われることがよくある。
e_banana <- expression( 100 * ( Y - X^2 )^2 + ( 1 - X )^2 ) gafit(e_banana, list( X=0, Y=0 ), step=0.1, maxiter=100) $X [1] 1.000984 $Y [1] 1.002066 attr(,"score") [1] 1.897574e-06
この例*3では良い結果を出すことができた。
他に 複素数も使うことができる。これは ?gafitしてExampleを参照下さい。
最適化
library(rgenoud)
最適化
library(seao)
まずは遺伝子の構成を決める。seaoの特徴的なところは 遺伝子の構成が上記のライブラリよりも自由に決定できることだと思う。
まずはgenomestruc (Genome Structureの略?)を使って 対話的に 遺伝子の構成を決めていく。
今回は 遺伝子の持つ意味を3種類とした。
> gs <- genomestruc(n.genes=3)
---- GENE 1 Name of the gene x Give the minimum value -10 Give the maximum value 10 Give the step (precision) 0.1 Maximum set to: 10 There are 201 alleles Optional: Give a list of names for the different alleles (seperated by a ';') WARNING: Number of allelenames not equal to number of alleles; namevector set to NULL
---- GENE 2 Name of the gene y Give the minimum value -10 Give the maximum value 10 Give the step (precision) 0.05 Maximum set to: 10 There are 401 alleles Optional: Give a list of names for the different alleles (seperated by a ';') WARNING: Number of allelenames not equal to number of alleles; namevector set to NULL ---- GENE 3 Name of the gene flag Give the minimum value 1 Give the maximum value 2 Give the step (precision) 1 Maximum set to: 2 There are 2 alleles Optional: Give a list of names for the different alleles (seperated by a ';') f1;f2
上記で定義した遺伝子構成を表示するには以下。
> str(gs) List of 1 $ genes:List of 3 ..$ :List of 5 .. ..$ name : chr "x" .. ..$ min : num -10 .. ..$ max : num 10 .. ..$ n.alleles: num 201 .. ..$ names : NULL ..$ :List of 5 .. ..$ name : chr "y" .. ..$ min : num -10 .. ..$ max : num 10 .. ..$ n.alleles: num 401 .. ..$ names : NULL ..$ :List of 5 .. ..$ name : chr "flag" .. ..$ min : num 1 .. ..$ max : num 2 .. ..$ n.alleles: num 2 .. ..$ names : chr [1:2] "f1" "f2"
Gene1 (x) は -10から10まで。分割は0.1毎なので 201種。
Gene2 (y) は -10から10まで。分割は0.05毎なので401種。
Gene3 (flag) は f1 (=1)もしくは f2 (=2)である。
この構成をもとに 10個体を生成する
> g00 <- newgen(gs, n.ind=10,method="random")
> g00
# 中略
$generations[[1]]$allele
[,1] [,2] [,3]
[1,] 8.5 -9.25 2
[2,] -8.4 9.00 1
[3,] 9.2 3.55 1
[4,] -4.9 -6.30 2
[5,] 3.3 -6.25 1
[6,] -0.5 -5.45 1
[7,] -6.7 1.05 2
[8,] -2.1 5.10 2
[9,] 1.8 -4.75 1
[10,] -8.8 7.25 2
$generations[[1]]$parents
[1] "-1" "random"
[,1]が x, [,2]がyそして、[,3]がflagになっていることがわかる。
なお、このseaoパッケージだが最適化の目的関数は自前で定義して generation (g00)に与えないといけない。
例えば
> g00$generations[[1]]$fit <- c(1,2,3,4,5,6,7,8,9,10)
以下、 g00(fit定義済)から発生させた g01と、その後 g01のfitを定義せずに g02を発生させた場合に示す挙動を書いておく。 fitがないと収束もへったくれもないので"怒る"ようである。
> g01 <- nextgen(g00,n.ind=10,gen.parent=NULL,selection=list(base="fit", rescale=0),
corate=90,mutation=list(base="unif", spread=1, rate=15))
> g02 <- nextgen(g01,n.ind=10,gen.parent=NULL,selection=list(base="fit", rescale=0),
corate=90,mutation=list(base="unif", spread=1, rate=15))
Error in if (d.maxfit - d.minfit < 1e-12) { :
missing value where TRUE/FALSE needed
>
例えば以下のような目的関数を与えて100世代の計算をさせてみた。
ちなみに、g00は Generation, fitが目的関数、
g00 <- newgen(gs, n.ind=10,method="random");
fit <- function(x) ( 100 * (x[2]-x[1]^2)^2 + (1-x[1])^2 ) * x[3] * -1 + 1000;
tmpfit = c(1:10); # 初期化のつもり
for(j in 1:100) { # 100世代繰り返す
for(i in 1:10){ # 10は個体数
tmpfit[i]=fit(g00$generations[[j]]$allele[i,]);
}
g00$generations[[j]]$fit <- tmpfit ;
g00 <- nextgen(g00,n.ind=10,gen.parent=NULL,selection=list(base="fit", rescale=0),
corate=100,mutation=list(base="unif", spread=1, rate=15));
}
こうして 100世代まで繰り返すと
> g00$generations[[2]]
$fit
[1] 984.724659 906.045219 984.724659 984.724659 984.724659 984.724659
[7] 812.090438 984.724659 842.146659 2.241219
$allele
[,1] [,2] [,3]
[1,] -0.89 0.45 1
[2,] -0.71 -0.45 1
[3,] -0.89 0.45 1
[4,] -0.89 0.45 1
[5,] -0.89 0.45 1
[6,] -0.89 0.45 1
[7,] -0.71 -0.45 2
[8,] -0.89 0.45 1
[9,] -0.89 -0.45 1
[10,] -0.71 -2.65 1
...
$generations[[100]]$fit
[1] -10395.395 999.830 999.830 999.830 661.280 999.830
[7] -5146.720 -9894.494 -3108.970 -437799.376
$generations[[100]]$allele
[,1] [,2] [,3]
[1,] 3.32 0.35 1
[2,] 0.60 0.35 1
[3,] 0.60 0.35 1
[4,] 0.60 0.35 1
[5,] 0.60 2.20 1
[6,] 0.60 0.35 1
[7,] 0.60 8.20 1
[8,] 2.78 0.35 2
[9,] 0.60 -6.05 1
[10,] -8.16 0.35 1
なんとなく最大値 1000に近い値を算出できていることがわかる。
library(subselect)
地図のラベル配置にＧＡを利用
