「Pythonによるデータ分析入門」を読む NumPy編
「Pythonによるデータ分析入門」を読んでいます。
pandasとNumpyについて書かれた本で、pandas目当てで購入したのですが、案外NumPy周りで知らないことが多く勉強になったのでメモです。
ここで挙げたもの以外にも、線形代数関連の関数や、ユニバーサル関数などがありますが、数が多いのでこの記事には書きません。
ぶっちゃけpandas使っているひとにとっては不要な情報かも。
なお、pandas編の予定は未定。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (7件) を見る
NumPyのデータ形式「ndarray」
NumPyのデータ形式の基本は、ndarrayです。
ndarrayは、「同じ型」のデータを連続するメモリ上に配置することで、効率的なデータアクセスを可能にします。
ndarrayの構成要素
ndarrayを構成する要素は、以下の4つです。
- データへのポインタ
- データ型（dtype）
- データの形状（shape）
- ストライド（stride）
以下の例では、2x10x5次元の1のみからなる3次元配列のdtype, shape, strideを見ています。
In [8]: arr = np.ones((2, 10, 5)) In [9]: arr Out[9]: array([[[ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.]], [[ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.]]]) In [10]: arr.dtype Out[10]: dtype('float64') In [11]: arr.shape Out[11]: (2, 10, 5) In [12]: arr.strides Out[12]: (400, 40, 8)
ストライドについては、少し解説が必要かもしれません。
ストライドとは、「『隣』の要素にアクセスするために、ポインタをいくら動かせばいいか」を表す数値だと思えばいいと思います。
float64型は8バイトですので、配列上の隣接要素はメモリ空間では以下のような隔たりがあります。
- arr[0][0][0]とarr[1][0][0]の間の距離はfloat64が10 x 5要素あるため、8 x 10 x 5 = 400バイト
- arr[0][0][0]とarr[0][1][0]の間の距離はfloat64が5要素分あるため、8 x 5 = 40バイト
- arr[0][0][0]とarr[0][0][1]の間の距離はfloat64が1要素分あるため、8バイト
C型配列とF型配列
上記の例では、arr[0][0][0]とarr[0][0][1]がメモリ上では隣接する要素でしたが、逆の場合も考えられます。つまり、arr[0][0][0]とarr[1][0][0]がメモリ上で隣接しているという形です。
前者のようなデータの保持方法を、C型順序といい、行優先ともいいます。「行優先」の由来は、2次元配列の場合だと、行ごとのデータがメモリ上に固まって格納されているためです。
一方で、後者のようなデータの保持方法を、F型順序といい、列優先ともいいます。
C型、F型という名称の由来は、C型がC言語で使われていた方式で、F型がFortranで使われていた方式だからだそうです。
In [61]: f_arr = np.ones((2,10,5), order='F') In [62]: f_arr.strides Out[62]: (8, 16, 160)
F型では、ストライドがC型とは異なるということがわかりますね。
C型とF型のパフォーマンス比較？
C型とF型の違いがもたらすパフォーマンス上の違いとは何でしょうか？
C型では、特定の行のデータはメモリ空間上のまとまった場所に格納されていますが、F型では、特定の行のデータはメモリ空間上に（規則的にではありますが、）分散しているということになります。
ということは、「各行の集計を行う」という動作は、分散している要素をかき集める分、F型のほうが処理に時間がかかりそうですね。
In [79]: c_arr = np.ones((1000, 1000), order='C') In [80]: f_arr = np.ones((1000, 1000), order='F') In [81]: %timeit c_arr.sum(1) # sum(1)は各行で総和をとる関数です 1000 loops, best of 3: 477 µs per loop In [82]: %timeit f_arr.sum(1) 1000 loops, best of 3: 453 µs per loop
あれ、全然変わんないや。。。
ちなみに本では明確にC型のほうが早い結果が出ています。（P.433）
本書に記載されている範囲では、メモリ上の連続した空間へのアクセスを高速にできる要因の一つとして、CPUのキャッシュ階層が挙げられています。（空間的局所性というやつですかね）
実験環境が、仮想マシンなのですが、なにか関係あるんでしょうかね。
配列の型
基本的な型
ndarrayを生成するタイミングでdtypeを指定してあげると、型を指定できます。
In [1]: import numpy as np In [2]: arr1 = np.array([1,2,3], dtype=np.float64) In [3]: arr1 Out[3]: array([ 1., 2., 3.]) In [4]: arr1.dtype Out[4]: dtype('float64') In [5]: arr2 = np.array([1,2,3], dtype=np.int32) In [6]: arr2 Out[6]: array([1, 2, 3], dtype=int32) In [7]: arr2.dtype Out[7]: dtype('int32')
上記では、float64とint32を試していますが、文字列や真偽値なんかも使えます。文字列の場合は固定長です。
型のキャスト
astype関数で型を変更できます。返り値として新しいndarrayオブジェクトが返ってきます。元のndarrayには変更が加えられません。
ndarrayはメモリ上に密に配置されているデータですので、型を変えるためには新しくメモリの領域を確保する必要があるため、新しいオブジェクトが返ってくるということですね。
In [8]: arr_a = np.array([True, False, True, False]) In [9]: arr_a Out[9]: array([ True, False, True, False], dtype=bool) In [10]: arr_b = arr_a.astype(np.float64) In [11]: arr_b Out[11]: array([ 1., 0., 1., 0.])
構造化配列
ndarrayの生成時に指定するdtypeに、「構造」をもたせることができます。
In [53]: dtype = [('name', 'U32'), ('age', 'u1')] In [54]: arr = np.array([(u'井上喜久子', 17), (u'武内駿輔', 17)], dtype=dtype) In [55]: arr Out[55]: array([('井上喜久子', 17), ('武内駿輔', 17)], dtype=[('name', '<U32'), ('age', 'u1')])
上記の例では、dtypeとしてname（U32型、つまりUnicode文字列で32文字）とage（u1型、つまりunsignedの1バイトInteger型）からなる構造を指定しています。
これによって、以下のようなアクセス方法が可能になります。
In [57]: arr[0]['name'] Out[57]: '井上喜久子' In [58]: arr['name'][0] Out[58]: '井上喜久子' In [59]: arr['name'] Out[59]: array(['井上喜久子', '武内駿輔'], dtype='<U32') In [60]: arr['age'] Out[60]: array([17, 17], dtype=uint8)
また、dtypeとして、以下のような階層を持った形も指定出来ます。
In [62]:# xは3つのint64からなり、yは1つのint32からなる In [63]: dtype = [('x', np.int64, 3), ('y', np.int32)] In [64]: arr = np.zeros(4, dtype=dtype) In [65]: arr Out[65]: array([([0, 0, 0], 0), ([0, 0, 0], 0), ([0, 0, 0], 0), ([0, 0, 0], 0)], dtype=[('x', '<i8', (3,)), ('y', '<i4')]) In [66]: # xはaとbからなり、aは倍精度浮動小数点数、bは単精度浮動小数点数、yはint32 In [67]: dtype = [('x', [('a', 'f8'), ('b', 'f4')]),('y', np.int)] In [68]: arr = np.array([((1, 2), 5), ((3, 4), 6)], dtype=dtype) In [69]: arr Out[69]: array([((1.0, 2.0), 5), ((3.0, 4.0), 6)], dtype=[('x', [('a', '<f8'), ('b', '<f4')]), ('y', '<i8')])
ソート
直接ソート（in-place）
ndarrayのソートは、既存のndarray内のデータが置換される形で行われます。新しいndarrayを作ったりはしません。
In [80]: values = np.array([5, 0, 1, 3, 2]) In [81]: values.sort() In [82]: values Out[82]: array([0, 1, 2, 3, 5])
多次元配列の場合は、sortの引数で軸を指定することができます。
In [108]: arr0 = np.random.rand(5, 3) In [109]: arr1 = arr0.copy() In [110]: arr0 Out[110]: array([[ 0.7135679 , 0.6001063 , 0.96356225], [ 0.94294085, 0.09982734, 0.11977476], [ 0.37965608, 0.46243769, 0.08483628], [ 0.04989755, 0.07350422, 0.2098485 ], [ 0.16695968, 0.00744647, 0.89886347]]) In [111]: arr0.sort(0) # 各列内でソート In [112]: arr0 Out[112]: array([[ 0.04989755, 0.00744647, 0.08483628], [ 0.16695968, 0.07350422, 0.11977476], [ 0.37965608, 0.09982734, 0.2098485 ], [ 0.7135679 , 0.46243769, 0.89886347], [ 0.94294085, 0.6001063 , 0.96356225]]) In [113]: arr1.sort(1) # 各行内でソート In [114]: arr1 Out[114]: array([[ 0.6001063 , 0.7135679 , 0.96356225], [ 0.09982734, 0.11977476, 0.94294085], [ 0.08483628, 0.37965608, 0.46243769], [ 0.04989755, 0.07350422, 0.2098485 ], [ 0.00744647, 0.16695968, 0.89886347]])
間接ソート
内部のデータの順番には影響を与えないまま、ソートしたいことがよくありますが、それを実現するのが間接ソートです。
argsort
argsortを用いると、ソートした時のインデクサー（arrayの各要素への参照）が得られますので、それを用いることで、ソートの結果を得ることができます。
In [117]: values = np.array([5, 0, 1, 3, 2]) In [118]: indexer = values.argsort() # 間接ソート In [119]: indexer Out[119]: array([1, 2, 4, 3, 0]) In [120]: values[indexer] # インデックスによるndarray参照 Out[120]: array([0, 1, 2, 3, 5])
多次元配列の時は、インデクサーの使い方に注意が必要です。1次元の時のようなわかりやすい参照を返してくれるわけではないのです。
In [145]: arr Out[145]: array([[ 0.43923064, 0.68340278, 0.58572889], [ 0.45396293, 0.61659209, 0.25095145], [ 0.58301305, 0.66421161, 0.44026581], [ 0.47286063, 0.65877653, 0.73912894], [ 0.06365072, 0.08286274, 0.38831605]]) In [146]: indexer0 = arr.argsort(0) In [147]: indexer0 # 各列内でソートされたインデックスがndarrayで返ってくる。 Out[147]: array([[4, 4, 1], [0, 1, 4], [1, 3, 2], [3, 2, 0], [2, 0, 3]]) In [148]: arr[indexer0] #そのまま使っても意味がない。 Out[148]: array([[[ 0.06365072, 0.08286274, 0.38831605], [ 0.06365072, 0.08286274, 0.38831605], [ 0.45396293, 0.61659209, 0.25095145]], [[ 0.43923064, 0.68340278, 0.58572889], [ 0.45396293, 0.61659209, 0.25095145], [ 0.06365072, 0.08286274, 0.38831605]], [[ 0.45396293, 0.61659209, 0.25095145], [ 0.47286063, 0.65877653, 0.73912894], [ 0.58301305, 0.66421161, 0.44026581]], [[ 0.47286063, 0.65877653, 0.73912894], [ 0.58301305, 0.66421161, 0.44026581], [ 0.43923064, 0.68340278, 0.58572889]], [[ 0.58301305, 0.66421161, 0.44026581], [ 0.43923064, 0.68340278, 0.58572889], [ 0.47286063, 0.65877653, 0.73912894]]]) In [153]: arr[:,1] #arrの2列目のみ抜粋 Out[153]: array([ 0.68340278, 0.61659209, 0.66421161, 0.65877653, 0.08286274]) In [154]: arr[:, 1][indexer0[:, 1]] # 2列目のみのソート結果を確認する Out[154]: array([ 0.08286274, 0.61659209, 0.65877653, 0.66421161, 0.68340278])
lexsort
あと、間接ソートにはlexsortなんてのもあります。
姓名のデータがあるときに、姓で優先的に並べ替えた後に、各姓の内部で名の並べ替えを行う、なんていう時に使えます。
In [165]: first_name = np.array(['Bob', 'Jane', 'Steve', 'Bill', 'Barbara']) In [166]: last_name = np.array(['Jones', 'Arnold', 'Arnold', 'Jones', 'Walters']) In [168]: sorter = np.lexsort((first_name, last_name)) In [173]: for name in zip(first_name[sorter], last_name[sorter]): .....: print (name) .....: ('Jane', 'Arnold') ('Steve', 'Arnold') ('Bill', 'Jones') ('Bob', 'Jones') ('Barbara', 'Walters')
last_nameでのソート結果は保持したまま、first_nameでもソートするという動作が実現できています。lexsortでは、後の方に渡したndarrayのほうが優先的にソートされます。
ファイルの入出力
バイナリ形式
すごく簡単にndarrayの保存と読み込みができます。
注意スべき点としては、勝手に.npyという拡張子がつくので、読み込みの際は忘れないようにする、ということでしょうか。
In [174]: arr = np.arange(10) In [175]: np.save('test', arr) # 保存 In [177]: new_arr = np.load('test.npy') # 読み込み In [178]: new_arr Out[178]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
テキスト形式
他のプログラムや言語とデータを共有して使いたいとき、csvファイルやタブ区切りファイルに一時保管して渡すという方法は、なんだかんだで一番使う手段だと思います。
In [183]: arr = np.random.rand(3, 5) In [184]: arr Out[184]: array([[ 0.75277139, 0.11091738, 0.68009413, 0.21621293, 0.72897406], [ 0.80775297, 0.39010844, 0.90030189, 0.25667751, 0.43382427], [ 0.64004754, 0.1824048 , 0.95000589, 0.29878554, 0.6953434 ]]) In [185]: np.savetxt('test.csv', arr, delimiter=',') # 保存 In [186]: %cat test.csv 7.527713922362844201e-01,1.109173805522273293e-01,6.800941327360305877e-01,2.162129304535805874e-01,7.289740567048890174e-01 8.077529740181056406e-01,3.901084449734447679e-01,9.003018943256323459e-01,2.566775088487044387e-01,4.338242701296414205e-01 6.400475420605206134e-01,1.824048033472692731e-01,9.500058880204274026e-01,2.987855382194802845e-01,6.953433979262562126e-01 In [187]: np.loadtxt('test.csv', delimiter=',') # 読み込み Out[187]: array([[ 0.75277139, 0.11091738, 0.68009413, 0.21621293, 0.72897406], [ 0.80775297, 0.39010844, 0.90030189, 0.25667751, 0.43382427], [ 0.64004754, 0.1824048 , 0.95000589, 0.29878554, 0.6953434 ]])
メモリマップファイル形式
非常に大きなデータファイルのために、メモリ上にすべてを展開できないとき用の形式です。
1TBとかのデータだとメモリに読み込めないので、普段はディスクにおいておいて、必要なときに必要な部分だけ取り出して処理をするということができるようです。
以下の例では、float64の3万 x 3万のndarrayとして扱えるメモリマップファイルを作成しています。
In [188]: mmap = np.memmap('mymmap', dtype='float64', mode='w+', shape=(30000, 30000)) # 作成 In [189]: %ls mymmap -lh -rw-rw-r-- 1 taromaru taromaru 6.8G 3月 1 13:49 mymmap In [190]: mmap Out[190]: memmap([[ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.], ..., [ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.], [ 0., 0., 0., ..., 0., 0., 0.]]) In [211]: mmap[:5] = 1 In [212]: mmap.flush() In [213]: del mmap # ipythonを再起動 In [3]: mmap = np.memmap('mymmap', dtype='float64', shape=(30000, 30000)) # 読み込み In [4]: mmap Out[4]: memmap([[ 1., 1., 1., ..., 1., 1., 1.], [ 1., 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., 0., ..., 0., 0., 0.]])
作成した時点で、6.8GBの巨大なファイルが作成されていることが確認できます。
dtypeで型が指定され、shapeも指定されているため、その分の領域をあらかじめ確保しているということです。
また、作成したmmapはflushで明示的にディスクへ書き込むことができます。
（裏でも、定期的に書き込みが発生しているようです。）
読み込みを行うときも、dtypeとshapeを指定しなくてはなりません。メモリマップファイルはただのバイナリファイルで、型や形状に関する情報を一切持っていません。
まとめ
ぼくがNumPyで普段使っているのは線形代数関連のほんの一部です。
ある程度体系的に勉強すると、こういうの学生時代に勉強したなーということが出てきたり、こんな便利な機能があったのかよ、などいろいろ発見があって、まあ楽しい。
