遺伝的アルゴリズムのメモ

最適化

基本

  • 1975にHollandが提案．
  • 最適化対象の関数の微分可能性や単峰性を要求しない（SAも同様）．
  • 但し，評価関数の順序性と，探索空間が位相を持っていることが要請される．
  • 解の遺伝子表現方法を決めるということは，表現型と遺伝子型の符号化方法(コーディング)を決めていると言える．
  • コーディング方法としては，バイナリコーディング（ナップザック問題のように，あるアイテムのあるなし，または個数を符号とする），順列コーディング（TSPやタスクスケジューリングなどのような，何かしらの並びを求める場合），実数値コーディング（実数を文字列としてコーディングする），また遺伝的プログラミングでは木（またはグラフ）としてコーディングする．
  • 探索の集中化（良い解の周辺を探索したい）と多様化（局所解からの脱出）をどう塩梅を取るか，という観点でSAやPSOなどと比べて見ると良い．
  • 集中化の前提としては，良い解は類似の構造を持つ，という近接最適性原理を仮定している．これが成り立たないような問題は，結局全探索するよりほか無いということを意味していて，工学の対象としては成立しない部類何だと思う．
  • 近傍解というもののしっかりとした定義は無い．そして，近傍解というものを適切に決めることは実は難しい．但し，近傍解の生成ルールを複数用意しておき，適用するルールを変化させることで局所構造への捕捉を回避することがよく行われる．
  • GAは理論的根拠が弱いとされているが，スキーマ定理と呼ばれる定理が一定の条件下では拠り所となる．スキーマとは部分構造のことであり，スキーマ定理はある部分構造の存在確率に関する定理であり，適応度の高い，ある一定の大きさの部分構造が残りやすい，ということを述べている．つまり，GAは交叉によりバンバンジャンプするけど，その進化は最適解が持っているであろう良い部分構造を内包するものが生き残り，最終的な解にたどり着く，と考えるんだ．だから，GAの適用が向く問題というのは，とても広い範囲をくまなく探索する必要な問題で，かつ，探索範囲は広いものの，部分構造としては最適解は准最適解類似の構造を内包している，というような仮定が成り立つ問題に向くんだと思う．

他の集団最適化アルゴリズムとの違い

  • DE(差分進化)は全体で最適解の情報を共有して，その全体の最適解に収束するような操作が入っており，探索序盤は広く，探索終盤には自然と最適解の周囲を集中的に探索したい，という目的を実現する意図が入っている．そして，それぞれは自分の周りをある程度しっかりと調べる．
  • 対してGAは大域的ジャンプをし続けてしまうというのが欠点と言える．本当に大域的な情報を調べたい，というときがGAの出番なのかな．DEは基本的には局所探索をバラバラにやっているイメージで，バンバンジャンプしていこう，ってイメージではないな．じんわりじんわりみんなで，って感じかな．

フロー

  1. 評価集団の生成
  2. 適応度計算
  3. 交叉候補の選択

　適応度に応じて交叉させる．交叉候補の選択方法は適応度に重みづけてルーレット選択する，もしくはN個をランダムに取ってきてトーナメントして強いものを残す，などが行われる．また，エリート保存と言い，最も良い上位何個かを必ず選択する，などもある． 4. 交叉 　　 　ある１点で切って，そこをスワップする１点交叉が基本．多点交叉もあるが，わざわざそうする理由があまり無いのか（１点交叉を何回も回すのと何が違う）？と思う．一様交叉というのは，どこかから後ろを全部入れ替えるとかではなく，マスクビットのようなもので，マスクしたビットを全部入れ替える，みたいなもの．これはマスクに知識を入れるなどである程度意味がありそうだ． 5. 突然変異 　 　交叉を大域ジャンプだとすると，突然変異は近傍解生成と言える．一般的には変数ごとに，今の変数の値を平均とするガウス分布に従うように乱数で変化させることが行われる．これをガウス突然変異というが，ガウス突然変異は裾の減速が激しいのでそれを緩めるためにコーシー分布に従う乱数で変化させるものもある．それがコーシー突然変異．

ヒストグラムの階級幅の決め方

統計・機械学習
  • 階級幅は基本的に等間隔．非等間隔だとどうしてもそこに恣意性が入る
  • 階級幅の決め方は，まず，幅に関する自明な知見，例えば学年などがあるならそれを採用する．
  • 次に，データから幅を決める場合，基本的には幅を決めるというよりはビン数を決める．データ全体の幅をR=Max-Minが求まるから，幅を一定だとすると，あとはビンの数が決まれば自動的に決まるのだ．
  • じゃあ，そのビン数はどう決めるか？実はしっかりとした理論的根拠のあるものはないようだ．殆どが経験則で，有名なものは下記がある．
  • データ数nに対してsqrt(n)
  • n<200だとsturgess法が有効と言われている．これは，幅hをh=1+ln(n)とする方法．
  • Freedmann and Diaconis法はh=2IQR(x)/n1/3, IQRは４分位範囲．
  • Scott法は3.5sigma/n1/3
  • pythonのmatplotlibのhistは内部ではnumpyのhistogramを使っていて，それは指定しないとビン数10として，幅を当分割する．

