Numpy

numpy(「ナンパイ」と読む)は、Pythonで科学技術計算のための中核となるライブラリです。これは、高性能の多次元配列オブジェクトを実装し、高速でメモリ効率のよい処理ツールを提供しています。MATLABに精通している人なら、 Matlabユーザーのための NumPyyという文書から有用な情報が得られるでしょう。</P>

Numpyを使うには、まず次のようにしてnumpyパッケージをimport(取り込む)必要があります:

In [6]:
import numpy as np

配列

numpy配列はndarrayと呼ばれ、この名称はN次元配列オブジェクト(N-dimensional array)に由来しています。ndarrayの要素は、すべて同じ型の値でなければなりません。これがリストとの大きな違いの一つです。

まずndarrayの算術演算を紹介します。ndarrayに対する算術操作は、その配列要素すべてに作用します。

In [7]:
data = np.array([[1,2,3],[-1,-2,-3]])
print(data)
data
[[ 1  2  3]
 [-1 -2 -3]]
Out[7]:
array([[ 1,  2,  3],
       [-1, -2, -3]])
In [8]:
print(data * 10)
print(data / 5.0)
print(data + 0.5)
[[ 10  20  30]
 [-10 -20 -30]]
[[ 0.2  0.4  0.6]
 [-0.2 -0.4 -0.6]]
[[ 1.5  2.5  3.5]
 [-0.5 -1.5 -2.5]]
In [9]:
print(data + 0.1*data)
[[ 1.1  2.2  3.3]
 [-1.1 -2.2 -3.3]]
In [10]:
print(data * data)
[[1 4 9]
 [1 4 9]]
In [11]:
print(1.0/ data)
[[ 1.          0.5         0.33333333]
 [-1.         -0.5        -0.33333333]]

統計用の関数も用意されています。meanは平均、varは分散、stdは標準偏差を計算します。これらの用語については「確率と確率分布」で述べます。

In [12]:
# 平均
print(data[0].mean())  # ベクトルに対し計算。 np.mean(data[0]) でもよい
print(data.mean())   # 配列に対し計算。 np.mean(data) でもよい

# 分散
print(data[0].var()) # ベクトルに対し計算。 np.var(data[0]) でもよい
print(data.var())   # 配列に対し計算。 np.var(data) でもよい

# 標準偏差
print(data[0].std()) # ベクトルに対し計算。 np.std(data[0]) でもよい
print(data.std())   # 配列に対し計算。 np.std(data) でもよい
2.0
0.0
0.666666666667
4.66666666667
0.816496580928
2.16024689947

ndarray配列にはshapedtypeという属性があり、それぞれの配列ごとに固有の値を持ちます。shapeは配列の次元数とそのサイズとのタプルです。dtypeは配列の要素として期待される「型」を示します。

In [13]:
print(data.shape)   # 2次元配列
print(data.dtype)   # 要素は整数

data2 = np.array([1.0, 2.0, 3.0])
print(data2.shape)  #  1次元配列
print(data2.dtype)  #  要素は浮動小数点数
(2, 3)
int32
(3,)
float64

ndarrayの生成

ndarray配列を作るにはいろいろな方法がありますが、一番簡単なのは今見たように、array関数を用いる方法です。 array関数は引数としてリストなど(専門用語を使えば、enumerate型のオブジェクト)を取ります。ただし、要素の型がすべて同じ場合に限るので、浮動小数点数と整数が混じっている場合には浮動小数点数に統一されます。

In [14]:
lst = [1, 2.3, 4.0, 5, 6]  # リスト
arr1 = np.array(lst)
print(arr1, arr1.shape, arr1.dtype)
[ 1.   2.3  4.   5.   6. ] (5,) float64

入れ子になったリストからもndarray配列が作れます。

In [15]:
lst2 = [[0,1,2], [3,4,5]]  
arr2 = np.array(lst2)
print(arr2, arr2.shape, arr2.dtype)
[[0 1 2]
 [3 4 5]] (2, 3) int32
In [16]:
lst3 = [(0,1,2), (3,4,5)]  
arr3 = np.array(lst3)
print(arr3, arr3.shape, arr3.dtype)
[[0 1 2]
 [3 4 5]] (2, 3) int32
In [17]:
lst4 = [(0,1,2), [3,4,5]]
arr4 = np.array(lst4)
print(arr4, arr4.shape, arr4.dtype)
[[0 1 2]
 [3 4 5]] (2, 3) int32

Numpyには、配列を作成するための関数がarray以外にもたくさんあります。主要なものを表にし、例を示します。

関数 説明
array 引数にリスト、タプル、配列など列挙型のデータをとり、ndarrayを生成する。要素の型は推測、もしくは明示的に指定されたもの
asarray arrayと同じくndarrayを生成するが、引数がndarrayの場合は新規には生成しない
arange Python組み込みの関数rangeと同じ動作でndarrayを生成する
ones, ones_like onesは指定されたサイズのndarrayを、指定されたdtypeで生成し、要素をすべて1とする。ones_likeは引数に別の列挙型オブジェクトを取り、それをテンプレートとして要素をすべて1としたndarray配列を生成する。
zeros, zeros_llike 要素を1ではなく0とすること以外は、onesおよびones_likeと同じ
empty, empty_like 要素を初期化せず不定であること以外は、onesおよびones_likeと同じ
eye, identity NxNの単位行列となるndarray配列を生成する。
In [18]:
a = np.zeros((2,2))  # すべての要素が0である配列を作る
print( a )
[[ 0.  0.]
 [ 0.  0.]]
