Numpy 入門
(機械学習勉強会 in 新潟 2019/03/23)
目次
1 はじめに
2 numpy とは何か
Python の数値計算のためのライブラリです。 numpy によってできることとメリットは以下になります
No. | ライブラリ | 説明 | できること、メリット |
---|---|---|---|
1 | numpy | 数値計算ライブラリ | 1. 行列の計算を行うコードを短く書ける |
2. 行列の計算、ベクトルの計算の高速化( 速度比較 参照) |
3 numpy モジュールのインポート方法について
numpy を使用する場合、numpy モジュールのインポートは以下のように行います。
コード1. 「モジュールのインポート」:
import numpy as np
これにより、'np' の2文字で numpy を使用することができます。
コード2. 「モジュールの使用方法」:
import numpy as np arr = np.array([1,2,3]) arr
Out: array([1, 2, 3])
使用する名前は numpy など np 以外の名前にもできますが、 np と書くのが一般的になっているためこちらを使用しましょう。
4 数値を作成する
はじめに、numpy で数値を作成しましょう。数値は、スカラ(Scalar)と呼ばれることがあります。 これらは、後述する配列の要素( 5.1.1 参照)となりますので、しっかり覚えましょう。 np.int32 を使用すると numpy の整数を作成できます。
コード3. 「numpy で数値の作成」:
import numpy as np answer = np.int32(100) answer
Out: 100
4.1 np.int32 とは
np.int32 はデータタイプ(dtype)として定義されています。 np.int32(100) を実行すると np.int32 データタイプのオブジェクトが作成されます。 np.int32 以外にも数値を表わすデータタイプが用意されています。
データタイプ | 数値の種類 | 用途 |
---|---|---|
np.uint8 | 8bit 符号無し整数 | グレースケール画像 |
np.int32 | 32bit 整数 | 画像処理 |
np.float32 | 32bit 浮動少数 | 画像処理、平均値、確率 |
np.float64 | 64bit 浮動少数 | 数値計算(精度が求められるもの) |
… |
真偽値(True/False)を現す場合は np.bool データタイプが用意されています。
5 「1次元配列」の数値計算をする
配列の基本となる1次元配について学習しましょう。
5.1 「1次元配列」を作成する
1次元の配列を作成しましょう。1次元配列とは以下のようなものです。配列は'要素'と呼ばれるもので構成されています。
値1 | 値2 | 値3 | ,,, |
5.1.1 配列の'要素'とは
一般的に、上記 5.1 のような配列の中にある値1や値2などのことを配列の”要素”(Element)と呼びます。 1次元配列だけでなく、2次元配列や3次元配列も同様です。配列は要素によって構成されています。
※ numpy では配列を作成する関数の引数に dtype を指定することで、その配列の要素のデータタイプを指定できます。 引数 dtype を指定しない場合は、自動的に配列のデータタイプが決定されます。 numpy では配列は同一のデータタイプの要素から構成されます。
5.1.2 np.array 関数
np.array 関数を使用すると、標準 Python のリストを配列に変換することができます.
コード4. 「標準 Python のリストを配列に変換する」:
import numpy as np arr = np.array([1, 2, 3]) arr
Out: array([1, 2, 3])
5.1.3 np.arange 関数
np.arange 関数を使用すると、数値の1次元配列を作成することができます。
コード5. 「np.arange 関数を使用し、数値の1次元配列を作成」:
import numpy as np arr = np.arange(3) arr
Out: array([0, 1, 2])
5.2 「1次元配列」の形状 shape について
配列には shape と呼ばれる、配列の形状を表す情報があります。 配列の形状 shape を調べるには以下のようにします。 この例の場合、配列は長さ3の1次元配列であることがわかります。
コード6. 「配列の形状 shape を調べる」:
import numpy as np arr = np.arange(3) answer = arr.shape answer
Out: (3,)
5.3 「1次元配列」 + 「1次元配列」
1次元配列と1次元配列の四則演算の基本を学習します。 ここでは例として、1次元配列と1次元配列の足し算を行います。 以下のように、1次元配列を二つ作成し、足し算記号(+) を使用して求めます。 計算は、要素ごと(element wise)に行われます。
コード7. 「1次元配列と1次元配列の足し算」:
import numpy as np arr_1 = np.array([1, 2, 3]) arr_2 = np.array([0, 0, 1]) answer = arr_1 + arr_2 answer
Out: array([1, 2, 4])
その他の演算子(-*/など)についても、同じように、配列の要素ごとに計算されます。 形状の異なる配列同士の計算は、エラーになるか、ならない場合があります。 形状の異なる配列同士の計算については 11 を参照してください。
5.4 配列の要素に数学関数を適用する
1次元配列の要素ごと(element wise)に数学関数を適用してみましょう。 ここでは例として平方根を適用してみましょう。
コード8. 「配列の各要素に平方根を適用」:
import numpy as np arr = np.array([1, 2, 3]) answer = np.sqrt(arr) answer
Out: array([1. , 1.41421356, 1.73205081])
np.sqrt 関数以外にも様々な数学関数が用意されています。 それらの関数も、上記のように使用することで、要素ごとに適用することができます。 以下に数学関数の一部を記載します。
No. | 関数 | 説明 |
1 | np.exp | 指数関数 |
2 | np.log | 対数関数 |
3 | np.sin | sin関数 |
4 | np.arcsin | arcsin関数 |
5.5 1次元配列の要素の合計を求める
合計を求める場合は、 np.sum 関数を使用します。
コード9. 「np.sum 関数を使用して1次元配列の要素の合計を求める」:
import numpy as np arr = np.arange(11) answer = np.sum(arr) answer
Out: 55
6 「2次元配列」の数値計算をする
2次元配列について学習しましょう。
6.1 「2次元配列」を作成する
2次元の配列を作成しましょう。2次元配列とは以下のようなものです。
値11 | 値12 | 値13 | ,,, | |
値21 | 値22 | 値23 | ,,, | |
値31 | 値32 | 値33 | ,,, | |
,,, |
6.1.1 np.array 関数
np.array 関数を使用すると、2次元配列を作成することができます。(1次元配列の時と同様)
コード10. 「np.array 関数を使用して2次元配列を作成する」:
import numpy as np arr = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) arr
Out: array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
6.1.2 reshape メソッド
配列の reshape メソッドを使用すると配列の形状を変えることができます。 ここでは、 1次元配列の形状を変えて、2次元配列を作成してみましょう。
コード11. 「reshape メソッドを使用して1次元配列の形状を変えて、2次元配列を作成する」:
import numpy as np arr = np.arange(9) answer = arr.reshape((3, 3)) answer
Out: array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
6.2 「2次元配列」 + 「2次元配列」
2次元配列同士の四則演算、関数を適用した時の動作は、1次元配列の時と同様、要素ごとに計算が行われます。 (そのため2次元配列の四則演算については説明を省略します。)2次元配列を使用すると、行列の掛け算を行うことができます。
行列の掛け算とは、具体的には以下のような計算のことです。 行列 A と行列 B の掛け算(積)は A・B と書き、A と B と A・B は以下のような関係になります。 ※説明の簡略化のため行列の大きさは2x2としています。
行列 A
a b c d 行列 B
e f g h 行列 A・B
a*e+b*g a*f+b*h c*e+d*g c*f+d*h
行列の掛け算を行うには np.dot 関数を使用します。
コード12. 「np.dot 関数を使用して行列の掛け算を行う」:
import numpy as np arr_1 = np.arange(9).reshape((3, 3)) arr_2 = np.arange(9).reshape((3, 3)) answer = np.dot(arr_1, arr_2) answer
Out: array([[ 15, 18, 21], [ 42, 54, 66], [ 69, 90, 111]])
配列オブジェクトの dot メソッドを使用して以下のように書くこともできます。
コード13. 「dot メソッドを使用して行列の掛け算を行う」:
import numpy as np arr_1 = np.arange(9).reshape((3, 3)) arr_2 = np.arange(9).reshape((3, 3)) answer = arr_1.dot(arr_2) answer
Out: array([[ 15, 18, 21], [ 42, 54, 66], [ 69, 90, 111]])
注意: np.dot は '∗'による掛け算と同じではありません。 '∗' による掛け算は要素ごと(element wise)に掛け算が行わます。 一方 np.dot は行列の掛け算ですので、積和演算(掛け算と足し算)が行われます。
6.3 「2次元配列」の要素の合計を求める
'すべての要素の合計'を求める場合は、1次元配列の時と変わりませんので省略しますが、 2次元配列では、さらに、'列の合計'と'行の合計'を求めることができます。
6.3.1 列の合計
np.sum 関数の引数で axis=0 とすることで列の合計を求めることができます。
コード14. 「列の合計を求める」:
import numpy as np arr = np.arange(9).reshape((3, 3)) answer = np.sum(arr, axis=0) answer
Out: array([ 9, 12, 15])
6.3.2 行の合計
np.sum 関数の引数で axis=1 とすることで行の合計を求めることができます。
コード15. 「行の合計を求める」:
import numpy as np arr = np.arange(9).reshape((3, 3)) answer = np.sum(arr, axis=1) answer
Out: array([ 3, 12, 21])
6.3.3 配列の次元数を保つには
6.3.1 、 6.3.2 で np.sum を使用して求めた配列を見ると、2次元配列の入力に対して出力が1次元配列となっており、次元数が変わってしまうことに気付くかと思います。 np.sum 関数の引数で keepdims=True とすることで、配列の次元数を保つことができます。
コード16. 「np.sum 関数の引数で keepdims=True とし配列の構造をなるべく保つ」:
import numpy as np arr = np.arange(9).reshape((3, 3)) answer = np.sum(arr, axis=1, keepdims=True) answer
Out: array([[ 3], [12], [21]])
6.4 行列計算用の関数
2次元配列には行列計算用の関数を使用することができます。 行列計算用の関数は np.linalg にあります。 その一部を紹介します。
関数 | 説明 |
---|---|
np.linalg.inv | 逆行列を求めるための関数 |
np.linalg.norm | 行列のノルムを求めるための関数 |
7 配列を作成する関数
7.1 固定された値を要素にもつ配列を作成する関数
No. | 関数 | 説明 |
---|---|---|
1 | np.zeros | 全ての要素が 0 の配列を作成します |
2 | np.ones | 全ての要素が 1 の配列を作成します |
3 | np.full | 全ての要素が x の配列を作成します(x は引数で指定) |
コード17. 「np.zeros 関数を使用し全ての要素が 0 の配列を作成」:
import numpy as np arr = np.zeros((3, 3)) arr
Out: array([[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]])
7.2 単位行列を作成する関数
7.3 ランダムな値を要素にもつ配列を作成する関数
No. | 関数 | 説明 |
---|---|---|
1 | np.random.randint | ランダムな整数の配列を作成します |
2 | np.random.random | ランダムな少数の配列を作成します |
コード19. 「np.random.randint 関数を使用しランダムな整数の配列を作成」:
import numpy as np arr = np.random.randint(0, 10, (3, 3)) arr
Out: array([[4, 0, 7], [2, 3, 4], [3, 0, 5]])
コード20. 「np.random.random 関数を使用しランダムな少数の配列を作成」:
import numpy as np arr = np.random.random((3, 3)) arr
Out: array([[0.23963676, 0.14041946, 0.67624721], [0.75209081, 0.3182057 , 0.18579186], [0.0888395 , 0.87384273, 0.24740674]])
7.4 似たような配列を作成する関数
配列の計算をプログラミングしているときに、ある配列と同じ形状をもつ配列だが、値は違うという配列を作成したい時があります。そのような場合は、以下の関数を使用します。
No. | 関数 | 説明 |
---|---|---|
1 | np.zeros_like | 引数で指定した配列と同じ形状 |
& | ||
全ての要素が 0 の配列を作成します | ||
2 | np.ones_like | 引数で指定した配列と同じ形状 |
& | ||
全ての要素が 1 の配列を作成します |
コード21. 「np.zeros_like 関数を使用し形状が同じ配列を作成」:
import numpy as np arr = np.random.randint(0, 10, (3, 3)) arr_2 = np.zeros_like(arr) arr_2
Out[4]: array([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
7.5 ある範囲の数値を要素にもつ配列を作成する関数
No. | 関数 | 説明 |
---|---|---|
1 | np.linspace | 数直線の範囲を等幅で分割した1次元配列を作成します |
1変数関数を調べる時などに使用できます | ||
2 | np.meshgrid | 2次元座標などの座標を配列を使い表現するときに使用します。 |
2変数関数を調べる時などに使用できます |
コード22. 「np.linspace 関数の使用例」:
import numpy as np arr = np.linspace(0.0, 0.9, 10) arr
Out: array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
コード23. 「np.meshgrid 関数の使用例」:
import numpy as np x = np.linspace(0, 0.4, 5) y = np.linspace(0, 0.4, 5) x_mesh, y_mesh = np.meshgrid(x, y) answer = x_mesh, y_mesh answer
Out: (array([[0. , 0.1, 0.2, 0.3, 0.4], [0. , 0.1, 0.2, 0.3, 0.4], [0. , 0.1, 0.2, 0.3, 0.4], [0. , 0.1, 0.2, 0.3, 0.4], [0. , 0.1, 0.2, 0.3, 0.4]]), array([[0. , 0. , 0. , 0. , 0. ], [0.1, 0.1, 0.1, 0.1, 0.1], [0.2, 0.2, 0.2, 0.2, 0.2], [0.3, 0.3, 0.3, 0.3, 0.3], [0.4, 0.4, 0.4, 0.4, 0.4]]))
8 配列の要素へのアクセス
8.1 配列の要素の取得
添字(インデックス)を使用することで、配列の要素を取得することができます。
8.1.1 要素の取得の基本
以下のようにすると、配列の要素の値を取得することができます。
コード24. 「1次元配列の配列の要素の取得の基本」:
import numpy as np arr = np.array([1, 3, 5]) arr[2]
Out: 5
コード25. 「2次元配列の配列の要素の取得の基本(1)」:
import numpy as np arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]]) arr[2][2]
Out: 23
上記コードは以下のようにも書くことができます。
コード26. 「2次元配列の配列の要素の取得の基本(2)」:
import numpy as np arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]]) arr[2, 2]
Out: 23
8.1.2 要素を範囲で指定して取得
以下のようにすると、指定した範囲にある配列の要素を取得することができます。
コード27. 「1次元配列の要素を範囲で指定して取得」:
import numpy as np arr = np.array([1, 3, 5]) arr[1:]
Out: array([3, 5])
コード28. 「2次元配列の要素を範囲で指定して取得」:
import numpy as np arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]]) arr[1:, 1:]
Out: array([[11, 13], [21, 23]])
8.1.3 逆順の配列の取得
8.2 配列の要素の書き換え
8.2.1 要素の書き換えの基本
配列の要素を添字(インデックス)で指定して別の値に書き換えることができます。 以下のようにすると、配列の要素の値を書き換えることができます。
コード30. 「1次元配列の配列の要素の書き換えの基本」:
import numpy as np arr = np.array([1, 3, 5]) arr[2] = 100 arr
Out: array([ 1, 3, 100])
コード31. 「2次元配列の配列の要素の書き換えの基本(1)」:
import numpy as np arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]]) arr[2][2] = 100 arr
Out: array([[ 1, 3, 5], [ 7, 11, 13], [ 17, 21, 100]])
上記コードは以下のようにも書くことができます。
コード32. 「2次元配列の配列の要素の書き換えの基本(2)」:
import numpy as np arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]]) arr[2, 2] = 100 arr
Out: array([[ 1, 3, 5], [ 7, 11, 13], [ 17, 21, 100]])
8.2.2 要素を範囲で指定して書き換え
以下のようにすると、指定した範囲にある配列の要素を書き換えることができます。
コード33. 「1次元配列の要素を範囲で指定して書き換え」:
import numpy as np arr = np.array([1, 3, 5]) arr[1:] = 100 arr
Out: array([ 1, 100, 100])
コード34. 「2次元配列の要素を範囲で指定して書き換え」:
import numpy as np arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]]) arr[1:, 1:] = 0 arr
Out: array([[ 1, 3, 5], [ 7, 0, 0], [17, 0, 0]])
上記コードの arr[1:] や arr[1:, 1:] は配列のコピーではないことに注意しましょう。 このようなものを numpy ではビュー(View)と呼びます。実態は arr のデータを参照しているものです。 このような ビュー arr[1:] や arr[1:, 1:] に対する値の代入は、ビューの参照先である arr に対する値の代入であることを意味します。
よくあるのが、値を保持しておきたい配列のビューに代入を行なってしまいその配列を破壊してしまうことです。 値を代入したい場合は、以下のように、配列のコピーを作成して、そのコピーに対して代入しましょう。
コード35. 「配列のコピーを作成して、そのコピーに対して代入する」:
import numpy as np arr_orig = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]]) arr = arr_orig.copy() arr[1:, 1:] = 0 arr_orig
Out: array([[ 1, 3, 5], [ 7, 11, 13], [17, 21, 23]])
9 真偽値(True/False)について
'==' や '<' などの比較演算子を使用すると配列の要素ごとに比較が行われます。
コード36. 「比較演算子による配列の比較」:
import numpy as np arr_1 = np.arange(9).reshape((3, 3))**2 answer = arr_1 >= 10 answer
Out: array([[False, False, False], [False, True, True], [ True, True, True]])
この真偽値を配列のインデックスに指定すると、 真偽値のTrue の箇所の要素だけ集めた1次元配列を求めることができます。
コード37. 「真偽値の配列による要素の収集」:
import numpy as np arr_1 = np.arange(9).reshape((3, 3))**2 answer = arr_1[arr_1 >= 10] answer
Out: array([16, 25, 36, 49, 64])
以下のように、 np.where を使用すると値が True のインデックスを求めることができます。 np.where の戻り値は、「インデックスを要素にもつ配列」を要素にもつタプル(tuple)です。
コード38. 「np.where による条件を満たす配列のインデックスの取得」:
import numpy as np arr_1 = np.arange(9).reshape((3, 3))**2 answer = np.where(arr_1 >= 10) answer
Out: (array([1, 1, 2, 2, 2]), array([1, 2, 0, 1, 2]))
10 配列の制約
11 配列の長さや、形状が異なる場合に行われる計算について
形状の異なる配列同士の演算を行う場合、それ専用の規則(broadcasting rule)に従い処理が行われます。 その規則に従いエラーになることもあれば計算結果が求まることもあります。
11.1 エラーになるか計算結果が求まるかの判定方法
配列 arr1 と arr2 のそれぞれの配列の shape を一番右から左にかけて以下条件を満たす時、計算結果が求まります。 条件を満たさない場合はエラーとなります。また、エラーにならず、形状が異なる場合は、お互いに異なる形状を補うよう値が補完され、計算されます(詳細は割愛します)。
(以下は「Broadcasting」(https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) より引用)
- 両方の値が同じ、または、
- 片方の数値が1
11.2 Broadcasting rule の適用例
OK
配列 shape arr1 10 arr2 1 結果 10 OK
配列 shape arr1 10 1 arr2 1 10 結果 10 10 OK
配列 shape arr1 10 10 1 arr2 1 1 3 結果 10 10 3 NG
配列 shape arr1 10 10 3 arr2 1 1 2 結果 エラー
12 速度比較
numpy を使用すると、どれくらい速くなるのでしょうか。 100x100の要素からなる2次元配列の各要素を +1 する処理を numpy を使用する前と後で処理時間を比較してみました。 処理時間の計測には、マジックコマンド %timeit を使用しています。 その結果、numpy の使用により処理速度は 100倍 になりました。また可読性も向上したと思います。
コード39. 「速度比較 numpy 使用前」:
def func(x): y = x.copy() for i in range(100): for j in range(100): y[i][j] += 1 return y # リストの初期化 x = [[0 for i in range(100)] for j in range(100)] # 処理時間計測 %timeit func(x)
785 µs ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
コード40. 「速度比較 numpy 使用後」:
import numpy as np def func(x): y = x.copy() y += 1 return y # 配列の初期化 x = np.zeros((100, 100)) # 処理時間計測 %timeit func(x)
7.11 µs ± 284 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
13 まとめ
配列の作成方法、配列同士の計算や、配列への関数の適用方法について学習しました。
numpy を使って得られる恩恵は以下の2つです。
- 行列計算のコード量が少なくなる
- 処理速度が早くなる
14 最後に
既存の行列の計算のソースコードの可読性の向上や、処理速度向上などを目的に、 numpy を使用したいと思われるかもしれません。 まず計算の対象が 10 の制約の範囲内で処理が可能であることを確認しましょう。 次に、以下の表のように可読性と、処理速度を天秤にかけ、 numpy で実装するか標準 Python で実装するか考えると良いのではないかと思います。 (numpy を使ったからといって常に「可読性が高い&処理速度が早い」コードになる訳ではありません。)
可読性 | 処理速度 | numpy を使用する? |
---|---|---|
高くなる | 高くなる | YES |
高くなる | 低くなる | YES or NO |
低くなる | 高くなる | YES or NO |
低くなる | 低くなる | NO |
numpy を使用してコーディングするとテクニカルで難解なコードになってしまうことがあります。 可読性重視であればテクニカルなコーディングはしないほうが良いでしょう。 numpy を使用してコーディングする際は、numpy を使用する当初の目的を見失わないようにしましょう。
それでは、よい numpy ライフを!
Created: 2019-03-23 Sat 11:25
SciPy Japan 2019 に参加します
Wikipediaの文書を使って文の長さの平均を求める
Wikipediaの文書を使って文の長さの平均値と最頻値を求めました。機械学習勉強会 in 新潟で発表した内容です。Wikipedia の約60億文字以上のXMLファイルから、独自に定義した文を抜き出し、その文に対して統計量を求めました。使用言語はいつも通りPythonです。
1. 動作環境について
コードは以下の環境で動作確認を行っています。
- CPU:6コア
- メモリ:256GB
- HDD:50GB
※環境はGCE(Google Compute Engine)で行っています。実行する場合はメモリ量に注意してください。
2. Wikipedia 文書のアーカイブについて
Wikipedia の文書は以下からダウンロードできます。 "jawiki-latest-pages-articles.xml.bz2"の中にあるXMLが今回解析に使用したファイルになります。
3. 発表スライド
発表スライドは、一番下に置いておきましたので、そちらをご覧ください。
4. やってみた感想
まず、このような分析は、自然言語処理の分野では頻度分析と呼ばれるらしいです。統計ですね。最近は統計のビジネス書をコツコツ読んで勉強していたたため、おかげでサンプルデータの分布を見た瞬間パラメタ推定できると気づくことができました。
パラメタ推定は本で読んでいたけど実際やってみると、ニューラルネットワークなど必要とせず非常にシンプルなのですが、強力でした。今回は、分布が対数正規分布にフィットしたのですがその理由を考察するときに、仮定した関数の性質が役に立ちました。結果、個人的に納得感のある考察を得ることができました(対数正規分布に従うのはなぜか?の考察です)。
最後ですが、Wikipediaを日常的に使用していたため、その中に法則性が存在するという事実に驚きました。おそらく解析しないと一生気付かなかったと思います。個人的に、その法則性の存在を客観的に示せたので非常に有意義だったと思います。また、求めた最頻値の値自体にも、価値があると思います。この値はWikipediaで遭遇する可能性が一番高い文の文字列の長さなので、それは、もしかすると我々が最も親しんでいる文字列長かもしれません。そのため24文字を目指して文を書くとわかりやすさの改善に繋がるつかもしれません。ビッグデータの時代ですが、手付かずのデータが山ほどあるとのことなので、まだ誰も知らない事実がデータの中に大量に眠っていると推測されます。今後それらのデータの中から、想像もしなかった事実が出てくることを期待しています。
5. 発表会で上がった質問や意見
発表会で上がった質問をここにメモしておきます。これらの質問や意見については、今後扱うかもしれませんし、扱わないかもしれません。個人的には、最頻値を求めて十分満足してしまったため、あとやったら面白そうなことは、コーパスで教師なし学習(クラスタリングなど)をさせてテキストマイニングするくらいしか思いつかないです。あと人工無脳。
1. 他の分布で推定してみたらどうなるのか?例えばポアソン分布でやってみたらどうなるか?
→ ググったら一番最初に出てきたキーワードが対数正規分布だったのでそれでやってみた、ポアソン分布ではパラメタ推定していない。(ただし、対数正規分布でパラメタ推定を行ったきっかけはどうであれ、実際は、考察に記載したように対数正規分布を仮定するもっともらしい理由があるため、ポアソン分布や他の分布でパラメタ推定は行う必要性を感じなかった)
2. 集めた文(コーパスとも言える)を今後どういったことに使えるのか?
→ オートエンコーダによる教師なし学習(クラスタリング)を行いテキストマイニングを行うことが考えられる。
3. Wikipediaではなくニュースサイトなどの文の場合は違う分布になるのか?
→ 本などをみるとコーパスの分野が変わると違う分布になったり、あるコーパスで学習した機械学習のモデルが、別の分野のコーパスを入力して予測すると精度が落ちるなどの事例がある様子。そのため、ニュースサイトの分布とは異なるものになる可能性が高い。
4. XMLの読み込みにかかる時間はどれくらいか
→ 約15分くらい
5. 正規表現の適用にかかる時間はどれくらいか
→ 6コアCPUで分散処理しているため早い、約5分くらい。CPUを増やせばほぼリニアに早くなる(←発表では10分くらいと言ったかもしれませんが訂正します)
発表スライド
発表したスライドになります。正規表現が冗長になっており(スクエアブラケットの中のバックスラッシュのこと)、もっとコンパクトにできますが、そこは認識しています。時間がないため修正せずにあげています。ご指摘ありましたらコメントをお願いいたします。
以上。
Pythonで2次元データ作成ーポリゴンに収まる点の集合
前回の続きです。今回は、前回に引き続き、ポリゴンに収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それがポリゴンの内側であれば点を残す、ポリゴンの外側であれば点を消すを繰り返すだけです。「ポリゴンの内側か?」を判定する関数を作成し、その関数に点の位置を次々と問い合わせればよいでしょう。
※ここでは「ポリゴン」は多角形の意味で使っています。
前回の記事:
keimina.hatenablog.jp
前々回の記事:
keimina.hatenablog.jp
N角形(N>=3)の領域内の点の集合を得ることができますので、前回や前々回の記事で作成した点の集合も、この記事の方法だけで作成することができるようになります。
目次
1. 点がポリゴンの領域の内側にある点か判定する方法
点がポリゴンの領域の内側にあるか判定する方法は以下に色々書かれていました。
上記の回答の一つを見たところ、点がポリゴンの内側にあるとき、その点とポリゴンを構成する点から算出される角度の総和は必ず360度になるらしいです。図を書いてイメージをつかみます。以下が、点A,点B,点C,,,点Jを結んだポリゴンの内部に点Pがあった時の図になります。点Bから点Cあたりでガタガタと移動する箇所があります。点Bから点Cに移動する時のように角度θiがマイナスになるときもあるのですが、その後のガタガタで相殺されます。θiの総和が360°のとき左回りに1回転できたということになります。※かなりざっくりした説明で厳密ではありません。
一方、点がポリゴンの外側にあるとき、360度にはならないと書かれています。これも図を書いてイメージをつかみます。θiの箇所についてみて見ると前半は右回りで、後半は左回りとなっています。右回りと、左回りで相殺されますので、θiの総和は0°になることがわかります。
ポリゴンのような多角形の図形は、円の一つだと考えられます。そのため多角形を円とみなして考えるとわかりやすいです。円形の壁の中に、自分がいたとして、ぐるっと見回して壁に囲まれていたら、円の中だとわかります。一方、ぐるっと見回して左には壁があったけど右には壁がなかったなどの場合は壁の中でないことがわかります。算出される角度の総和が360°だったというのは、壁をたどると左周りに1回転できたということと同じです。1回転できなかった場合は、壁を辿っても1回転できなかったということですので、壁がつながっていない、すなわち壁の内側ではないということがわかります。
今回は本アルゴリズムを使用します。(アルゴリズムがちゃんと機能するのかも確認します)
2. 角度の求め方
角度を求めるには以下の数学公式を用います。
cosθ = (PA・PB)/(|PA|*|PB|)
θ = arccos(cosθ)
3. 角度を求めるための関数を作成する
numpyで角度を求める関数を実装すると以下になります。
コード:
import numpy as np from numpy import inner from numpy.linalg import norm def angle(P, A, B): ''' args ---- P: np.array(x, y)) A: np.array(x, y)) B: np.array(x, y)) returns ------- theta: x (degree) ''' # 全体をシフトさせる PA = A - P PB = B - P # cosθを求める costheta = inner(PA, PB)/norm(PA)/norm(PB) # 角度を求める(度) theta = np.arccos(costheta)/np.pi * 180 return theta # テスト P = np.array((0,0)) A = np.array((1,1)) B = np.array((1,0)) print(angle(P, A, B))
出力:
45.00000000000001
4. 時計回りか、反時計回りかを判定する
3. によって角度の絶対値はもとめることができました。しかし、反時計回りか時計回りか区別できません。時計回りか、反時計回りか区別できないと、角度の総和が常に360°になってしまうため、区別できるようにする必要があります。そこで、点A(a, b)と点B(c, d)と点P(0, 0)によってできるベクトルPAとPBのなす角θについて、時計回りと反時計回りを以下のように定義します。
(i) d > (b/a) * c
ならば、反時計回り方向の角度
(ii) d < (b/a) * c
ならば時計回り方向の角度
※説明を簡単にするため点P(0, 0)としています。点Pが(0, 0)以外のときでも、全体をシフトさせて点Pが(0, 0)になるように点Aと点Bを移動させても、角度に関する情報は失われないので、P(0, 0)でも問題ありません。
また(i)(ii)の不等号の式の両辺にaをかけて、数式を変形すると以下のように書き直せます。
(i)' a*d - b*c > 0
(ii)' a*d - b*c < 0
図で表すと以下になります。点Bが緑色の領域にあればPABの角度は反時計回り方向。緑色の領域にあれば、時計回り方向というふうに定義しています。
なぜ、この定義をするのかは、説明しづらいですが、図1、図2を実際に手でなぞって見ると、まず回転方向が必要なことがわかります。また0°<θ<180°のときは足し算で-180°<θ<0°のときは引き算のように計算方法を変えなければいけないこともわかると思います。そのため、このような定義が必要であるということがわかります。
5. 角度を求めるための関数2(反時計回り、時計回りの情報を追加)
3.で作成したangle関数を、4.の定義を使用して、反時計周りならプラスを、時計回りならマイナスの角度を返すように修正します。
コード:
import numpy as np from numpy import inner from numpy.linalg import norm def angle2(P, A, B): ''' args ---- P: np.array(x, y)) A: np.array(x, y)) B: np.array(x, y)) returns ------- theta: x (degree) ''' # 全体をシフトさせる PA = A - P PB = B - P # cosθを求める costheta = inner(PA, PB)/norm(PA)/norm(PB) # 角度を求める(ラジアン) theta = np.arccos(costheta) # 角度を求める(度) theta = np.arccos(costheta)/np.pi * 180 # if PA[0]*PB[1] - PA[1]*PB[0] > 0: # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス) pass else: # 条件(ii)'に該当するときはthetaをマイナスにする theta *= -1.0 return theta # テスト P = np.array((0,0)) A = np.array((1,1)) B = np.array((1,0)) print(angle2(P, A, B))
出力:
-45.00000000000001
6. 点がポリゴンの領域の内側にある点かどうかを判定する関数を作成しテストする
準備が整いましたので、点がポリゴンの領域の内側にある点かどうかを判定する関数を作成し、機能するかテストします。
点(x, y)が1.のアルゴリズムを使用しポリゴンの領域内の点であるかどうかを判定するis_inside_polygon関数を実装します。以下のようになります。
この関数は[(x1, y1), (x2, y2),,,]のような点(x, y)の集合pointsと点p(x, y)を受け取り、pointsからなるポリゴンの領域内に点pが含まれるか判定します。
コード:
import random import matplotlib.pyplot as plt import numpy as np from numpy import inner from numpy.linalg import norm def angle2(P, A, B): ''' args ---- P: np.array((x, y)) A: np.array((x, y)) B: np.array((x, y)) returns ------- theta: float (dgree) ''' # 全体をシフトさせる PA = A - P PB = B - P # cosθを求める costheta = inner(PA, PB)/norm(PA)/norm(PB) # 角度を求める(度) theta = np.arccos(costheta)/np.pi * 180 # 時計回りか、反時計回りか判定し、thetaをプラスかマイナスにして返す if PA[0]*PB[1] - PA[1]*PB[0] > 0: # 条件(i)'に該当するときはなにもしない(ここではthetaはプラス) pass else: # 条件(ii)'に該当するときはthetaをマイナスにする theta *= -1.0 return theta def is_inside_polygon(points, P): ''' args ---- points: list of np.array((x, y)) P: np.array((x, y)) returns ------- theta: x (degree) ''' points = [np.array(p) for p in points] P = np.array(P) angle_total = 0.0 for A, B in zip(points, points[1:] + points[0:1]): # Notice! "points[1:] + points[0:1]" means something like "points[1:].extend(points[0:1])" angle_total += angle2(P, A, B) if np.isclose(abs(angle_total), 360): # 360°の場合はポリゴンの内側なのでTrueを返す return True else: # ポリゴンの外側 return False # 三角形 #points = ((0,0), (1,0), (0, 1)) # リボン? #points = ((1.0, 0.0), (0.5, 0.5), (1.0, 1.0), (-1.0, 1.0), (-0.5, 0.5), (-1.0, 0.0)) # 十字 points = ((1.0, -3.0), (1.0, -1.0), (3.0, -1.0), (3.0, 1.0), (1.0, 1.0), (1.0, 3.0), (-1.0, 3.0), (-1.0, 1.0), (-3.0, 1.0), (-3.0, -1.0), (-1.0, -1.0), (-1.0, -3.0)) points = [(i[0]/3, i[1]/3) for i in points] N = 5000 for _ in range(N): # 点を打つ場所をランダムに決める(ただし、-1<x<1, -1<y<1とする) x = (random.random() - 0.5) * 2 y = (random.random() - 0.5) * 2 p = (x, y) # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ if is_inside_polygon(points, p): plt.plot(x, y, color='k', marker='.') # グラフを表示する plt.gca().set_xlim(-2,2) plt.gca().set_ylim(-2,2) plt.gca().set_aspect('equal') plt.show()
7. 出力結果
上から順に三角形、リボン、十字形になります。ランダムに点を打っているのですが、pointsで定義された多角形の内側にある点だけ残ることがわかります。
※形を変えるには、多角形の定義(#points... の箇所)を変えてください。
8. まとめ
ポリゴンに収まる点の集合を作成できました。
以上です。
Pythonで2次元データ作成ー三角形に収まる点の集合
前回の続きです。今回は、前回に引き続き、Pythonで三角形に収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それが三角形の内側であれば点を残す、三角形の外側であれば点を消すを繰り返すだけです。「三角形の内側か?」を判定する関数を作成し、その関数に点の位置を次々と問い合わせればよいでしょう。ここでは三角形は直角二等辺三角形ということにします。
前回の記事:
keimina.hatenablog.jp
目次
1. 三角形の領域の定義
実装する前に使用する三角形の形を明確にします。以下のように使用する三角形の領域を定義します。
三角形の領域の定義:
x-y座標上の点(x, y)で「x > 0 かつ y > 0 かつ y < -x + 1」を満たす点の集合
この条件を満たす点(x, y)は3つの点(0, 0)と(1, 0)と(0, 1)を結ぶとできる三角形の領域を構成します。なぜそうなるかについては、以下の①②③を参考にしてください。
①「x > 0 かつ y > 0」 の領域は x-y座標の右上の領域すべてを表します。
②「y < -x + 1』の領域 は 直線 y = -x + 1 より下側の領域をすべてを表します。
③「x > 0 かつ y > 0 かつ y < -x + 1」 の領域は上記①②の領域が重なった領域であり、三角形になります。
①②③を図で表すと以下になります。
2. 三角形の領域内の点であるかどうかを判定する関数を作成する
点(x, y)が上記1.で定義したの三角形の領域内の点であるかどうかを判定する関数を実装します。以下のようになります。
def is_inside_triangle(x, y): if x >= 0 and \ y >= 0 and \ y <= -x + 1: # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」 # すなわち点(x,y)が三角形の内側ならTrue return True else: # 三角形の外側ならFalse return False
is_inside_triangle関数は引数で与えられた点(x, y)の情報が1.で定義した三角形の内側にあるかどうかを判定します。三角形の内側であればTrue、それ以外はFalseを戻します。
3. コード1(ライブラリを極力使わない場合)
2.で作成した関数を使用することで、三角形に収まる点の集合を求めることができます。ライブラリを極力つかわないで全体のコードを書くとこんな感じになるかと思います。前回の記事のコードにある、is_inside_square関数 を 今回作成したis_inside_triangle関数 に置き換えただけです。
※matplotlibは可視化のため使用
import random import matplotlib.pyplot as plt def is_inside_triangle(x, y): if x >= 0 and \ y >= 0 and \ y <= -x + 1: # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」 # すなわち点(x,y)が三角形の内側ならTrue return True else: # 三角形の外側ならFalse return False N = 1000 for _ in range(N): # 点を打つ場所をランダムに決める x = random.random() y = random.random() # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ if is_inside_triangle(x, y): plt.plot(x, y, color='k', marker='.') # グラフを表示する plt.gca().set_xlim(-2,2) plt.gca().set_ylim(-2,2) plt.gca().set_aspect('equal') plt.show()
4. コード2(Numupyを使用して書いた場合)
Numpyを使用して書くと以下のようになると思います。これも、前回の記事のコードにある、is_inside_square関数 を 今回作成したis_inside_triangle関数 に置き換えただけです。
import numpy as np import matplotlib.pyplot as plt def is_inside_triangle(x, y): if x >= 0 and \ y >= 0 and \ y <= -x + 1: # 「xが0以上の領域」かつ「yが0以上の領域」かつ「y=-x+1の直線より下の領域」 # すなわち点(x,y)が三角形の内側ならTrue return True else: # 三角形の外側ならFalse return False # 関数をnumpy化する is_inside_triangle = np.vectorize(is_inside_triangle) N = 1000 # ランダムに点を打つ場所を決める x = np.random.random(N) y = np.random.random(N) # 三角形の内側にある点だけを抜き出す conditions = is_inside_triangle(x, y) x = x[conditions] y = y[conditions] # グラフに点を打つ plt.plot(x, y, color='k', marker='.') # グラフを表示する plt.gca().set_xlim(-2,2) plt.gca().set_ylim(-2,2) plt.gca().set_aspect('equal') plt.show()
Pythonで2次元データ作成ー四角形に収まる点の集合
Pythonで四角形に収まる点の集合を作成します。アルゴリズムは簡単で、ランダムに点を打ち、それが四角形の内側であれば点を残す、四角形の外側であれば点を消すを繰り返すだけです。「四角形の内側か?」を判定する関数を作成し、その関数に様々な点を与え、戻り値がTrueのときの点を集めればよいということになります。ここでは四角形は正方形ということにします。一つ目のソースコードは以下になります。
1. ライブラリを極力つかわないで書いた場合
ライブラリを極力つかわないで書くとこんな感じになるかと思います。(matplotlibは可視化のため使用)
コード1:
import random import matplotlib.pyplot as plt def is_inside_square(x, y): if 0 <= x <= 1 and 0 <= y <= 1: # 「xが0以上1以下の領域」かつ「yが0以上1以下の領域」すなわち点(x,y)が四角形の内側ならTrue return True else: # 四角形の外側ならFalse return False N = 1000 for _ in range(N): # 点を打つ場所をランダムに決める x = random.random() y = random.random() # 決めた点(x, y)が三角形の内側にあればグラフに点を打つ if is_inside_square(x, y): plt.plot(x, y, color='k', marker='.') # グラフに点を打つ # グラフを表示する plt.gca().set_xlim(-2,2) plt.gca().set_ylim(-2,2) plt.gca().set_aspect('equal') plt.show()
random.random()を使用すると0から1までの少数をランダムに選んでくれます。これを利用することで、点をランダムに打つことができます。基本的には、ランダムに打った点の位置がなんらかの形状の内側にあるか、外側にあるか判定できれば、その形状を点で表すことができると思います。
2. 行数を短く書いた場合
上記コードに加え、以下の二つ目のコードも載せておきます。結果は同じです。random.random()の戻り値が[0, 1)の範囲内の値なので、四角形の内側であるかの判定は、実はこの問題では必要ないのです。またnumpyを使いfor文をなくしています。
コード2:
import numpy as np import matplotlib.pyplot as plt N = 1000 # ランダムに点を打つ場所を決める x = np.random.random(N) y = np.random.random(N) # グラフに点を打つ plt.plot(x, y, color='k', marker='.') # グラフを表示する plt.gca().set_xlim(-2,2) plt.gca().set_ylim(-2,2) plt.gca().set_aspect('equal') plt.show()
3. 最終的なコード
上記のコードをみると短くはなりましたが、四角形の内側かを判定するis_inside_square関数があったほうがよいと思いました。なぜなら、関数があったほうが何をしたいのかがわかるのと、四角形ではなく三角形や他の形で同様のことをしたいときに役立つと思ったためです。
そのため最終的なコードは以下のようになりました。
コード3:
import numpy as np import matplotlib.pyplot as plt def is_inside_square(x, y): if 0 <= x <= 1 and 0 <= y <= 1: # 「xが0以上1以下の領域」かつ「yが0以上1以下の領域」すなわち点(x,y)が四角形の内側ならTrue return True else: # 四角形の外側ならFalse return False # 関数をnumpy化する is_inside_square = np.vectorize(is_inside_square) N = 1000 # ランダムに点を打つ場所を決める x = np.random.random(N) y = np.random.random(N) # 四角形の内側にある点だけを抜き出す conditions = is_inside_square(x, y) x = x[conditions] y = y[conditions] # グラフに点を打つ plt.plot(x, y, color='k', marker='.') # グラフを表示する plt.gca().set_xlim(-2,2) plt.gca().set_ylim(-2,2) plt.gca().set_aspect('equal') plt.show()
関数をnumpy化するのにnp.vectorizeを使用しています。
これでPythonで四角形に収まる点の集合を作成できました。同様の方法で三角形や円の内側に収まる点の集合も作成できます。
以上です。おやすみなさい。
numpyの2次元配列の非ゼロ領域を囲む四角形の情報を取得する方法について理解する
本記事は、Stackoverflowにに記載されているnumpyの2次元配列の非ゼロ領域を囲む四角形の情報を取得する方法について理解を深めるための記事です。
2. numpyの2次元配列の非ゼロ領域を囲む四角形の情報を取得する方法
以下のURLを参照してください。関数自体は非常にコンパクトで5行くらいで記載されています。本記事では、この5行くらいのコードについて理解を深めるために解説を行います。
上記によると、2次元配列の非ゼロ領域を囲む四角形の情報を取得するコードは以下になります。
def bbox2(img): rows = np.any(img, axis=1) cols = np.any(img, axis=0) rmin, rmax = np.where(rows)[0][[0, -1]] cmin, cmax = np.where(cols)[0][[0, -1]] return rmin, rmax, cmin, cmax
3. 解説
なぜ、上記のコードで非ゼロ領域を囲む四角形の情報が得られるかを解説します。
以下の図1にアルゴリズムの動作をグラフで表しました。np.anyを使って行と列を圧縮しています。np.anyを使うと非ゼロをみつけたら計算が打ち切られるため限界まで高速化できていると思います。行方向に圧縮したものは行のサマリ(要約)に、列方向に圧縮したものは列のサマリであると解釈するとわかりやすいかもしれません。そして、その行のサマリの左端と右端、列のサマリの上端と下端の位置が我々が欲する四角形の情報です。それらをnp.whereを使用して求めています。ここではnp.whereは引数で与えられた1次元配列のTrueのインデックス番号を取得するために使用しています。np.whereの戻り値として2次元配列が戻ってきますが、インデックス番号0しか使わないため
np.where(x)[0]
のように0を指定しています。また、
np.where(x)[0][0] np.where(x)[0][-1]
のようにインデックスに0と-1を指定しているのは左端と右端、上端と下端を得るためです。ちなみに、ここで使用しているnp.whereはnp.nonzeroと同じです*1。最終的に図1の右下にある式が得られます。
図1の右下にある式とStackoverflowにある式とでは見た目が違うので(実際の処理内容は同じですが)式を変形して同じになるか確認します。以下の図2のように、どんどんワンライナーでコードの行数を圧縮します。すると最終的にStackoverflowの回答と同じような形の式が得られます。
4. 注意点
np.whereで得ようとしている配列(2.に記載した関数bbox2の引数img)がすべて0やFalseのときインデックスエラーになるはずですので注意が必要かと思います。
以上です。おやすみなさい。