Machine Learning: An Algorithmic Perspective 12章
遺伝アルゴリズム
Machine Learning: An Algorithmic Perspective (Chapman & Hall/Crc Machine Learning & Patrtern Recognition)Stephen Marsland
Chapman & Hall / ￥ 6,529 ()
通常2～3週間以内に発送
Four Peaks Problemってのがあるらしい。目的関数が
- 最初のビットからの連続した1の長さか最後のビットから連続した０の長さの大きいほう
- 但し、両端の1,0の連続した数が10以上の場合は100ポイント獲得する。
目的関数はこんな感じ。
#!/usr/bin/env python
# -*- encoding:utf-8 -*-
from itertools import takewhile
def o(bits):
return len(list(takewhile(lambda x: x == 1,bits)))
def z(bits):
return len(list(takewhile(lambda x: x == 0,reversed(bits))))
def f(bits):
reward = 100 if o(bits) > 10 and z(bits) > 10 else 0
return max(o(bits),z(bits)) + reward
if __name__ == "__main__":
bits1 = [1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
bits2 = [1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
for bits in [bits1,bits2]:
print "score: %d %s" % (f(bits),bits)
#score: 20 [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
#score: 114 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
要するに素直に探索していくとローカルミニマムに落ちるようになっていて、ピークの数が4つあるのでFour Peaks Problem
連続した1,0の長さでcontour plotを描いた。
#!/usr/bin/env python
# -*- encoding:utf-8 -*-
def sim_score():
for x in range(100):
for y in range(100):
if x + y > 100:
yield 0
else:
score = four_peak(x,y)
yield score
def four_peak(x,y):
reward = 100 if (x > 10 and y > 10) else 0
score = max(x,y) + reward
return score
if __name__ == "__main__":
from pylab import *
delta = 1
x = arange(0, 100, delta)
y = arange(0, 100, delta)
X, Y = meshgrid(x, y)
Z = array([z for z in sim_score()])
Z.shape = 100,100
im = imshow(Z,origin='lower' ,alpha=.9)
colorbar(im)
cset = contour(X,Y,Z)
clabel(cset,inline=1,fmt='%1.1f',fontsize=10)
hot()
savefig('4peaks.png')
LeadOptimizationをシミュレーションしたい
ふと、つぶやいた。
例えば
- プロジェクトで500検体合成
- ひとつの親あたり5-25化合物くらいの子を合成する
- プロジェクトの人員の制限により、一度に10ラインまでしか同時に走らない
- 11番目以降はそれ以降のオプティマイゼーションは行われない
その時知りたいことは、
- {計画、合成、アッセイ、解析}というサイクルは何回まわるのか、
- 有望そうな化合物だけど、人的リソースの関係で選択されなかった歴史はどこにあるか？
- バックアップ候補の化合物はどこら辺で現れるか（初期？中期？）
あたり。
現実の系に近づけるために、親化合物にはポテンシャルがあって、近傍探索すると、ポテンシャルの近くで活性が変化するモデル。あと同じライン（ブランチ）を継続して合成すると、ポテンシャル曲線が底に近づくのと、合成もSARに関する知識がついてくるので、活性は上がりやすくなって、かつ変動幅も小さくなるようにweightはだんだん小さくなるようにしてみた。
#!/usr/bin/env python
# -*- encoding:utf-8 -*-
# kzfm <kerolinq@gmail.com>
from random import random
import networkx as nx
import matplotlib.pyplot as plt
class Branch(object):
"""LeadOptimization flow"""
count = 0
@classmethod
def inc_count(cls):
cls.count = cls.count + 1
return cls.count
@classmethod
def get_count(cls): return cls.count
def __cmp__(self, other):
return cmp(self.act, other.act)
def __init__(self,potency,weight):
self.num = Branch.inc_count()
self.potency = potency
self.weight = weight
self.act = self.potency + self.weight * random()
def make_child(self,num_childs,potency,weight):
return [Branch(self.potency + self.weight * random(),self.weight * 0.7) for i in range(num_childs)]
if __name__ == "__main__":
max_comps = 500
initial_branches = 3
nodesize = []
nrange = []
erange = []
generation = 1
heads = [Branch(0.3,1) for i in range(initial_branches)]
G=nx.Graph()
for h in heads:
nodesize.append(h.act * 30)
nrange.append(generation)
while Branch.get_count() < max_comps:
new_heads = []
generation += 1
for branch in heads[:10]:
for new_comp in branch.make_child(int(5+20*random()),branch.potency,branch.weight):
nodesize.append(new_comp.act * 30)
nrange.append(generation)
erange.append(generation)
G.add_edge(branch.num,new_comp.num)
if new_comp.act > 0.75:
new_heads.append(new_comp)
heads = new_heads
heads.sort()
nx.draw_spring(G, node_size=nodesize, node_color=nrange, edge_color=erange,alpha=0.7, \
cmap=plt.cm.hot, edge_cmap=plt.cm.hot, with_labels=False)
plt.savefig("proj.png")
んーなんかいまいち。一応活性の強弱にあわせてノードのサイズが変わるようにしてみたんだけど、大きさがわからん。
matplotlibの本
Pythonの描画ライブラリであるmatplotlibの本が出てた。
Matplotlib for Python Developers欲しいのう
Box-MullerをPythonで
PRML11章。一様乱数から正規乱数を発生させる。
from math import sin,cos,pi,log,sqrt
from random import random
def box_muller():
z1 = random()
z2 = random()
return sqrt(-2*log(z2)) * sin(2*pi*z1)
if __name__ == "__main__":
from pylab import *
x = [box_muller() for i in range(10000)]
hist(x,20)
savefig("box_muller.png")
パターン認識と機械学習 下 - ベイズ理論による統計的予測macbookのmatplotlibでLookupError
macbookでmatplotlibを使っていたら出たエラー。
LookupError: unknown encoding: X-MAC-JAPANESE
/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/matplotlib/cbook.pyというファイルを書き換えた
def unicode_safe(s):
if preferredencoding is None: return unicode(s)
elif preferredencoding == 'X-MAC-JAPANESE': return unicode(s)
else: return unicode(s, preferredencoding)
matplotlibでgumowski
ローレンツのついでにgumowski
class Gumowski:
"Gumowski-Mira Map"
def __init__(self, x0, y0, a):
self.x0 = x0; self.y0 = y0
self.a = a;
def f(self,x):
return self.a*x + 2*(1-self.a)*x**2/(1+x**2)
def calc(self, n):
dat = []
x,y = self.x0, self.y0
for i in range(n):
nx = y + self.f(x)
ny = -x + self.f(nx)
dat.append([nx,ny])
x,y = nx,ny
return dat
if __name__ == "__main__":
from pylab import *
from random import random
n = 3
for i in range(n*n):
p = str(n) + str(n) + str(i+1)
subplot(p)
a = random() * 0.1 + 0.3
gumow = Gumowski(20, 20, a)
dat = gumow.calc(1000)
x = [d[0] for d in dat]
y = [d[1] for d in dat]
plot(x,y,'bo', alpha=0.1,ms=5)
savefig("gumow.png")
matplotlibでローレンツアトラクタ
どう書くより。
python+matplotlibで。
class Lorenz:
"Lorenz attractor"
def __init__(self, x0, y0, z0, p, r, b, dt):
self.x0 = x0; self.y0 = y0; self.z0 = z0
self.p = p; self.r = r; self.b = b
self.dt = dt
def calc(self, n):
dat = []
x,y,z = self.x0, self.y0, self.z0
for i in range(n):
dx = (-1 * self.p * x + self.p * y) * self.dt
dy = (-x * z + self.r * x - y) * self.dt
dz = (x * y - self.b * z) * self.dt
x += dx; y += dy; z += dz
dat.append([x,y,z])
return dat
if __name__ == "__main__":
from pylab import *
lorenz = Lorenz(1.0, 1.0, 1.0, 10.0, 28.0, 8.0/3.0, 0.01)
dat = lorenz.calc(5000)
x = [d[0] for d in dat]
y = [d[1] for d in dat]
plot(x,y)
savefig("lorenz.png")
Beginning Python Visualization
matplotlibがメインなのかな？目次だけだとよくわからん。
Beginning Python Visualization: Crafting Visual Transformation Scripts (Books for Professionals by Professionals)Shai Vaingast
Apress / ￥ 4,123 ()
在庫あり。
ビジュアライジング・データ読んどけば良い気もするが。
ビジュアライジング・データ ―Processingによる情報視覚化手法matplotlibで格子のパターンをつくる
ちょっと描きたいものがあったので調べたらpcolormeshを使えばいいだけだった。
from numpy import *
from pylab import *
a = rand(100,100)
pcolormesh(a)
colorbar()
savefig('colour.png')
matplotlibでレーダーチャート
五角形にしたかったのでrule of fiveにPolar Surface Areaを加えておいた。
#!/usr/bin/env python
from matplotlib.projections.polar import PolarAxes
from matplotlib.projections import register_projection
from pylab import *
class RadarAxes(PolarAxes):
"""Class for creating a radar chart (a.k.a. a spider or star chart)
http://en.wikipedia.org/wiki/Radar_chart
"""
name = 'radar'
# use 1 line segment to connect specified points
RESOLUTION = 1
def draw_frame(self, x0, y0, r):
verts = [(r*cos(t) + x0, r*sin(t) + y0) for t in theta]
return Polygon(verts, closed=True)
def set_varlabels(self, labels):
self.set_thetagrids(theta * 180/pi, labels)
def get_axes_patch(self):
x0, y0 = (0.5, 0.5)
r = 0.5
return self.draw_frame(x0, y0, r)
if __name__ == '__main__':
register_projection(RadarAxes)
N = 5
theta = 2*pi * linspace(0, 1, N+1)[:-1]
theta += pi/2
labels = ['HBA', 'HBD', 'cLogP', 'MWT', 'PSA']
rule_of_five = [10, 5, 5, 500, 140]
desc = [12, 3, 3.6, 532, 160]
desc_rate = [100*desc[i]/float(v) for (i,v) in enumerate(rule_of_five)]
ax = subplot(111, projection='radar')
ax.fill(theta, [100]*N)
ax.fill(theta, desc_rate)
for patch in ax.patches:
patch.set_alpha(0.5)
ax.set_varlabels(labels)
rgrids((20, 40, 60, 80, 100))
grid(True)
show()
rule of fiveのようにある範囲内に収まっていること(超えるとリスク)というような指標を表すのにレーダーチャートは適しているんだろうか。つまり充足している事を示すような面積の表現はいいのかなぁ。あと、レンジが負になったりするのでそれもどうかと思う。
再考の余地はあるな。
TODO:多変量がある決まったレンジ内に収まっているかどうかを視覚的に捉えやすい表現手段を探す。