In [19]:
b = np.ones((1,2))   # すべての要素が1である配列を作る
print (b)
[[ 1.  1.]]
In [20]:
c = np.full((2,2), 7) # 定数配列を作る
print (c)
[[ 7.  7.]
 [ 7.  7.]]
C:\Anaconda3\lib\site-packages\numpy\core\numeric.py:294: FutureWarning: in the future, full((2, 2), 7) will return an array of dtype('int32')
  format(shape, fill_value, array(fill_value).dtype), FutureWarning)
In [21]:
d = np.eye(2)        #  2x2 の単位行列を作る
print (d)
[[ 1.  0.]
 [ 0.  1.]]
In [33]:
e = np.random.random((2,2)) # 配列を作り、乱数で要素を埋める
print (e)
[[ 0.44860158  0.2469678 ]
 [ 0.28257495  0.68850618]]

配列の別な作成方法については このドキュメントを参照してください。

なお、乱数生成には numpyのrandomモジュールの方が randomモジュールよりもいろいろな関数が用意されています。

In [25]:
import random
print(len(dir(random)))
import numpy as np
print(len(dir(np.random)))
dir(np.random)
help(np.random.random)
58
74
Help on built-in function random_sample:

random_sample(...) method of mtrand.RandomState instance
    random_sample(size=None)
    
    Return random floats in the half-open interval [0.0, 1.0).
    
    Results are from the "continuous uniform" distribution over the
    stated interval.  To sample :math:`Unif[a, b), b > a` multiply
    the output of `random_sample` by `(b-a)` and add `a`::
    
      (b - a) * random_sample() + a
    
    Parameters
    ----------
    size : int or tuple of ints, optional
        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
        ``m * n * k`` samples are drawn.  Default is None, in which case a
        single value is returned.
    
    Returns
    -------
    out : float or ndarray of floats
        Array of random floats of shape `size` (unless ``size=None``, in which
        case a single float is returned).
    
    Examples
    --------
    >>> np.random.random_sample()
    0.47108547995356098
    >>> type(np.random.random_sample())
    <type 'float'>
    >>> np.random.random_sample((5,))
    array([ 0.30220482,  0.86820401,  0.1654503 ,  0.11659149,  0.54323428])
    
    Three-by-two array of random numbers from [-5, 0):
    
    >>> 5 * np.random.random_sample((3, 2)) - 5
    array([[-3.99149989, -0.52338984],
           [-2.99091858, -0.79479508],
           [-1.23204345, -1.75224494]])

配列のインデックス

Numpyでは、配列の要素を参照する方法がいろいろあります。これは奥が深いテーマです。

まず1次元ndarrayについてインデックス参照を見ていきましょう。これはPythonのリストの場合とよく似ています。

In [33]:
arr = np.arange(10)   # array([0,1,2,3,4,5,6,7,8,9])
print(arr)
arr
[0 1 2 3 4 5 6 7 8 9]
Out[33]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [27]:
print(arr[5])  
5

ndarrayから切り出した一部(スライス)を取り出すことができます。

In [28]:
print(arr[5:8])
[5 6 7]

次はarr[5], arr[6], arr[7]の3つの要素それぞれに12という値を代入するものです。これをブロードキャストと呼びます。

In [29]:
arr[5:8] = 12
print(arr)
[ 0  1  2  3  4 12 12 12  8  9]

Pythonのリストとndarrayの違いのひとつは、ndarrayのスライスは元のデータのコピーではなく元のデータを指している、ということです。つまりndarray配列のスライスに対する変更は元のndarray配列のに対する変更になります。これをリストと比較する形で示しましょう。

In [30]:
# Pythonのリスト
lst = list(range(10))
lst_slice = lst[5:8]
lst_slice[1] = 12
print(lst)
# ndarray
arr = np.array(np.arange(10))
arr_slice = arr[5:8]
arr_slice[1] = 12
print(arr)
arr_slice[:]=100
print(arr)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
[ 0  1  2  3  4  5 12  7  8  9]
[  0   1   2   3   4 100 100 100   8   9]

参考: コピーを作るには 次のようにcopy 関数を用います

arr[5:8].copy()

2次元以上の多次元配列の場合、インデックス参照されるものがスカラー値ではなく配列になります。そのためインデックス参照の方法が増えます。

In [63]:
arr2 = np.array([[0,1,2], [3,4,5], [6,7,8]] )
print(arr2[1])
[3 4 5]

したがって個々の要素を取り出すには、arr2[1][2]のように階層的にアクセする必要があります。 また次に見るように、もう少し簡単に表す方法もあります。

In [64]:
print(arr2[1][2], arr2[1,2], arr2[(1,2)])
(5, 5, 5)

ndarray配列はスライス記法で部分配列を取り出せますが、配列のそれぞれの次元にスライスを指定することができます。</p>

In [34]:
# 次のような階数2で形状(3,4)の配列を作る
# [[ 1  2  3  4]
#  [ 5  6  7  8]
#  [ 9 10 11 12]]
a = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])