pandas@python

統計・機械学習 python

概要

  • pandasというのをよく目にする．Rっぽくデータ処理するためのライブラリみたいだ．
  • DataFrameに2次元配列を渡してデータを生成する．タプルの配列や，ハッシュ配列でも行ける．
  • 行はindex, 列はcolumnsで名前を付ければ，名前でアクセス可能．行は意味が無いなら名前つけない．その方がスライシングが分かりやすい．
  • ix[row, col]でアクセスするのが一番汎用性が高い．
  • 連結はjoinするような場合以外はconcatが柔軟で良い．
  • CSVは特に悩むところはないのでWeb参考にする．
import pandas as pd

##############################
# data setting
## from array (numpy array is also o.k.)
df = pd.DataFrame( [[0,1,2], [3,4,5], [6,7,8], [9,10,11]] ) # 4x3 table
df.columns = ['c0', 'c1', 'c2']
df.index = ['r0', 'r1','r2','r3']
df.shape   #=> (4,3)

## from array of hash
### The 'NaN' is filled automatically
# c0  c1  c2
# 1   2   nan
# nan 1   2
# 2   nan 1
df1 = pd.DataFrame( [{'c0':1, 'c1':2}, {'c1':1, 'c2':2}, {'c0':1, 'c2':2}] )

##############################
# row/column process
## reference
df[1:3]          # row slicing
df['c0']         # column exraction
df[['c0', 'c1']] # multi column extraction
df.ix[:,:]       # ix[row, col] 
## filter
df['c0'] > 0     # [False, True, True, True]
df[df['c0'] > 0] # filter by the above boolean array
## add/delete row/column
df['c3'] = [6,6,6,6]        # add column
print df.drop('c3', axis=1) # delete column(non-destructive)

##############################
# NaN process (non-destructive)
df1.fillna(0)      # NaN->0 
df1.interpolate()  # interpolatraion
df1.dropna(axis=0) # axis=0:row, axis=1 column 
print df1

データフレームの生成

  • 配列の配列として渡すと，行と列が逆になってしまう，などがあって，そのたびに調べてしまう。個人的に一番覚えやすいのは，下記の例のように辞書として渡す方法。これだと，列がキーになって明確。下の例はキーが'x', 'y'の列が出来る。
x = [1, 2, 3]
y = ['a', 'b', 'c']
df = pd.DataFrame({'x': x, 'y': y})

多重のソート

  • あるキーでソートして，その結果の中で別のあるキーでソートする，と言う場合。ソートしたいキーをリストで渡すだけ。便利だ。
import pandas as pd

x = [1, 3, 1, 2, 2]
y = [5, 1, 2, 4, 3]
names = list('abced')
df = pd.DataFrame({'x':x, 'y':y, 'name':names})
df.sort_values(['x'])
df.sort_values(['x', 'y'])

よくハマること

文字列の数値への変換

  • apply(pd.to_numeric, args=('coerce',))だと変換できなかったセルにはNaNが入る．で．それに気付かずにdropnaすると行全体がドロップしてしまう．
  • convert_objects(convert_numeric=True)だとそれはないんだけど，古いAPIだぞ，ってワーニングが出る．

dropnaに関して

  • dropnaするとインデックスが歯抜けになる．それを忘れて，乱数でアクセス行を作ってサンプリングしちゃうとNaNな行にアクセスしてしまう．

コピーに関して

  • データフレームの代入は浅いコピーがされる．つまりdf=df_orgとかして，dfをいじるとdf_orgも弄られてしまうので注意．

scikit-learnのメモ

全体通して

付属のデータセットに関して

  • 分類問題が2つ(iris:虫, digits：手書き文字), 回帰問題が3つ(boston：家の価格, diabetes：糖尿病, linnerud：心理状態)入っている．
  • それぞれが，dataでx(サンプル数x特徴量数の2次元配列，目標値は入らない)を, targetでyを取り出せる．digitsなどはimageで画像としても取れる．

パラメタ探索

  • 各アルゴリズムは設定パラメタを持っている．grid searchやcross validationで機械的に適切なパラメタを設定することができることもある．

Numpyのメモ

Numpyだけでも色々とメモが多いので別にする．

よく忘れること

  • scipyとの役割分担は？というのは，BLASとLAPACの関係．Numpyが基本的な線形代数演算(BLAS)で，それを使って科学技術計算(FFTとか)をするのがscipy(LAPACK)．
