●基本的な線形代数
- 線形代数関連の関数はモジュール numpy.linalg に用意されている
- ここでは基本的な使い方を簡単に説明する
- 詳細はリファレンスマニュアル Linear algebra を参照
- 関数 norm(V) はベクトル V のノルム (大きさ) を求める
- 行列 \(A\) のトレース (trace, 跡) は numpy.trace(A)
- トレースは行列の対角線上にある要素の和のこと (\(a_{11} + a_{22} + \cdots + a_{nn}\))
- 行列 \(A\) の逆行列 \(A^{-1}\) は inv(A) で、行列式は det(A) で求めることができる
- matrix の場合、逆行列は A.I で求めることができる
- 行列 \(A\) のランク (rank) は matrix_rank(A)
>>> a = np.array([1.0, 1.0])
>>> b = np.array([1.0, 1.0, 1.0])
>>> np.linalg.norm(b)
1.7320508075688772
>>> c = np.array([[2, 3], [1, 2]])
>>> c
array([[2, 3],
[1, 2]])
>>> np.linalg.inv(c)
array([[ 2., -3.],
[-1., 2.]])
>>> c @ np.linalg.inv(c)
array([[1., 0.],
[0., 1.]])
>>> np.linalg.det(c)
1.0
>>> np.linalg.matrix_rank(c)
2
>>> np.trace(c)
4
- 連立一次方程式は solve() で解くことができる
- 鶴亀算
- 鶴と亀、合わせて 100 匹いる。足の合計が 272 本のとき、鶴と亀はそれぞれ何匹ずついるか
- 鶴と亀とトンボが合わせて 10 匹いる。足の合計が 38 本で羽の合計が 14 枚であるとき、鶴と亀とトンボはそれぞれ何匹ずついるか。(トンボの足は 6 本で羽は 4 枚)
- 問題 1 は次の連立方程式を解けば求めることができる
x + y = 100 2x + 4y = 272
>>> a = np.array([[1, 1], [2, 4]]) >>> b = np.array([100, 272]) >>> np.linalg.solve(a, b) array([64., 36.]) >>> np.linalg.inv(a) @ b array([64., 36.])
x + y + z = 10 2x + 4y + 6z = 38 2x + 4z = 14
>>> a = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]]) >>> b = np.array([10, 38, 14]) >>> np.linalg.solve(a, b) array([3., 5., 2.]) >>> np.linalg.inv(a) @ b array([3., 5., 2.])
- a, b がベクトルの場合、内積と同じ
- a, b が行列の場合、a @ b.T と同じ
- a が多次元配列、b がベクトルの場合、a @ b と同じ
- a, b のどちらかがスカラー (数値) の場合、a * b と同じ
- 引数が [a1, a2, ..., an], [b1, b2, ..., bn] とすると [[a1b1, a1b2, ... a1bm], ..., [anb1, anb2, ... anbm]] を求める
- 引数 a, b が多次元配列の場合、平坦化してから計算する
>>> a = np.array([1, 2, 3])
>>> b = np.array([4, 5, 6])
>>> np.inner(a, b)
32
>>> a @ b
32
>>> c = np.array([[1, 2], [3, 4]])
>>> d = np.array([[5, 6], [7, 8]])
>>> np.inner(c, d)
array([[17, 23],
[39, 53]])
>>> c @ d.T
array([[17, 23],
[39, 53]])
>>> c @ d
array([[19, 22],
[43, 50]])
>>> e = np.arange(9).reshape((3, 3))
>>> e
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> np.inner(e, a)
array([ 8, 26, 44])
>>> e @ a
array([ 8, 26, 44])
>>> np.inner(c, 10)
array([[10, 20],
[30, 40]])
>>> np.inner(10, c)
array([[10, 20],
[30, 40]])
>>> np.outer(a, b)
array([[ 4, 5, 6],
[ 8, 10, 12],
[12, 15, 18]])
>>> np.outer(c, d)
array([[ 5, 6, 7, 8],
[10, 12, 14, 16],
[15, 18, 21, 24],
[20, 24, 28, 32]])
>>> np.linalg.matrix_power(np.array([[1, 0], [0, 1]]), 2)
array([[1, 0],
[0, 1]])
>>> np.linalg.matrix_power(np.array([[1, 0], [0, 1]]), 10)
array([[1, 0],
[0, 1]])
>>> np.linalg.matrix_power(np.array([[1, 1], [1, 0]]), 2)
array([[2, 1],
[1, 1]])
>>> np.linalg.matrix_power(np.array([[1, 1], [1, 0]]), 3)
array([[3, 2],
[2, 1]])
>>> np.linalg.matrix_power(np.array([[1, 1], [1, 0]]), 4)
array([[5, 3],
[3, 2]])
>>> np.linalg.matrix_power(np.array([[1, 1], [1, 0]]), 10)
array([[89, 55],
[55, 34]])
>>> np.linalg.matrix_power(np.array([[1, 1], [1, 0]]), 40)
array([[165580141, 102334155],
[102334155, 63245986]])
●配列の操作
- 配列の連結
- concatenate((a1, a2, ...), axis=0), 指定した軸 (axis) で配列を連結する
- たとえば、shape が (x1, y1, z1) と (x2, y2, z2) の配列を連結する場合
- axis = 0: (x1 + x2, y1, z1), 条件 y1 == y2 and z1 == z2
- axis = 1: (x1, y1 + y2, z1), 条件 x1 == x2 and z1 == z2
- axis = 0: (x1, y1, z1 + z2), 条件 x1 == x2 and y1 == y2
- 条件を満たさない場合はエラー
- hstack((a1, a2, ...)), 軸 0 で配列を連結
- vstack((a1, a2, ...)), 軸 1 で配列を連結
- dstack((a1, a2, ...)), 軸 2 で配列を連結
- ベクトルを vstack() で連結すると行列になる
- ベクトルまたは行列を dstack() で連結すると三次元配列になる
>>> a = np.arange(9).reshape((3, 3))
>>> a
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> b = np.arange(10, 19).reshape((3, 3))
>>> b
array([[10, 11, 12],
[13, 14, 15],
[16, 17, 18]])
>>> np.concatenate((a, b), axis=0)
array([[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[10, 11, 12],
[13, 14, 15],
[16, 17, 18]])
>>> np.concatenate((a, b), axis=1)
array([[ 0, 1, 2, 10, 11, 12],
[ 3, 4, 5, 13, 14, 15],
[ 6, 7, 8, 16, 17, 18]])
>>> np.hstack((a, b))
array([[ 0, 1, 2, 10, 11, 12],
[ 3, 4, 5, 13, 14, 15],
[ 6, 7, 8, 16, 17, 18]])
>>> np.vstack((a, b))
array([[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[10, 11, 12],
[13, 14, 15],
[16, 17, 18]])
>>> np.dstack((a, b))
array([[[ 0, 10],
[ 1, 11],
[ 2, 12]],
[[ 3, 13],
[ 4, 14],
[ 5, 15]],
[[ 6, 16],
[ 7, 17],
[ 8, 18]]])
>>> c = np.arange(1, 6)
>>> d = np.arange(6, 11)
>>> c
array([1, 2, 3, 4, 5])
>>> d
array([6, 7, 8, 9, 10])
>>> np.hstack((c, d))
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
>>> np.vstack((c, d))
array([[ 1, 2, 3, 4, 5],
[ 6, 7, 8, 9, 10]])
>>> np.dstack((c, d))
array([[[ 1, 6],
[ 2, 7],
[ 3, 8],
[ 4, 9],
[ 5, 10]]])
- 配列の分割
- split(a, indices_or_sections, axis=0), 指定した軸で配列を分割
- hsplit(a, indices_or_sections), 軸 0 で配列を分割
- vsplit(a, indices_or_sections), 軸 1 で配列を分割
- dsplit(a, indices_or_sections), 軸 2 で配列を分割
- 引数 indices_or_sections が整数 n の場合、配列を n 個に分割する
- たとえば、shape が (x, y, z) で axis=0 の場合、(x / n, y, z) の配列が n 個できる
- split() の場合、x / n が割り切れないとエラー
- array_split() はなるべく大きさが等しくなるように分割する
- 位置を指定して配列を分割する場合、位置を格納したリストを渡す
>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> np.split(a, 5)
[array([0, 1]), array([2, 3]), array([4, 5]), array([6, 7]), array([8, 9])]
>>> np.array_split(a, 3)
[array([0, 1, 2, 3]), array([4, 5, 6]), array([7, 8, 9])]
>>> np.split(a, [4])
[array([0, 1, 2, 3]), array([4, 5, 6, 7, 8, 9])]
>>> np.split(a, [2, 7])
[array([0, 1]), array([2, 3, 4, 5, 6]), array([7, 8, 9])]
>>> b = np.arange(1, 9).reshape(2, 4)
>>> b
array([[1, 2, 3, 4],
[5, 6, 7, 8]])
>>> np.hsplit(b, 2)
[array([[1, 2],
[5, 6]]), array([[3, 4],
[7, 8]])]
>>> np.vsplit(b, 2)
[array([[1, 2, 3, 4]]), array([[5, 6, 7, 8]])]
>>> c = np.arange(1, 9).reshape((2, 2, 2))
>>> c
array([[[1, 2],
[3, 4]],
[[5, 6],
[7, 8]]])
>>> np.dsplit(c, 2)
[array([[[1],
[3]],
[[5],
[7]]]), array([[[2],
[4]],
[[6],
[8]]])]
- 配列の削除は delete(), 挿入は insert() を使う
delete(a, obj, axis=None) insert(a, obj, values, axis=None)
>>> a = np.arange(10, 26).reshape((4, 4))
>>> a
array([[10, 11, 12, 13],
[14, 15, 16, 17],
[18, 19, 20, 21],
[22, 23, 24, 25]])
>>> np.delete(a, 1, 0)
array([[10, 11, 12, 13],
[18, 19, 20, 21],
[22, 23, 24, 25]])
>>> np.delete(a, 1, 1)
array([[10, 12, 13],
[14, 16, 17],
[18, 20, 21],
[22, 24, 25]])
>>> np.delete(a, [0, 2], 0)
array([[14, 15, 16, 17],
[22, 23, 24, 25]])
>>> np.delete(a, [0, 2, 4, 6], None)
array([11, 13, 15, 17, 18, 19, 20, 21, 22, 23, 24, 25])
>>> np.insert(a, 0, 99, 0)
array([[99, 99, 99, 99],
[10, 11, 12, 13],
[14, 15, 16, 17],
[18, 19, 20, 21],
[22, 23, 24, 25]])
>>> np.insert(a, 0, 99, 1)
array([[99, 10, 11, 12, 13],
[99, 14, 15, 16, 17],
[99, 18, 19, 20, 21],
[99, 22, 23, 24, 25]])
>>> np.insert(a, 0, 99, None)
array([99, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25])
>>> np.insert(a, 2, [1, 2, 3, 4], 0)
array([[10, 11, 12, 13],
[14, 15, 16, 17],
[ 1, 2, 3, 4],
[18, 19, 20, 21],
[22, 23, 24, 25]])
>>> np.insert(a, 2, [1, 2, 3, 4], 1)
array([[10, 11, 1, 12, 13],
[14, 15, 2, 16, 17],
[18, 19, 3, 20, 21],
[22, 23, 4, 24, 25]])
>>> np.insert(a, 2, [1, 2, 3, 4], None)
array([10, 11, 1, 2, 3, 4, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23, 24, 25])
- swapaxes(a, n, m) は配列 a の軸 n と m を交換する
>>> a = np.arange(8).reshape((2, 4))
>>> a
array([[0, 1, 2, 3],
[4, 5, 6, 7]])
>>> np.swapaxes(a, 0, 1)
array([[0, 4],
[1, 5],
[2, 6],
[3, 7]])
>>> b = np.arange(8).reshape((2,2,2))
>>> b
array([[[0, 1],
[2, 3]],
[[4, 5],
[6, 7]]])
>>> np.swapaxes(b, 1, 0)
array([[[0, 1],
[4, 5]],
[[2, 3],
[6, 7]]])
>>> np.swapaxes(b, 0, 2)
array([[[0, 4],
[2, 6]],
[[1, 5],
[3, 7]]])
- 関数 numpy.sort(A, axis=-1) は配列 A をソートした新しい配列を返す
- メソッド A.sort(axis=-1) は配列 A を in-place でソートする (破壊的)
- numpy.sort() の場合、axis に None を指定すると、配列を平坦化してソートする
- axis に -1 を指定すると、末尾の軸がソート対象になる
- 二次元であれば axis = 1, 三次元であれば axis = 2 と同じ
>>> a = np.arange(9, 0, -1).reshape((3, 3))
>>> a
array([[9, 8, 7],
[6, 5, 4],
[3, 2, 1]])
>>> np.sort(a)
array([[7, 8, 9],
[4, 5, 6],
[1, 2, 3]])
>>> a
array([[9, 8, 7],
[6, 5, 4],
[3, 2, 1]])
>>> np.sort(a, axis=0)
array([[3, 2, 1],
[6, 5, 4],
[9, 8, 7]])
>>> np.sort(a, axis=None)
array([1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a.sort(axis=0)
>>> a
array([[3, 2, 1],
[6, 5, 4],
[9, 8, 7]])
>>> a.sort(axis=1)
>>> a
array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
●ブロードキャスト
- NumPy では、サイズや形状 (shape) が異なる配列でも条件によっては演算できることがある
- これを「ブロードキャスト (Broadcasting)」と呼ぶ
- ブロードキャストのルール
- 次元数 (ndim) が異なる場合、shape の先頭に 1 を追加することで調整する
- たとえば、(4,) と (1, 4) であれば、(4,) を (1, 4) に拡張する
- 各次元の要素数が、最大値と等しいか 1 であること
- たとえば、(x, y) と演算できるのは (1, y), (x, 1), (x, y)
- 1 のルールを適用すれば、(y, ) も演算できる
- 出力される配列の次元の要素数は各次元の最大値となる
- たとえば、(x, 1) と (1, y) の演算結果は (x, y) になる
- 次元の要素数が 1 の場合、その軸で得られる値が繰り返し使用される
- たとえば、A が (1, y), B が (x, y) の場合、A[0,:] の値を x 回繰り返して (x, y) の行列を作る
- A が (x, 1), B が (x, y) であれば、A[:, 0] の値を y 回繰り返して (x, y) の行列を作る
- A が (x, 1), B が (1, y) であれば、A[:, 0] を y 回、B[0, :] を x 回繰り返して行列を作る
- つまり、[1, 2] + [[3], [4]] => [[1, 2], [1, 2]] + [[3, 3], [4, 4]] = [[4, 5], [5, 6]]
>>> a = np.array([[1, 2, 3, 4]])
>>> a
array([[1, 2, 3, 4]])
>>> a.shape
(1, 4)
>>> b = np.array([5, 6, 7, 8])
>>> b
array([5, 6, 7, 8])
>>> b.shape
(4,)
>>> a + b
array([[ 6, 8, 10, 12]])
>>> b + a
array([[ 6, 8, 10, 12]])
>>> c = np.arange(16).reshape((4, 4))
>>> c
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
>>> c + a
array([[ 1, 3, 5, 7],
[ 5, 7, 9, 11],
[ 9, 11, 13, 15],
[13, 15, 17, 19]])
>>> c + b
array([[ 5, 7, 9, 11],
[ 9, 11, 13, 15],
[13, 15, 17, 19],
[17, 19, 21, 23]])
>>> a.T
array([[1],
[2],
[3],
[4]])
>>> a + a.T
array([[2, 3, 4, 5],
[3, 4, 5, 6],
[4, 5, 6, 7],
[5, 6, 7, 8]])
>>> c + a.T
array([[ 1, 2, 3, 4],
[ 6, 7, 8, 9],
[11, 12, 13, 14],
[16, 17, 18, 19]])
●ファイル入出力
- 配列のセーブとロードは次の関数で行うことができる
- numpy.savetxt(filename, A), 配列 A を ASCII 形式でファイル filename に保存
- numpy.loadtxt(filename), ファイル filename から ASCII 形式のデータを配列に読み込む
- numpy.save(filename, A), 配列 A をバイナリ形式でファイル filename に保存
- numpy.load(filename), ファイル filename からバイナリ形式のデータを配列に読み込む
- save() の場合、拡張子を省略すると .npy になる
- load() の場合、拡張子は省略できない
>>> a = np.arange(16).reshape((4, 4))
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
>>> np.savetxt('a.txt', a)
>>> np.loadtxt('a.txt')
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]])
0.000000000000000000e+00 1.000000000000000000e+00 2.000000000000000000e+00 3.000000000000000000e+00
4.000000000000000000e+00 5.000000000000000000e+00 6.000000000000000000e+00 7.000000000000000000e+00
8.000000000000000000e+00 9.000000000000000000e+00 1.000000000000000000e+01 1.100000000000000000e+01
1.200000000000000000e+01 1.300000000000000000e+01 1.400000000000000000e+01 1.500000000000000000e+01
a.txt の内容
>>> np.save('a', a)
>>> b = np.load('a.npy')
>>> b
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
>>> b.dtype
dtype('int64')
>>> a.astype(np.float_)
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]])
>>> np.save('b', a.astype(np.float_))
>>> b = np.load('b.npy')
>>> b
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]])
>>> b.dtype
dtype('float64')
- 複数の配列を一つのファイルに保存するときは savez(filename, A, B, ...) を使う
- 拡張子を省略すると .npz になる
- npz ファイルの読み込みは load() で行う
- この場合、load() は複数の配列を格納したデータ構造 (NpzFile) を返す
- 最初の配列はキー 'arr_0' で, 次の配列は 'arr_1', n 番目の配列は 'arr_n' でアクセスできる
- 配列を name_a = A のようにキーワード引数で渡すと、キーワード name_a でアクセスできる
- これらのキーは属性 files に格納されている
>>> a = np.arange(9)
>>> b = np.arange(9).reshape((3, 3))
>>> c = np.arange(8).reshape((2, 2, 2))
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> b
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> c
array([[[0, 1],
[2, 3]],
[[4, 5],
[6, 7]]])
>>> np.savez('test', a, b, c)
>>> npz = np.load('test.npz')
>>> type(npz)
<class 'numpy.lib.npyio.NpzFile'>
>>> npz.files
['arr_0', 'arr_1', 'arr_2']
>>> npz['arr_0']
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> npz['arr_1']
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> npz['arr_2']
array([[[0, 1],
[2, 3]],
[[4, 5],
[6, 7]]])
>>> np.savez('test', a = a, b = b, c = c)
>>> npz = np.load('test.npz')
>>> npz.files
['a', 'b', 'c']
>>> npz['a']
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> npz['b']
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> npz['c']
array([[[0, 1],
[2, 3]],
[[4, 5],
[6, 7]]])
●配列の操作 (2)
- 多次元配列の場合、numpy.r_[a, b, ...] の基本動作は vstack() と同じ
- 同じく numpy.c_[a, b, ...] の基本動作は hstack() と同じ
>>> a = np.arange(6).reshape(2, 3)
>>> a
array([[0, 1, 2],
[3, 4, 5]])
>>> b = np.arange(10, 16).reshape(2, 3)
>>> b
array([[10, 11, 12],
[13, 14, 15]])
>>> np.r_[a, b]
array([[ 0, 1, 2],
[ 3, 4, 5],
[10, 11, 12],
[13, 14, 15]])
>>> np.vstack((a, b))
array([[ 0, 1, 2],
[ 3, 4, 5],
[10, 11, 12],
[13, 14, 15]])
>>> np.c_[a, b]
array([[ 0, 1, 2, 10, 11, 12],
[ 3, 4, 5, 13, 14, 15]])
>>> np.hstack((a, b))
array([[ 0, 1, 2, 10, 11, 12],
[ 3, 4, 5, 13, 14, 15]])
- ベクトルの場合、r_[a, b, ...] の動作は hstack() と同じ
- このとき、要素の値、Python のリスト、スライス表記を含めることができる
- c_[a, b, ...] の動作は、ベクトルの shape (n, ) を (1, n) に変更して hstack() したものと同じになる
- c_ も Python のリストやスライス表記を使用できる
>>> np.r_[1, 2, 3, [4, 5, 6], [0] * 3]
array([1, 2, 3, 4, 5, 6, 0, 0, 0])
>>> np.r_[0:10:2]
array([0, 2, 4, 6, 8])
>>> np.c_[[1, 2, 3], [10, 20, 30]]
array([[ 1, 10],
[ 2, 20],
[ 3, 30]])
>>> x = np.array([[1], [2], [3]])
>>> y = np.array([[10], [20], [30]])
>>> x
array([[1],
[2],
[3]])
>>> y
array([[10],
[20],
[30]])
>>> np.c_[x, y]
array([[ 1, 10],
[ 2, 20],
[ 3, 30]])
>>> np.hstack((x, y))
array([[ 1, 10],
[ 2, 20],
[ 3, 30]])
- r_ の場合、第 1 引数に文字列 'a, b, c' を渡すことで軸などの指定もできる
- a, b, c は整数値
- a, 軸 (axis)
- b, 最低の次元数 (ndim)
- c, shape を揃えるときに挿入する 1 の位置
- デフォルトは 0, 0, -1
- -1 は末尾を表す
- r_['-1, 2, 0', ...] は c_ と同じ動作になる
>>> np.r_['0', a, b]
array([[ 0, 1, 2],
[ 3, 4, 5],
[10, 11, 12],
[13, 14, 15]])
>>> np.r_['1', a, b]
array([[ 0, 1, 2, 10, 11, 12],
[ 3, 4, 5, 13, 14, 15]])
>>> np.r_['0, 2', [1, 2, 3], [4, 5, 6]]
array([[1, 2, 3],
[4, 5, 6]])
>>> np.r_['1, 2', [1, 2, 3], [4, 5, 6]]
array([[1, 2, 3, 4, 5, 6]])
>>> np.r_['-1, 2, 0', [1, 2, 3], [4, 5, 6]]
array([[1, 4],
[2, 5],
[3, 6]])
- 配列の添字にインデックス配列を指定すると、その要素を格納した一次元配列を返す
- 多次元配列を返したい場合は numpy.ix_() を使うと便利
- 二次元配列の場合は次のようになる
numpy.ix_([x0, x1, ..., xn], [y0, y1, ..., ym]) => インデックス配列
>>> a = np.arange(25).reshape((5, 5))
>>> a
array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14],
[15, 16, 17, 18, 19],
[20, 21, 22, 23, 24]])
>>> b = np.ix_([1, 3], [0, 2, 4])
>>> b
(array([[1],
[3]]), array([[0, 2, 4]]))
>>> c = a[b]
>>> c
array([[ 5, 7, 9],
[15, 17, 19]])
>>> a[b] = 0
>>> a
array([[ 0, 1, 2, 3, 4],
[ 0, 6, 0, 8, 0],
[10, 11, 12, 13, 14],
[ 0, 16, 0, 18, 0],
[20, 21, 22, 23, 24]])
>>> c
array([[ 5, 7, 9],
[15, 17, 19]])
- 多次元配列のイテレータは、最初の軸 (axis = 0) が対象となる
- 値を順番に取り出したい場合は属性 flat を使うと簡単
- 添字が必要な場合は numpy.ndenumerate(A) を使う
- ndindex() はインデックスを返すイテレータを生成する
- 引数は配列の形状 (shape)
>>> a = np.array([[1, 2, 3],[4, 5, 6]])
>>> a
array([[1, 2, 3],
[4, 5, 6]])
>>> for x in a: print(x)
...
[1 2 3]
[4 5 6]
>>> for x in a:
... for y in x: print(y)
...
1
2
3
4
5
6
>>> for x in a.flat: print(x)
...
1
2
3
4
5
6
>>> for i, x in np.ndenumerate(a): print(i, x)
...
(0, 0) 1
(0, 1) 2
(0, 2) 3
(1, 0) 4
(1, 1) 5
(1, 2) 6
>>> for i in np.ndindex((2, 3)): print(i)
...
(0, 0)
(0, 1)
(0, 2)
(1, 0)
(1, 1)
(1, 2)
>>> for i in np.ndindex((2, 3)): print(i, a[i])
...
(0, 0) 1
(0, 1) 2
(0, 2) 3
(1, 0) 4
(1, 1) 5
(1, 2) 6
●集合関数
- NumPy には配列を集合として扱う関数が用意されている
- unique(xs), xs の重複要素を取り除いた集合 (配列) を返す
- union1d(xs, ys), 和集合
- intersect1d(xs, ys), 積集合
- setdiff1d(xs, ys), 差集合
- setxor1d(xs, ys), 排他的論理和
- in1d(xs, ys), xs の要素が ys に含まれているか (真偽値の配列を返す)
- これらの関数は引数に多次元配列を渡すと、それらを平坦化してから演算を行う
>>> a = np.array([1, 2, 3, 4]) >>> b = np.array([3, 4, 5, 6]) >>> np.union1d(a, b) array([1, 2, 3, 4, 5, 6]) >>> np.intersect1d(a, b) array([3, 4]) >>> np.setdiff1d(a, b) array([1, 2]) >>> np.setxor1d(a, b) array([1, 2, 5, 6]) >>> np.in1d(a, b) array([False, False, True, True]) >>> np.in1d(b, a) array([ True, True, False, False]) >>> np.in1d(a, a) array([ True, True, True, True]) >>> c = np.array([5, 4, 3, 2, 1, 3, 5, 7, 0]) >>> c array([5, 4, 3, 2, 1, 3, 5, 7, 0]) >>> np.unique(c) array([0, 1, 2, 3, 4, 5, 7]) >>> c array([5, 4, 3, 2, 1, 3, 5, 7, 0]) >>> np.unique(np.array([[3, 2, 1], [6, 5, 4], [2, 4, 6]])) array([1, 2, 3, 4, 5, 6])
- 部分集合を判定する関数 issubset() は all() と in1d() を使って簡単に定義できる
>>> def issubset(xs, ys): return np.all(np.in1d(xs, ys)) ... >>> issubset(a, b) False >>> issubset(a, np.union1d(a, b)) True >>> issubset(a, a) True
●数学関数
- NumPy には多数の数学関数が用意されている
- 数学関数の一覧はリファレンスマニュアル Mathematical functions を参照
- ここでは、面白そうな関数の使い方を簡単に説明する
- sum(), cumsum(), 和と累積和
- prod(), cumprod(), 累積和と累積積
- diff(), ediff1d(), 階差 (差分)
- diff() のキーワード引数 n は階差を表す (デフォルトは 1)
- たとえば、n = 2 とすると、要素の差分の数列を求め、さらにその数列の差分を求める
- ediff1d() は引数の配列を平坦化して差分を求め、結果を 1 次元配列に格納して返す
- キーワード引数 to_begin が指定されている場合、その値を配列の先頭に追加する
- キーワード引数 to_end が指定されている場合、その値を配列の先頭に追加する
>>> a = np.arange(1, 11)
>>> a
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
>>> np.sum(a)
55
>>> np.cumsum(a)
array([ 1, 3, 6, 10, 15, 21, 28, 36, 45, 55], dtype=int32)
>>> np.prod(a)
3628800
>>> np.cumprod(a)
array([ 1, 2, 6, 24, 120, 720, 5040,
40320, 362880, 3628800], dtype=int32)
>>> np.diff(a)
array([1, 1, 1, 1, 1, 1, 1, 1, 1])
>>> np.diff(a, n=2)
array([0, 0, 0, 0, 0, 0, 0, 0])
>>> np.diff(a, n=0)
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
>>> np.ediff1d(a)
array([1, 1, 1, 1, 1, 1, 1, 1, 1])
>>> np.ediff1d(a, to_end=2, to_begin=0)
array([0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2])
- gradient(), 勾配
- 第 1 引数は関数 func(x, y, ...) を適用した結果を格納した配列
- 第 2 引数以降に func の引数を格納した配列 xs, ys, ... を渡す
- 第 2 引数が数値の場合、差分 (xs[i+1] と xs[i] の差) を表す
- 第 2 引数が省略された場合、差分は 1 になる
- たとえば、func の引数が 1 つの場合、gradient(fs, xs) の返り値 zs は次のように計算される
- zs[i] = (fs[i + 1] - fs[i - 1]) / (xs[i + 1] - xs[i - 1])
- zs[0] = (fs[1] - fs[0]) / (xs[1] - xs[0])
- zs[-1] = (fs[-1] - fs[-2]) / (xs[-1] - xs[-2])
- それぞれ、数値微分の中央差分、前進差分、後退差分に対応している
- キーワード引数 edge_order に 2 を指定すると、両端の精度がよくなるようだ
- edge_order のデフォルト値は 1
>>> xs = np.arange(11) >>> xs array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) >>> fs = xs ** 2 >>> fs array([ 0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100], dtype=int32) >>> zs1 = np.gradient(fs, xs) >>> zs1 array([ 1., 2., 4., 6., 8., 10., 12., 14., 16., 18., 19.]) >>> np.gradient(zs1, xs) array([1. , 1.5, 2. , 2. , 2. , 2. , 2. , 2. , 2. , 1.5, 1. ]) >>> zs2 = np.gradient(fs, xs, edge_order = 2) >>> zs2 array([ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20.]) >>> np.gradient(zs2, xs) array([2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.])
リスト : gradient() の簡単な使用例 (sample08.py)
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-3, 3, 6001)
y = x ** 3 - x * 3 + 1 # x3 - 3x + 1
z1 = np.gradient(y, x, edge_order = 2) # 3x2 - 3
z2 = np.gradient(z1, x, edge_order = 2) # 6x
plt.plot(x, y)
plt.plot(x, z1)
plt.plot(x, z2)
# plt.savefig('grad.png', dpi=80)
plt.show()
- trapz(), 台形則で数値積分を行う
リスト : 数値積分 (円周率を求める, sample09.py)
import numpy as np
n = 10
for _ in range(5):
# 中点則
w = 1 / n
a = (np.arange(1, n + 1) - 0.5) * w
b = 4.0 / (1.0 + a ** 2)
p1 = np.sum(b) * w
print(n, p1, np.pi - p1)
# 台形則
c = np.linspace(0, 1, n + 1)
d = 4.0 / (1.0 + c ** 2)
p2 = np.trapz(d) * w
print(n, p2, np.pi - p2)
# シンプソン則
p3 = (p1 * 2 + p2) / 3
print(n, p3, np.pi - p3)
n *= 10
$ python3 sample09.py 10 3.1424259850010987 -0.0008333314113055934 10 3.1399259889071587 0.0016666646826344333 10 3.1415926529697855 6.200076008155975e-10 100 3.1416009869231245 -8.333333331389525e-06 100 3.1415759869231277 1.6666666665443586e-05 100 3.141592653589792 1.3322676295501878e-15 1000 3.1415927369231267 -8.33333335670261e-08 1000 3.141592486923127 1.6666666624587378e-07 1000 3.141592653589793 0.0 10000 3.1415926544231265 -8.333334022836425e-10 10000 3.1415926519231268 1.6666663604780751e-09 10000 3.141592653589793 0.0 100000 3.1415926535981273 -8.334222201256125e-12 100000 3.1415926535731273 1.666577986725315e-11 100000 3.1415926535897944 -1.3322676295501878e-15
- 誤差の絶対値は台形則が中点則の 2 倍程度、誤差の符号は逆になる
- 両者を 1 : 2 で割合で重みづけして平均値をとれば、より正確な値を求めることができる
- これを「シンプソン則」という
- シンプソン則で数値積分を行う関数は NumPy には用意されていない
- NumPy で見つからない場合はライブラリ SciPy を探してみるとよい
- SciPy には高度な関数が多数用意されている
- convolve(a, v, mode='full'), 畳み込み積分
- 畳み込み積分については 畳み込み - Wikipedia を参照
- 引数 a, v は配列
- a, v の大きさが n, m とすると convolve() は次の計算結果を配列に格納して返す
- a[i] * v[0] + a[i+1] * v[1] + ... + a[i+m-1] * v[m-1]
- mode が full の場合、a の先頭と v の終端が一致するところから始まる
- そして、a の終端と v の先頭が一致したところで終了する
- mode が same の場合、a と同じ大きさになるように開始位置と終了位置を調整する
- mode が valid の場合、a と v の先頭が一致したところから始まり、a と v の終端が一致したところで終了する
- a と v が重ならない (a の要素がない) 場合は 0 が補填される
>>> a = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0]) >>> np.convolve(a, [1, 1, 1, 1, 1]) array([0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 4, 3, 2, 1, 0, 0, 0, 0, 0]) >>> np.convolve(a, [1, 1, 1, 1, 1], mode='same') array([0, 0, 0, 1, 2, 3, 4, 5, 4, 3, 2, 1, 0, 0, 0]) >>> np.convolve(a, [1, 1, 1, 1, 1], mode='valid') array([0, 1, 2, 3, 4, 5, 4, 3, 2, 1, 0])
- 簡単な例題として、移動平均法のプログラムを示す
- 移動平均法については拙作のページ 統計学の基礎知識 [3] を参照
リスト : 移動平均法 (sample10.py)
import numpy as np
import matplotlib.pyplot as plt
# 東京の年平均気温 (1975 - 2014)
data = np.array([
15.6, 15.0, 15.8, 16.1, 16.9, 15.4, 15.0, 16.0, 15.7, 14.9,
15.7, 15.2, 16.3, 15.4, 16.4, 17.0, 16.4, 16.0, 15.5, 16.9,
16.3, 15.8, 16.7, 16.7, 17.0, 16.9, 16.5, 16.7, 16.0, 17.3,
16.2, 16.4, 17.0, 16.4, 16.7, 16.9, 16.5, 16.3, 17.1, 16.6
])
x = np.arange(1975, 2015)
# 移動平均法
data1 = np.convolve(data, np.ones(5) / 5.0, mode = 'valid')
x1 = np.arange(1977, 2013)
# 散布図
plt.plot(x, data, 'o')
plt.plot(x1, data1, 'o')
# 回帰直線
a, b = np.polyfit(x, data, 1)
a1, b1 = np.polyfit(x1, data1, 1)
plt.plot(x, x * a + b)
plt.plot(x1, x1 * a1 + b1)
# plt.savefig('tokyo1.png', dpi=80)
plt.show()
●多項式
- NumPy には多項式 (polynomial) を扱う関数が用意されている
- poly1d([An, An-1, ..., A1, A0]) => p, 多項式 AnXn + An-1Xn-1 + ... + A1X + A0 を生成する
- 多項式 p は p(x) で計算できる
- 引数 x は配列 (ndarray) や他の多項式でも OK
- polyval(xs, x) は p = poly1d(xs), p(x) と同じ
- 引数 xs は多項式でもよい
- poly(xs), 配列 xs の要素が根となる多項式の係数を求める
- roots(p), 多項式 p の根を求める
- p.r でも根を求めることができる
- polyder(p), 導関数
- polyint(p), 不定積分
- polyadd, polysub, polymul, polydiv, 多項式の加減乗除
- 演算子 +, -, *, / も使用できる
- polyfit(xs, ys, n), n 次式で 2 変数の回帰分析
>>> p1 = np.poly1d([1, 0, -1])
>>> p1
poly1d([ 1, 0, -1])
>>> print(p1)
2
1 x - 1
>>> p1(np.linspace(-2, 2, 21))
array([ 3. , 2.24, 1.56, 0.96, 0.44, 0. , -0.36, -0.64, -0.84,
-0.96, -1. , -0.96, -0.84, -0.64, -0.36, 0. , 0.44, 0.96,
1.56, 2.24, 3. ])
>>> p1.r
array([-1., 1.])
>>> np.poly([1, -1])
array([ 1., 0., -1.])
>>> np.polyder(p1)
poly1d([2, 0])
>>> np.polyint(p1)
poly1d([ 0.33333333, 0. , -1. , 0. ])
>>> p2 = np.poly1d([1, 1])
>>> print(p2)
1 x + 1
>>> p1(p2)
poly1d([1., 2., 0.])
>>> p1 + p2
poly1d([1, 1, 0])
>>> p1 - p2
poly1d([ 1, -1, -2])
>>> p1 * p2
poly1d([ 1, 1, -1, -1])
>>> p1 / p2
(poly1d([ 1., -1.]), poly1d([0.]))
リスト : 多項式の簡単な使用例 (sample11.py)
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-3, 3, 31)
p = np.poly1d([1, 1, -4, 1])
y = p(x) + np.random.randn(31) * 4
plt.plot(x, y, 'o')
p1 = np.poly1d(np.polyfit(x, y, 3))
print(p1)
plt.plot(x, p1(x))
# plt.savefig('poly01.png', dpi=80)
plt.show()
3 2
0.8696 x + 0.5349 x - 3.553 x + 2.93
●meshgrid
- 格子状の座標を生成したい場合、meshgrid() が役に立つ
meshgrid(xs, ys, zs, ...) => [Xs, Ys, Zs, ...]
>>> np.meshgrid([1, 2, 3], [10, 20, 30])
[array([[1, 2, 3],
[1, 2, 3],
[1, 2, 3]]), array([[10, 10, 10],
[20, 20, 20],
[30, 30, 30]])]
>>> a, b = np.meshgrid([1, 2, 3], [10, 20, 30])
>>> a
array([[1, 2, 3],
[1, 2, 3],
[1, 2, 3]])
>>> b
array([[10, 10, 10],
[20, 20, 20],
[30, 30, 30]])
>>> a + b
array([[11, 12, 13],
[21, 22, 23],
[31, 32, 33]])
>>> x, y, z = np.meshgrid([1, 2], [10, 20], [100, 200])
>>> x
array([[[1, 1],
[2, 2]],
[[1, 1],
[2, 2]]])
>>> y
array([[[10, 10],
[10, 10]],
[[20, 20],
[20, 20]]])
>>> z
array([[[100, 200],
[100, 200]],
[[100, 200],
[100, 200]]])
>>> x + y + z
array([[[111, 211],
[112, 212]],
[[121, 221],
[122, 222]]])
- numpy.mgrid[...] でも同様なことができる
- mgrid はスライス表記を使用することができる
>>> x, y = np.mgrid[-1:1.5:0.5, -1:1.5:0.5]
>>> x
array([[-1. , -1. , -1. , -1. , -1. ],
[-0.5, -0.5, -0.5, -0.5, -0.5],
[ 0. , 0. , 0. , 0. , 0. ],
[ 0.5, 0.5, 0.5, 0.5, 0.5],
[ 1. , 1. , 1. , 1. , 1. ]])
>>> y
array([[-1. , -0.5, 0. , 0.5, 1. ],
[-1. , -0.5, 0. , 0.5, 1. ],
[-1. , -0.5, 0. , 0.5, 1. ],
[-1. , -0.5, 0. , 0.5, 1. ],
[-1. , -0.5, 0. , 0.5, 1. ]])
>>> np.sqrt(x**2 + y**2)
array([[1.41421356, 1.11803399, 1. , 1.11803399, 1.41421356],
[1.11803399, 0.70710678, 0.5 , 0.70710678, 1.11803399],
[1. , 0.5 , 0. , 0.5 , 1. ],
[1.11803399, 0.70710678, 0.5 , 0.70710678, 1.11803399],
[1.41421356, 1.11803399, 1. , 1.11803399, 1.41421356]])
- meshgrid() は matplotlib でグラフを描画するとき役に立つ
- matplotlib.contour() を使うと等高線を簡単に描画することができる
リスト : 等高線の表示 (sample12.py)
import numpy as np
import matplotlib.pyplot as plt
x, y = np.mgrid[-1:1.05:0.05, -1:1.05:0.05]
z = np.sqrt(x ** 2 + y ** 2)
c = plt.contour(x, y, z)
c.clabel(fmt='%1.1f', fontsize=14)
plt.gca().set_aspect('equal')
# plt.savefig('mesh.png', dpi=80)
plt.show()
- mpl_toolkits.mplot3d をインポートすると 3D 表示ができる
- ドラッグで視点を変えることもできる
リスト : 等高線の 3D 表示 (sample13.py)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
x, y = np.mgrid[-1:1.05:0.05, -1:1.05:0.05]
z = np.sqrt(x ** 2 + y ** 2)
ax = plt.figure().add_subplot(projection='3d')
ax.contour3D(x, y, z)
# plt.savefig('mesh1.png', dpi=80)
plt.show()
- contour3D() を plot_wireframe() に変えると、ワイヤーフレームで表示される
リスト : ワイヤーフレーム (sample14.py)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
x, y = np.mgrid[-1:1.05:0.1, -1:1.05:0.1]
z = np.sqrt(x ** 2 + y ** 2)
ax = plt.figure().add_subplot(projection='3d')
ax.plot_wireframe(x, y, z)
# plt.savefig('mesh2.png', dpi=80)
plt.show()
- 散布図 (3D) の作成は scatter3D() を使う
- 引数は 1 次元配列に平坦化してから渡す
リスト : 散布図 (3D, sample15.py)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
x, y = np.mgrid[-1:1.05:0.1, -1:1.05:0.1]
z = np.sqrt(x ** 2 + y ** 2)
ax = plt.figure().add_subplot(projection='3d')
ax.scatter3D(x.flat, y.flat, z.flat)
# plt.savefig('mesh3.png', dpi=80)
plt.show()