# スライスを用いて最初の2行、1列目と2列目の要素からなる部分配列を抜き出す
#  b はそのような形状(2,2)の配列
# [[2 3]
#  [6 7]]
b = a[:2, 1:3]#a=1,2行取り出し、b=2~3列取り出し
print (b)
[[2 3]
 [6 7]]

配列の真ん中の行にあるデータにアクセスする2通リの方法を見ます。整数インデックスとスライスを混在させることで低い階数の配列ができます。一方、スライスだけを使うと、元の配列と同じ階数の配列ができます。

In [25]:
row_r1 = a[1, :]    # 配列aの2行目に対する階数 1 のスライス
row_r2 = a[1:2, :]  # 配列aの2行目に対する階数 2 のスライス
row_r3 = a[[1], :]  # 配列aの2行目に対する階数 2のスライス
print (row_r1, row_r1.shape )
print (row_r2, row_r2.shape)
print (row_r3, row_r3.shape)
(array([5, 6, 7, 8]), (4,))
(array([[5, 6, 7, 8]]), (1, 4))
(array([[5, 6, 7, 8]]), (1, 4))
In [26]:
# 配列の列に対しても同じことが起こります
col_r1 = a[:, 1]
col_r2 = a[:, 1:2]
print (col_r1, col_r1.shape)

print (col_r2, col_r2.shape)
(array([ 2,  6, 10]), (3,))
(array([[ 2],
       [ 6],
       [10]]), (3, 1))

整数インデックス:

スライスを使用してnumpyの配列にインデックスした結果は、常に元の配列の部分配列になります。これと対照的に、整数インデックスは、別の配列のデータを使用して任意の配列を構築することを可能にします。次に例を示します。</p>

In [28]:
a = np.array([[1,2], [3, 4], [5, 6]])

# 整数インデックスの例
# 返される配列は形状(3, )
print (a[[0, 1, 2], [0, 1, 0]] )

# 上の例は次のものと等価
print ( np.array([a[0, 0], a[1, 1], a[2, 0]]) )
[1 4 5]
[1 4 5]
In [27]:
# 整数インデックスを使えば、同じ配列の要素を使い回すことができる
print (a[[0, 0], [1, 1]])

# これは前のものと同じ
print (np.array([a[0, 1], a[0, 1]]))
[2 2]
[2 2]

整数配列のインデックスが便利なのは、行列の各行から1つの要素を選択したり変更したりする場合です。

In [30]:
# 要素を選んで配列を作る
a = np.array([[1,2,3], [4,5,6], [7,8,9], [10, 11, 12]])
print (a)
[[ 1  2  3]
 [ 4  5  6]
 [ 7  8  9]
 [10 11 12]]
In [31]:
# インデックスの配列を作る
b = np.array([0, 2, 0, 1])

# bのインデックスを用いてそれぞれの行からひとつずつ要素を選ぶ
print (a[np.arange(4), b] ) # "[ 1  6  7 11]"と出力
[ 1  6  7 11]
In [32]:
# bのインデックスを用いてそれぞれの行からひとつずつ要素の値を変化させる
a[np.arange(4), b] += 10
print (a)
[[11  2  3]
 [ 4  5 16]
 [17  8  9]
 [10 21 12]]

ブール配列のインデックス

ブール配列インデックスを使用すると、配列の任意の要素を取り出すことができます。このタイプのインデックス参照は、ある条件を満たす配列の要素を選択するために頻繁に使用されます。次に例を示します。</p>

In [7]:
import numpy as np

a = np.array([[1,2], [3, 4], [5, 6]])

bool_idx = (a > 2)  # 2よりも大きな要素を見つける;
                    # これによりaと同じ形状のブール値のNumpy配列が返る
                    #  bool_idx のそれぞれの要素は
                    # 対応するaの要素が2より大きいかどうかの情報を与える

print (bool_idx)
[[False False]
 [ True  True]
 [ True  True]]
In [8]:
# ブール配列インデックスbool_idxを用いて階数1の配列を作る
# bool_idxの要素がTrueである要素に対応するaの要素からなる
print (a[bool_idx])

# 以上の事を一行で完結に書くことができる:
print (a[a > 2])
[3 4 5 6]
[3 4 5 6]

話を簡潔にするため、numpy配列のインデックス参照の詳細は省略しました。もっと知りたい場合は このドキュメントを読みましょう。

データ型

numpy配列はすべて、同じ型の要素から構成されます。Numpyは、配列を構築するために使用できる多数の数値データ型を提供しています。Numpyは配列を作成するときにデータ型を推測しようとしますが、配列を構築する関数は通常、明示的にデータ型を指定するオプションの引数も含みます。次に例を示します。

In [9]:
x = np.array([1, 2])  # Let numpy choose the datatype
y = np.array([1.0, 2.0])  # Let numpy choose the datatype
z = np.array([1, 2], dtype=np.int64)  # Force a particular datatype

print (x.dtype, y.dtype, z.dtype)
int32 float64 int64

numpyのデータ型について詳しく知りたい場合は、 このドキュメントのを読むことをお勧めします。

配列の数学

基本的な数学関数は配列の要素ごとに作用し、暗黙的に拡張された演算子や、numpyモジュールの関数によって利用できます。

In [13]:
x = np.array([[1,2],[3,4]], dtype=np.float64)
y = np.array([[5,6],[7,8]], dtype=np.float64)