# ゼロ要素を除いたデータの和
d = numpy.array([1,2,0,3,4,0,5,0,6,7,8,0,9])
dd = numpy.delete(d, numpy.where(d == 0))
dd.sum()

NULL, NA, NaNの整理

基本事項

  • NULL, NA, NaNはそれぞれ意味が違う．NULLは無いこと（未定義），NAは使えないこと(Not Available，欠損値), NaNは数値でないこと(Not a Number，非数)だということを踏まえておく．

各言語での扱い

R

  • 全てあるので特に悩まない．
  • 0/0はNaN．
x = NULL
is.null(x) #=> True
x = Nan
is.nan(x)  #=> True
x = NA
is.na(x)   #=> True

python

  • まとめると，普通はNone，numpy,pandas使う時はnanを使っても良い，という感じかな？
  • 標準ライブラリではNoneがNullに相当する．Sqlite3にNULLを入れたいときはNoneを渡す．
  • numpyにはnanがある．masked arrayというのもあるみたいだけど使っているのかな？
  • scikit-learnには何もないので，別途自分で整理する必要がある
  • pandasはNaNがある．
# sqlite3使う場合には，タプルで渡すことに注意
db.execute("INSERT INTO table_name VALUES('hoge', ?, 1)", (None,))

SQLite3

  • NULLを使う．pythonから叩く場合にはNone(python)とNULL(sqlite)が相互に変換される．

scipy.spatial.Voronoiの使い方

問題設定

点群(points cloud)が与えられた時に，各点のボロノイ多面体の体積を求めたい．

pythonの例

  • scipyにあるspatial.Voronoiを使うとボロノイ多面体を簡単に計算してくれる．spatial.VoronoiはQhullという凸解析ライブラリのラッパのようだ．
  • 但し，spatial.Voronoiには各ボロノイ多面体の体積を計算するアルゴリズムは無いようだ．

多面体の体積の計算

  • ボロノイ多面体の体積を計算するためのアルゴリズムは無いが，各ボロノイ領域の頂点集合は計算されているため，それらの頂点集合を多面体の頂点集合として多面体の体積が求まればボロノイ多面体の体積を計算できる．
  • stackoverflowにstackoverflow.com同様の質問があったので参考にする．計算方法はドロネー三角分割で，四面体に分割して，各四面体の体積の和を取るというもの（たぶん）．
# 
import numpy as np
from scipy.spatial import Delaunay
from scipy.spatial import ConvexHull


def tetrahedron_volume(a, b, c, d):
    return np.abs(np.einsum('ij,ij->i', a-d, np.cross(b-d, c-d))) / 6

def convex_hull_volume(pts):
    ch = ConvexHull(pts)
    dt = Delaunay(pts[ch.vertices])
    tets = dt.points[dt.simplices]
    return np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1],
                                     tets[:, 2], tets[:, 3]))


#pts = np.random.rand(10, 3)
pts = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0],\
                [0,0,1], [1,0,1], [1,1,1], [0,1,1],\
                [0.5,0.5,2.0] ]           )

dt = Delaunay(pts)
tets = dt.points[dt.simplices]
print dt.simplices
print tets
vol = np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1], 
                                tets[:, 2], tets[:, 3]))
print pts
print vol
print convex_hull_volume(pts)

from scipy.spatial import Voronoi,Delaunay
import numpy as np
import matplotlib.pyplot as plt

def tetravol(a,b,c,d):
 '''Calculates the volume of a tetrahedron, given vertices a,b,c and d (triplets)'''
 tetravol=abs(np.dot((a-d),np.cross((b-d),(c-d))))/6
 return tetravol

def vol(vor,p):
 '''Calculate volume of 3d Voronoi cell based on point p. Voronoi diagram is passed in v.'''
 dpoints=[]
 vol=0
 for v in vor.regions[vor.point_region[p]]:
  dpoints.append(list(vor.vertices[v]))
 tri=Delaunay(np.array(dpoints))
 for simplex in tri.simplices:
  vol+=tetravol(np.array(dpoints[simplex[0]]),np.array(dpoints[simplex[1]]),np.array(dpoints[simplex[2]]),np.array(dpoints[simplex[3]]))
 return vol

x= [np.random.random() for i in xrange(50)]
y= [np.random.random() for i in xrange(50)]
z= [np.random.random() for i in xrange(50)]
dpoints=[]
points=zip(x,y,z)
print points
vor=Voronoi(points)
vtot=0


for i,p in enumerate(vor.points):
 out=False
 for v in vor.regions[vor.point_region[i]]:
  if v<=-1: #a point index of -1 is returned if the vertex is outside the Vornoi diagram, in this application these should be ignorable edge-cases
   out=True
 if not out:
  pvol=vol(vor,i)
  vtot+=pvol
  print "point "+str(i)+" with coordinates "+str(p)+" has volume "+str(pvol)

print "total volume= "+str(vtot)