# 要素ごとの和;両方とも配列を返す
print (x + y)
print (np.add(x, y))
#help(np.add)
[[  6.   8.]
 [ 10.  12.]]
[[  6.   8.]
 [ 10.  12.]]
In [11]:
# 要素ごとの差; どちらも配列を返す
print (x - y)
print (np.subtract(x, y))
[[-4. -4.]
 [-4. -4.]]
[[-4. -4.]
 [-4. -4.]]
In [12]:
# 要素ごとの積; どちらも配列を返す
print (x * y)
print (np.multiply(x, y))
[[  5.  12.]
 [ 21.  32.]]
[[  5.  12.]
 [ 21.  32.]]
In [14]:
# 要素ごとの商; どちらも配列を返す
# [[ 0.2         0.33333333]
#  [ 0.42857143  0.5       ]]
print (x / y)
print (np.divide(x, y))
[[ 0.2         0.33333333]
 [ 0.42857143  0.5       ]]
[[ 0.2         0.33333333]
 [ 0.42857143  0.5       ]]
In [15]:
print(np.sqrt(2))
1.41421356237
In [38]:
# 要素ごとの平方根; 配列を返す
# [[ 1.          1.41421356]
#  [ 1.73205081  2.        ]]
print (np.sqrt(x))
[[ 1.          1.41421356]
 [ 1.73205081  2.        ]]

*はMATLABとは異なり、行列乗算ではなく、要素単位の乗算です。numpyでは、dot関数を使って、行列の積を計算したり、行列とベクトルの積を求めたり、ベクトルの内積を計算します。dotは、numpyモジュールの関数としても、配列オブジェクトのインスタンス・メソッドとしても利用できます。

In [39]:
x = np.array([[1,2],[3,4]])
y = np.array([[5,6],[7,8]])

v = np.array([9,10])
w = np.array([11, 12])

# ベクトルの内積計算; どちらも 219 を返す
print (v.dot(w))
print (np.dot(v, w))
219
219
In [40]:
# 行列/ ベクトルの積; どちらも階数1 の配列 [29 67] を返す
print (x.dot(v))
print (np.dot(x, v))
[29 67]
[29 67]
In [41]:
# 行列 / 行列の積; どちらも階数 2 の配列を返す
# [[19 22]
#  [43 50]]
print (x.dot(y))
print (np.dot(x, y))
[[19 22]
 [43 50]]
[[19 22]
 [43 50]]

Numpyは、配列の計算を実行するための多くの便利な関数を提供しています。sumは最も有用な関数の一つです。

In [5]:
x = np.array([[1,2],[3,4]])

print (np.sum(x))   # すべての要素の総和を計算; "10"と表示する
print (np.sum(x, axis=0)) # 各列の総和を計算;  "[4 6]"と表示する
print (np.sum(x, axis=1))  # 各行の総和を計算; "[3 7]"と表示する
10
[4 6]
[3 7]

numpyが提供する数学関数の完全なリストは このドキュメントにあります。

配列を使用して数学関数を計算することとは別に、配列内のデータを再形成または操作する必要があることがよくあります。このタイプの操作の最も単純な例は行列を転置することですが、 これには単に配列オブジェクトのT属性を使えばできます:

In [9]:
print (x)
print (x.T)
print(x.transpose())  # transpose関数を使ってもよい
[[1 2]
 [3 4]]
[[1 3]
 [2 4]]
[[1 3]
 [2 4]]
In [42]:
v = np.array([[1,2,3]])
print ( v ) 
print (v.T)
[[1 2 3]]
[[1]
 [2]
 [3]]

もっとも行列クラスを使うと、転置行列だけではなく、逆行列を求めることもできる。

In [20]:
a = np.matrix('1 2 ; 3 4')    # 行列クラスのインスタンスを作成
#  a=np.matrix([[1, 2], [3, 4]])   でも同じ行列が得られる

print(type(a)) # matrixクラス
print(a)

print(a.T) # 転置

inv_a = a.I
print(inv_a) # 逆行列

print(a*inv_a) # 行列の乗算 --- 単位行列になるはず
<class 'numpy.matrixlib.defmatrix.matrix'>
[[1 2]
 [3 4]]
[[1 3]
 [2 4]]
[[-2.   1. ]
 [ 1.5 -0.5]]
[[  1.00000000e+00   1.11022302e-16]
 [  0.00000000e+00   1.00000000e+00]]

ブロードキャスト

ブロードキャストはNumpyがいろいろな形状の配列に対し算術演算を許すような強力なメカニズムです。よくあるケースは、小さな配列と大きな配列があり、小さな配列を何回か大きな配列に対して演算する、という場合です。

例えば、行列のそれぞれの行にある定数ベクトルを加算したいとしましょう。それにはこのようにすればできます:

In [16]:
# ベクトルvを行列xの各行に加算する
# その結果を行列yとして記憶する
x = np.array([[1,2,3], [4,5,6], [7,8,9], [10, 11, 12]])
v = np.array([1, 0, 1])
y = np.empty_like(x)   # xと同じサイズの空の行列を作る

# ループを用いてvを行列xの各行に加える
for i in range(4):
    y[i, :] = x[i, :] + v

print( y )
[[ 2  2  4]
 [ 5  5  7]
 [ 8  8 10]
 [11 11 13]]

これはこれで動くのですが、行列 x がかなり大きい場合、 Python の機能を使ってこのような繰り返しをやるととても遅くなります。ここで、ベクトル v を行列 x の各行に足すというのは、 v のコピーを垂直にたくさん積み上げた行列 vv を作り、それから行列 xと行列vv を要素ごとに足すことに等しいことに注目すると、この操作を次のように実現することができます:

In [17]:
vv = np.tile(v, (4, 1))  # vを4つ積み上げる
print (vv)               # 出力は: "[[1 0 1]
                         #          [1 0 1]
                         #          [1 0 1]
                         #          [1 0 1]]"
[[1 0 1]
 [1 0 1]
 [1 0 1]
 [1 0 1]]
In [18]:
y = x + vv  # xとvvを要素ごとに加算する
print (y)
[[ 2  2  4]
 [ 5  5  7]
 [ 8  8 10]
 [11 11 13]]

Numpyのブロードキャストにより、 このような計算を、実際にはvのコピーを作らずにすませることができます。 このような場合にブロードキャストを考えましょう:

In [5]:
import numpy as np

# ベクトルvを行列xの各行に加え
# その結果を行列yで記憶する
x = np.array([[1,2,3], [4,5,6], [7,8,9], [10, 11, 12]])
v = np.array([1, 0, 1])
y = x + v  # ブロードキャストを使ってvをxのそれぞれの行に足す
print (y)
[[ 2  2  4]
 [ 5  5  7]
 [ 8  8 10]
 [11 11 13]]

y = x + vは、 たとえx(4, 3)という形状であり、 v(3,) という形状であっても、ブロードキャストのおかげで計算できます。 あたかもvが実際に(4, 3)という 形状をもっていて、各行がみな vのコピーであるかのように、和が要素ごとに計算されるのです。

2つの配列に対するブロードキャストは次の規則に従います:

  1. もしも配列が同じ階数でなければ、小さい階数の配列の形状に対し、両者の形状が同じ長さになるまで1を付け加える。
  2. 2つの配列がある次元において「両立可能 compatible」 であるとは、 それらがその次元において同じサイズを持つか、どちらかがその次元のサイズが1の場合をいう。
  3. どの次元においても両立可能ならば、配列はともにブロードキャストが可能である。
  4. ブロードキャストにより、2つの入力配列は、それらのうち、次元ごとに最大の形状をもっているように振る舞う
  5. 一方の配列がサイズ1で、もう一方の配列のサイズが1よりも大きい場合、前者の配列はその次元にそってコピーされたかのように振る舞う。

もしもこの説明で意味がわからなければ、 この解説か、 こちらの解説を読んでみよう。

ブロードキャストができる関数は 汎化関数と呼ばれている。どんなものがそうかは、 このドキュメント を参照すること.

以下に、ブロードキャストの使用例をいくつか示す:

In [19]:
import numpy as np

# ベクトルの外積計算
v = np.array([1,2,3])  # v の形状 (3,)
w = np.array([4,5])    # w の形状 (2,)
# 外積を計算するには、まずvの形状を変えて、
# 形状 (3, 1)の列ベクトルにする; 次にこれをwに対してブロードキャストし、
# 形状 (3, 2)の出力を得る。これが v と wの外積である
# [[ 4  5]
#  [ 8 10]
#  [12 15]]
print (np.reshape(v, (3, 1)) * w)

# ベクトルを行列のそれぞれの「列」に加える
x = np.array([[1,2,3], [4,5,6]])
# x は形状 (2, 3) で v は形状 (3,) であるから、それらをブロードキャストし
# 次のような形状 (2, 3)の行列を得る
# [[2 4 6]
#  [5 7 9]]
print (x + v)

# ベクトルを行列の各「行」に加える
# 行列x は形状 (2, 3) で、 w は形状 (2,)である。
# 行列 x を転置すると形状(3, 2) となり、これをwに対しブロードキャストすることで
# 形状(3, 2)の行列を得る;この結果を転置することで
# 形状 (2, 3)の行列を最終結果として得る。これは行列xの各列に
# ベクトル w を加えたものになっており、結果は:
# [[ 5  6  7]
#  [ 9 10 11]]
print ( (x.T + w).T )
# 別解:  wの形状を形状 (2, 1)のベクトルに変えると
# これを行列 x にブロードキャストして、同じ結果を得る
print (x + np.reshape(w, (2, 1)))

# 行列の定数倍:
# x は形状(2, 3) の行列とする. Numpy ではスカラーは形状 ()のベクトル扱いである
# これらは一緒にブロードキャストされ、形状(2, 3)の次の行列を得る:
# [[ 2  4  6]
#  [ 8 10 12]]
print (x * 2)
[[ 4  5]
 [ 8 10]
 [12 15]]
[[2 4 6]
 [5 7 9]]
[[ 5  6  7]
 [ 9 10 11]]
[[ 5  6  7]
 [ 9 10 11]]
[[ 2  4  6]
 [ 8 10 12]]

ブロードキャストを使うと、一般にコードが簡潔でしかも処理速度が速くなります。 ですから、可能な限りこれを使うよう、努めてください。

Numpyのドキュメント

この短い文書では、知らなければならない重要なことのホンのさわりしか紹介できていないので、 numpyの文献 を読んで、Numpyについて、よりよく知る努力をしてください。

SciPy

Numpyは高性能、多次元の配列を提供し、かつその配列の基本ツールや処理機能を提供している。 SciPyはNumpyをベースに、 その配列を用いた多くの関数や、いろいろな科学技術応用を提供している。

SciPyをよりよく知るには このドキュメントを見てほしい。 ここではSciPyについて役に立つほんの一部について紹介する。

画像処理

SciPyは画像を扱う基本関数を提供している。たとえば、 ディスクから画像を読み込みNumpyの配列に取り込んだり、その配列を画像としてディスクに書き込んだり、 画像のサイズを変えたりする関数がある。ここではそのような関数のいくつかを紹介しよう。

In [20]:
from scipy.misc import imread, imsave, imresize

# JPEG 画像ファイルを読み込み、numpy の配列として記憶する
img = imread('assets/cat.jpg')
print( img.dtype, img.shape ) # "uint8 (400, 248, 3)" と出力する

# 色のチャネルごとにいろいろな定数を掛けてスケーリングすることで
# 画像に色合いをつけることができる。この画像の形状は(400, 248, 3)である
# これに形状(3,) の定数ベクトル[1, 0.95, 0.9] をNumpyのブロードキャスト
# する。つまり赤色(第1チャネル)は変化させないが、緑(第2チャネル)と青
# (第3チャネル)の大きさをそれぞれ0.95倍、および0.9倍する

img_tinted = img * [1, 0.95, 0.9]

# 結果の画像を 300 x 300 ピクセルになるようサイズを変更する
img_tinted = imresize(img_tinted, (300, 300))

# この画像をディスクに書き込む
imsave('assets/cat_tinted.jpg', img_tinted)
uint8 (400, 248, 3)

MATLAB ファイル

scipy.io.loadmatscipy.io.savemat関数を使うと、 MATLABのファイルの読み書きができます。 これについては このドキュメントをご覧ください.

座標間の距離

SciPyは、点(座標)の集合の間の距離を計算するための有用な関数を用意しています、

関数 scipy.spatial.distance.pdistは、 集合中のすべての点と点のペアの間の距離を計算します

In [21]:
import numpy as np
from scipy.spatial.distance import pdist, squareform

# 各行が2D空間における点となっている配列を作る:
# [[0 1]
#  [1 0]
#  [2 0]]
x = np.array([[0, 1], [1, 0], [2, 0]])
print (x)

# xのすべての行の間のユークリッド距離を計算する 
# d[i, j] は x[i, :] と x[j, :] の間の距離となる
# 結果は次:
# [[ 0.          1.41421356  2.23606798]
#  [ 1.41421356  0.          1.        ]
#  [ 2.23606798  1.          0.        ]]
d = squareform(pdist(x, 'euclidean'))
print (d)
[[0 1]
 [1 0]
 [2 0]]
[[ 0.          1.41421356  2.23606798]
 [ 1.41421356  0.          1.        ]
 [ 2.23606798  1.          0.        ]]

この関数については この文書 で詳細を調べてください。

似たような関数として (scipy.spatial.distance.cdist) がある。 これは2つの点集合のすべての組み合わせに対して距離を計算する。これについては この文書 を参照のこと.

Matplotlib

Matplotlibは、描画のためのライブラリです。 ここでは、 matplotlib.pyplotモジュールについて簡単に紹介します。 これは、MATLABの描画機能と類似のシステムを提供しています。

プロット

matplotlibで最も重要な関数が plotです。これにより 2Dデータをプロットすることができます。例を示します:

In [46]:
%matplotlib inline
# 上により、matplotlibで作成した図がブラウザ上に表示されるようになる

import numpy as np
import matplotlib.pyplot as plt

# x座標に対するsin関数の値を計算しy座標の値とする
x = np.arange(0, 3 * np.pi, 0.05)
y = np.sin(x)

# x軸を引く
plt.axhline(y=0, xmin=min(x), xmax=max(x), color='green',linestyle='dashed')

# matplotlibを使ってプロットする
plt.plot(x, y)
plt.show()  # 本来はこの plt.show() を呼ばなければ図が見えない
In [44]:
print(x)
[ 0.    0.05  0.1   0.15  0.2   0.25  0.3   0.35  0.4   0.45  0.5   0.55
  0.6   0.65  0.7   0.75  0.8   0.85  0.9   0.95  1.    1.05  1.1   1.15
  1.2   1.25  1.3   1.35  1.4   1.45  1.5   1.55  1.6   1.65  1.7   1.75
  1.8   1.85  1.9   1.95  2.    2.05  2.1   2.15  2.2   2.25  2.3   2.35
  2.4   2.45  2.5   2.55  2.6   2.65  2.7   2.75  2.8   2.85  2.9   2.95
  3.    3.05  3.1   3.15  3.2   3.25  3.3   3.35  3.4   3.45  3.5   3.55
  3.6   3.65  3.7   3.75  3.8   3.85  3.9   3.95  4.    4.05  4.1   4.15
  4.2   4.25  4.3   4.35  4.4   4.45  4.5   4.55  4.6   4.65  4.7   4.75
  4.8   4.85  4.9   4.95  5.    5.05  5.1   5.15  5.2   5.25  5.3   5.35
  5.4   5.45  5.5   5.55  5.6   5.65  5.7   5.75  5.8   5.85  5.9   5.95
  6.    6.05  6.1   6.15  6.2   6.25  6.3   6.35  6.4   6.45  6.5   6.55
  6.6   6.65  6.7   6.75  6.8   6.85  6.9   6.95  7.    7.05  7.1   7.15
  7.2   7.25  7.3   7.35  7.4   7.45  7.5   7.55  7.6   7.65  7.7   7.75
  7.8   7.85  7.9   7.95  8.    8.05  8.1   8.15  8.2   8.25  8.3   8.35
  8.4   8.45  8.5   8.55  8.6   8.65  8.7   8.75  8.8   8.85  8.9   8.95
  9.    9.05  9.1   9.15  9.2   9.25  9.3   9.35  9.4 ]
In [51]:
# sin関数を散布図として表示する
x = np.arange(0, 3 * np.pi, 0.05)
y = np.sin(x)
plt.plot(x,y,'b.') # 第3引数に色(Blue)とマーカーの種類を指定すると散布図となる

plt.grid()     # 枠線を書く
plt.show()

ちょっとだけ余計に作業することで、簡単にいくつもの線を同時に描画したり、題目をつけたり、 凡例をつけたり、軸のラベルをつけたりできます:

In [20]:
import numpy as np
import matplotlib.pyplot as plt

# x座標に対するsin関数の値を計算しy座標の値とする
x = np.arange(0, 3 * np.pi, 0.1)
y_sin = np.sin(x)
y_cos = np.cos(x)

# matplotlibを使ってプロットする
plt.plot(x, y_sin)
plt.plot(x, y_cos)
plt.xlabel('x axis label')  # ラベルに日本語は使えないと思って良い
plt.ylabel('y axis label')
plt.title('Sine and Cosine')
plt.legend(['Sine', 'Cosine'])
plt.grid()
plt.show()     #  これを最後に呼ぶこと

plot 関数について、詳しくは この文書をご覧ください

Subplot

いろいろな図を同じ画面に表示するにはsubplot関数を使います。 次の例をみてください:

In [22]:
import numpy as np
import matplotlib.pyplot as plt

# x座標に対応するsinとcos関数の値を計算する 
x = np.arange(0, 3 * np.pi, 0.1)
y_sin = np.sin(x)
y_cos = np.cos(x)

# subplotの枠を設定:高さ2、幅1とする,
# そして、最初の枠をactiveにする
plt.subplot(2, 1, 1)

# 最初の図を描画する
plt.plot(x, y_sin)
plt.title('Sine')
plt.grid()

# 2番めの枠をactiveにして、別な図を描画する
plt.subplot(2, 1, 2)
plt.plot(x, y_cos)
plt.title('Cosine')
plt.grid()

# 図を表示する(これを最後に呼ぶこと)
plt.show()

subplot 関数について、詳しくは この文書をご覧ください

画像

imshow 関数を使って画像を表示できます。次の例をみてください:

In [28]:
import numpy as np
from scipy.misc import imread, imresize
import matplotlib.pyplot as plt

img = imread('assets/cat.jpg')
img_tinted = img * [0.5, 0.5, 0.5]

# 原画像を表示
plt.subplot(1, 2, 1)
plt.imshow(img)

# 色を変えた画像を表示
plt.subplot(1, 2, 2)

# A slight gotcha with imshow is that it might give strange results
# if presented with data that is not uint8. To work around this, we
# explicitly cast the image to uint8 before displaying it.
plt.imshow(np.uint8(img_tinted))
plt.show()

練習問題


1. 乱数を用いて、3x3の行列を作れ。要素はすべて乱数によって生成された値とする。この行列を変数Aの値とせよ。 つぎに、np.linalg.matrix_rank 関数にて行列Aの階数を求めよ。ここでもしも階数が3でなければAが正則行列でない(つまり逆行列が存在しない)ので、最初からやり直しすること。

In [65]:
# 問題1
import numpy as np
A = np.random.random((3,3))
print(A)
print( np.linalg.matrix_rank(A)  )
[[ 0.75078998  0.66571317  0.2969129 ]
 [ 0.96454726  0.67272454  0.81742854]
 [ 0.7499657   0.21511015  0.49049703]]
3


2. 前問Aの逆行列を求め、変数Ainvの値とせよ。そして、AとAinvの(行列)積が単位行列になることを確かめよ。

In [66]:
# 問題2
aMatrix = np.matrix(A)
Ainv = aMatrix.I
print(aMatrix.dot(Ainv))
[[  1.00000000e+00  -1.43551457e-16   1.02024117e-16]
 [  1.48718626e-16   1.00000000e+00   1.94088126e-17]
 [  2.02352779e-16   6.39657000e-17   1.00000000e+00]]


3. 行列の固有値と固有ベクトルは次のようにして求めることができる(参考:固有値、固有ベクトルを可視化する)

# Mを行列とする
eigenVals, eigenVecs = np.linalg.eig(M)   # eigenValsが固有値、eigenVecsが固有ベクトルの行列

これに基づき、Aの固有値と固有ベクトルを求めよ。(注意: 固有ベクトルは、固有ベクトルの行列の「列」ベクトルである)

In [67]:
# 問題3
eigenVals, eigenVecs = np.linalg.eig(A) 
print(eigenVals)
print(eigenVecs )
[ 1.86144624+0.j          0.02628266+0.25326414j  0.02628266-0.25326414j]
[[ 0.54688314+0.j          0.39220388-0.37259025j  0.39220388+0.37259025j]
 [ 0.72802476+0.j         -0.00314007+0.55470668j -0.00314007-0.55470668j]
 [ 0.41339906+0.j         -0.63217344+0.j         -0.63217344-0.j        ]]


4. 3x3の行列に対し、固有値とそれに対応する固有ベクトルは三組ある。元の行列を$\mathbf{A}$、 固有値と対応する固有ベクトルをそれぞれ $\lambda$と$\mathbf{x}$ と書くと、 $\mathbf{Ax}=\lambda\mathbf{x}$ が成り立つことを確認せよ。また、固有ベクトルがみな正則化されている(長さが1)であることも確認せよ。

In [68]:
# 問題4
for n in range(3):
   print(eigenVecs[:,n])
   print(aMatrix.dot(eigenVecs[:,n]))
   print(eigenVals[n] * eigenVecs[:,n])
   print(np.sum(eigenVecs[:,n]**2))
[ 0.54688314+0.j  0.72802476+0.j  0.41339906+0.j]
[[ 1.01799356+0.j  1.35517895+0.j  0.76952012+0.j]]
[ 1.01799356+0.j  1.35517895+0.j  0.76952012+0.j]
(1+0j)
[ 0.39220388-0.37259025j -0.00314007+0.55470668j -0.63217344+0.j        ]
[[ 0.10467191+0.08953852j -0.14056984+0.0137839j  -0.01661520-0.16010687j]]
[ 0.10467191+0.08953852j -0.14056984+0.0137839j  -0.01661520-0.16010687j]
(0.106954016781-0.295746322229j)
[ 0.39220388+0.37259025j -0.00314007-0.55470668j -0.63217344-0.j        ]
[[ 0.10467191-0.08953852j -0.14056984-0.0137839j  -0.01661520+0.16010687j]]
[ 0.10467191-0.08953852j -0.14056984-0.0137839j  -0.01661520+0.16010687j]
(0.106954016781+0.295746322229j)


5. 次を実行し、散布図をmatplotlibを用いて以下のような図を表示せよ。

x=np.array([2.8, 2.9, 3.0, 3.1, 3.2, 3.2, 3.2, 3.3, 3.4])
y=np.array([30,26,33,31,33,35,37,36,33])

In [3]:
# 問題5
%matplotlib inline
import matplotlib.pyplot as plt
x=np.array([2.8, 2.9, 3.0, 3.1, 3.2, 3.2, 3.2, 3.3, 3.4])
y=np.array([30,26,33,31,33,35,37,36,33])
plt.scatter(x,y)
Out[3]:
<matplotlib.collections.PathCollection at 0x4fb1f10>


6. 最小二乗法を用いて、上記のxとyの関係を近似する線形関数$y = ax + b$を求めよう(つまり、傾き$a$と切片$b$を求める)。
これは、上記のデータを$x_i, y_i$ (i=1,..,9)とおくと、 $ \left( \begin{array}{r} a \\ b \end{array} \right) = \left( \begin{array}{rr} \sum_i x_i^2 & \sum_i x_i \\ \sum_i x_i & \sum_i 1 \\ \end{array} \right)^{-1} \left( \begin{array}{r} \sum_i x_i y_i \\ \sum_i y_i \end{array} \right) $ という式から得られる。これによって近似線形関数を求め、上の図にその直線を書き込め。次のような図が得られれば良い。 注意: 上記を満たす関数は後で学ぶ。ここでは上記の方法にしたがって解いてみよ。

In [4]:
# 問題6
xSq = np.sum(x**2)
xSum = np.sum(x)
n = len(x)
xySum = np.sum(x*y)
ySum = np.sum(y)
A = (np.matrix([[xSq, xSum],[xSum,n]]).I)
B = np.matrix([xySum,ySum]).T
print(A)
print(B)
C =np.array(A.dot(B))
a = C[0][0] ; b = C[1][0]
print(a,b)

plt.scatter(x,y)
plt.plot([min(x),max(x)],[a*min(x)+b, a*max(x)+b],color="g")
[[  3.38345865 -10.56390977]
 [-10.56390977  33.09398496]]
[[ 921.5]
 [ 294. ]]
12.0676691729 -5.01127819549
Out[4]:
[<matplotlib.lines.Line2D at 0x5113fd0>]
In [44]:
# 問題6の追加
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def lsm(x,y):
    xxsum=np.sum(pow(x,2))
    xsum=np.sum(x)
    xysum=np.sum(x*y)
    ysum=np.sum(y)
    n=len(x)
    x=np.matrix([[xxsum,xsum],[xsum,n]])
    x=x.I
    y=np.matrix([[xysum,ysum]])
    y=y.T
    ab=np.array(x.dot(y))
    return ab
    

x=np.random.randint(1, 100, 10)
y=np.random.randint(1, 100, 10)
plt.scatter(x,y)
print(x,y)

ab=lsm(x,y)
a= ab[0]
b= ab[1]
print(a,b)

plt.scatter(x,y)
plt.plot([min(x),max(x)],[a*min(x)+b, a*max(x)+b],color="g")
[30 96 74 43 32 19 66 94 12 92] [33 63 39 25 32 58 81 66 51 83]
[ 0.34703254] [ 33.73558422]
Out[44]:
[<matplotlib.lines.Line2D at 0x6a50c70>]