NumPy & SciPy メモ

2018年12月24日

はじめに

Python モジュール Numpy & SciPy のメモ。

準備

NumPy のインポート

In [1]:
import numpy as np

NumPy の配列

NumPy には多次元配列を扱うための ndarray というクラスがある。ndarray のオブジェクトは array() 関数により生成する。1次元の配列を生成する場合には次のようにする。

In [2]:
np.array([1., 2., 3.])
Out[2]:
array([1., 2., 3.])

2 次元の配列を生成する場合には次のようにする。

In [3]:
np.array([[1., 2.], [3., 4.]])
Out[3]:
array([[1., 2.],
       [3., 4.]])

2 次元の配列の特殊な場合として、matrix クラスが用意されている。

In [4]:
np.matrix([[1., 2.], [3., 4.]])
Out[4]:
matrix([[1., 2.],
        [3., 4.]])

結合

リストや配列どうしを結合して配列を作ることができる。

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

縦に結合することもできる。

In [7]:
a = np.array([[1.], [2.]])
b = np.array([[3.], [4.]])
np.vstack((a, b))
Out[7]:
array([[1.],
       [2.],
       [3.],
       [4.]])

ベクトル

定義

ベクトルは 1 次元の配列で表す。

In [8]:
a = np.array([1., 2., 3.])
a
Out[8]:
array([1., 2., 3.])

行列で表してもよい。

In [9]:
a = np.matrix([1., 2., 3.]).T
a
Out[9]:
matrix([[1.],
        [2.],
        [3.]])

行列の要素は文字列で指定することもできる。

In [10]:
a = np.matrix('1., 2., 3.').T
a
Out[10]:
matrix([[1.],
        [2.],
        [3.]])

np.linspace() を使えば、次のように 0 から 1 まで 11 個の値をもつ配列といったものを作ることができる。

In [11]:
a = np.linspace(0, 1, 11)
a
Out[11]:
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])

ある範囲の整数の要素を持つ配列を作る場合は、np.arange() が使える。

In [12]:
a = np.arange(5, 11)
a
Out[12]:
array([ 5,  6,  7,  8,  9, 10])

成分

In [13]:
a = np.array([1., 2., 3.])
a[1]
Out[13]:
2.0

添え字の範囲を指定して配列の部分を取り出すことができる。

In [14]:
a[1:3]
Out[14]:
array([2., 3.])

省略もできる。

In [15]:
a[:2]
Out[15]:
array([1., 2.])
In [16]:
a[:]
Out[16]:
array([1., 2., 3.])

代入

In [17]:
a = np.array([1., 2., 3.])
b = np.array(a)
b
Out[17]:
array([1., 2., 3.])
In [18]:
b = a.copy()
b
Out[18]:
array([1., 2., 3.])

スカラーとの和・差

In [19]:
a = np.array([1., 2., 3.])
a + 1
Out[19]:
array([2., 3., 4.])
In [20]:
a - 1
Out[20]:
array([0., 1., 2.])

スカラーとの積・除算

In [21]:
a = np.array([1., 2., 3.])
a*2
Out[21]:
array([2., 4., 6.])
In [22]:
a/2.
Out[22]:
array([0.5, 1. , 1.5])

ベクトルどうしの和・差

In [23]:
a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])
a + b
Out[23]:
array([5., 7., 9.])
In [24]:
b - a
Out[24]:
array([3., 3., 3.])

ベクトルの成分どうしの積・除算

In [25]:
a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])
a*b
Out[25]:
array([ 4., 10., 18.])
In [26]:
b/a
Out[26]:
array([4. , 2.5, 2. ])

ノルム

In [27]:
a = np.array([1., 2., 3.])
np.linalg.norm(a)
Out[27]:
3.7416573867739413

内積

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

次のようにしてもよい。

In [29]:
a.dot(b)
Out[29]:
32.0

内積を使うと、ノルムを次のように計算できる。

In [30]:
np.sqrt(np.dot(a, a))
Out[30]:
3.7416573867739413

外積

In [31]:
a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])
np.cross(a, b)
Out[31]:
array([-3.,  6., -3.])

総和・平均

In [32]:
a = np.array([1., 2., 3.])
np.sum(a)
Out[32]:
6.0
In [33]:
a.sum()
Out[33]:
6.0
In [34]:
np.average(a)
Out[34]:
2.0

最大・最小

In [35]:
np.max(a)
Out[35]:
3.0
In [36]:
np.min(a)
Out[36]:
1.0
In [37]:
a.max()
Out[37]:
3.0
In [38]:
a.min()
Out[38]:
1.0

論理

配列に比較演算子を適用することができる。

In [39]:
a = np.linspace(1, 10, 10)
a % 2 == 0
Out[39]:
array([False,  True, False,  True, False,  True, False,  True, False,
        True])

これを用いて、条件を満たす要素だけを取り出すことができる。

In [40]:
a[a % 2 == 0]
Out[40]:
array([ 2.,  4.,  6.,  8., 10.])

and などは直接使えないので、条件を段階的に適用する必要がある。たとえば、5 以上の偶数を取り出すには次のようにする。

In [41]:
a2 = a[a > 5]
a2[a2 % 2 == 0]
Out[41]:
array([ 6.,  8., 10.])

np.all() ですべて真かどうか、np.any() で真があるかどうかを調べることができる。

In [42]:
np.all(a % 2 == 0)
Out[42]:
False
In [43]:
np.any(a % 2 == 0)
Out[43]:
True

行列

定義

In [44]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
A
Out[44]:
matrix([[1., 2., 3.],
        [4., 5., 6.],
        [7., 8., 9.]])

行列の成分を文字列で指定することもできる。

In [45]:
A = np.matrix('''
    1., 2., 3.;
    4., 5., 6.;
    7., 8., 9.
''')
A
Out[45]:
matrix([[1., 2., 3.],
        [4., 5., 6.],
        [7., 8., 9.]])

matrix の代りに短く mat でもよい。

In [46]:
A = np.mat([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
A
Out[46]:
matrix([[1., 2., 3.],
        [4., 5., 6.],
        [7., 8., 9.]])

単位行列

In [47]:
np.eye(3)
Out[47]:
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
In [48]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
np.eye(3)*A
Out[48]:
matrix([[1., 2., 3.],
        [4., 5., 6.],
        [7., 8., 9.]])

ゼロ行列

In [49]:
np.zeros([3, 3])
Out[49]:
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
In [50]:
np.zeros([2, 3])
Out[50]:
array([[0., 0., 0.],
       [0., 0., 0.]])

行列のサイズは shape() で得られる。

In [51]:
A.shape
Out[51]:
(3, 3)

成分

In [52]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
A[:, 1]
Out[52]:
matrix([[2.],
        [5.],
        [8.]])

添え字の範囲を指定して行列の部分を取り出すことができる。

In [53]:
A[0:2, 1:3]
Out[53]:
matrix([[2., 3.],
        [5., 6.]])

添え字のリストを指定することもできる。

In [54]:
A[:, [0, 2]]
Out[54]:
matrix([[1., 3.],
        [4., 6.],
        [7., 9.]])

代入

In [55]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
B = np.matrix(A)
B
Out[55]:
matrix([[1., 2., 3.],
        [4., 5., 6.],
        [7., 8., 9.]])
In [56]:
B = A.copy()
B
Out[56]:
matrix([[1., 2., 3.],
        [4., 5., 6.],
        [7., 8., 9.]])

スカラーとの和・差

In [57]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
A + 1
Out[57]:
matrix([[ 2.,  3.,  4.],
        [ 5.,  6.,  7.],
        [ 8.,  9., 10.]])
In [58]:
A - 1
Out[58]:
matrix([[0., 1., 2.],
        [3., 4., 5.],
        [6., 7., 8.]])

スカラーとの積・除算

In [59]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
A*2.
Out[59]:
matrix([[ 2.,  4.,  6.],
        [ 8., 10., 12.],
        [14., 16., 18.]])
In [60]:
A/2.
Out[60]:
matrix([[0.5, 1. , 1.5],
        [2. , 2.5, 3. ],
        [3.5, 4. , 4.5]])

行列どうしの和・差

In [61]:
A = np.array([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
B = np.array([
    [10., 11., 12.],
    [13., 14., 15.],
    [16., 17., 18.]
])
A + B
Out[61]:
array([[11., 13., 15.],
       [17., 19., 21.],
       [23., 25., 27.]])
In [62]:
B - A
Out[62]:
array([[9., 9., 9.],
       [9., 9., 9.],
       [9., 9., 9.]])

行列の成分どうしの積・除算

array どうしの積は成分どうしの積になる。

In [63]:
A = np.array([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
B = np.array([
    [10., 11., 12.],
    [13., 14., 15.],
    [16., 17., 18.]
])
A*B
Out[63]:
array([[ 10.,  22.,  36.],
       [ 52.,  70.,  90.],
       [112., 136., 162.]])

matrix どうしの積は行列の積になる。

In [64]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
B = np.matrix([
    [10., 11., 12.],
    [13., 14., 15.],
    [16., 17., 18.]
])
A*B
Out[64]:
matrix([[ 84.,  90.,  96.],
        [201., 216., 231.],
        [318., 342., 366.]])

matrix の成分どうしで積を取りたい場合は、array に変換してやればよい。

In [65]:
np.matrix(np.array(A)*np.array(B))
Out[65]:
matrix([[ 10.,  22.,  36.],
        [ 52.,  70.,  90.],
        [112., 136., 162.]])

除算は普通にできる。

In [66]:
B/A
Out[66]:
matrix([[10.        ,  5.5       ,  4.        ],
        [ 3.25      ,  2.8       ,  2.5       ],
        [ 2.28571429,  2.125     ,  2.        ]])

行列の積

In [67]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
B = np.matrix([
    [10., 11., 12.],
    [13., 14., 15.],
    [16., 17., 18.]
])
A*B
Out[67]:
matrix([[ 84.,  90.,  96.],
        [201., 216., 231.],
        [318., 342., 366.]])

dot() 関数を使ってもよい。

In [68]:
np.dot(A, B)
Out[68]:
matrix([[ 84.,  90.,  96.],
        [201., 216., 231.],
        [318., 342., 366.]])

次のような書き方もできる。

In [69]:
A.dot(B)
Out[69]:
matrix([[ 84.,  90.,  96.],
        [201., 216., 231.],
        [318., 342., 366.]])

array の場合、"*" の代わりに "@" を使えば行列の積になる。

In [70]:
a = np.array(A)
b = np.array(B)
a*b
Out[70]:
array([[ 10.,  22.,  36.],
       [ 52.,  70.,  90.],
       [112., 136., 162.]])
In [71]:
a@b
Out[71]:
array([[ 84.,  90.,  96.],
       [201., 216., 231.],
       [318., 342., 366.]])

dot() 関数はそのまま使える。

In [72]:
np.dot(a, b)
Out[72]:
array([[ 84.,  90.,  96.],
       [201., 216., 231.],
       [318., 342., 366.]])

行列の累乗

In [73]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
A**2
Out[73]:
matrix([[ 30.,  36.,  42.],
        [ 66.,  81.,  96.],
        [102., 126., 150.]])

転置

In [74]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
A.T
Out[74]:
matrix([[1., 4., 7.],
        [2., 5., 8.],
        [3., 6., 9.]])

トレース

In [75]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
np.trace(A)
Out[75]:
15.0

対角成分

In [76]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 9.]
])
np.diag(A)
Out[76]:
array([1., 5., 9.])

対角の和と取るとトレースになる。

In [77]:
np.sum(np.diag(A))
Out[77]:
15.0

行列式

In [78]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 10.]
])
np.linalg.det(A)
Out[78]:
-3.000000000000001

逆行列

In [79]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 10.]
])
np.linalg.inv(A)
Out[79]:
matrix([[-0.66666667, -1.33333333,  1.        ],
        [-0.66666667,  3.66666667, -2.        ],
        [ 1.        , -2.        ,  1.        ]])
In [80]:
A*np.linalg.inv(A)
Out[80]:
matrix([[ 1.00000000e+00, -4.44089210e-16, -1.11022302e-16],
        [ 4.44089210e-16,  1.00000000e+00, -2.22044605e-16],
        [ 4.44089210e-16,  8.88178420e-16,  1.00000000e+00]])

余因子行列を用いた逆行列の計算

In [81]:
def cofactor(A):
    def M(A, i, j):
        idx = np.arange(0, A.shape[0])
        return A[idx != i, :][:, idx != j]

    A_tilde = np.matrix(np.zeros(A.shape))
    for i in range(A.shape[0]):
        for j in range(A.shape[1]):
            A_tilde[j, i] = (-1)**(i+j)*np.linalg.det(M(A, i, j))
    return A_tilde

cofactor(A)/np.linalg.det(A)
Out[81]:
matrix([[-0.66666667, -1.33333333,  1.        ],
        [-0.66666667,  3.66666667, -2.        ],
        [ 1.        , -2.        ,  1.        ]])

連立方程式の解

In [82]:
b = np.array([1., 2., 3.])
x = np.linalg.solve(A, b)
x
Out[82]:
array([-3.33333333e-01,  6.66666667e-01,  3.17206578e-17])
In [83]:
A.dot(x)
Out[83]:
matrix([[1., 2., 3.]])

固有値・固有ベクトル

In [84]:
A = np.matrix([
    [1., 2., 3.],
    [4., 5., 6.],
    [7., 8., 10.]
])
lam, P = np.linalg.eig(A)
P, lam
Out[84]:
(matrix([[-0.22351336, -0.86584578,  0.27829649],
         [-0.50394563,  0.0856512 , -0.8318468 ],
         [-0.83431444,  0.4929249 ,  0.48018951]]),
 array([16.70749332, -0.90574018,  0.19824686]))

対角化

In [85]:
np.linalg.inv(P)*A*P
Out[85]:
matrix([[ 1.67074933e+01, -7.20607628e-15, -1.59878008e-15],
        [-3.28519015e-15, -9.05740180e-01, -1.83986877e-16],
        [ 3.54501614e-15, -1.06062765e-15,  1.98246863e-01]])

最適化

非線形方程式を解く

scipy.optimize の newton() で Newton 法により非線形方程式を解くことができる。

In [86]:
from scipy.optimize import newton

newton(lambda x: x**2 - 2., 1.)
Out[86]:
1.4142135623730947

最小化

scipy.optimize の fmin() 関数を最小化する変数を求めることができる。

In [87]:
from scipy.optimize import fmin

fmin(lambda x: abs(x**2 - 2.), 1.)

Optimization terminated successfully.
         Current function value: 0.000125
         Iterations: 14
         Function evaluations: 28
Out[87]:
array([1.41425781])

カーブフィッティング

scipy.optimize の curve_fit() で任意の関数の係数を決定することができる。

In [88]:
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

from scipy.optimize import curve_fit

point = [
    [10, 1.],
    [50, 10.],
    [100, 30.],
    [200, 100.],
]

data = np.array(point)
data_x = data[:, 0]
data_y = data[:, 1]

def f(x, a, b):
    return  a*x**b

opt, cov = curve_fit(f, data_x, data_y)
a, b = opt
print("a = {}, b = {}".format(a, b))

x = np.linspace(0, 200, 100)
plt.plot(data_x, data_y, "o")
plt.plot(x, f(x, a, b), label="{0:.3e} x^{1:.3f}".format(a, b))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

a = 0.011393733856615188, b = 1.7135856775825407

多項式補間

np.polyfit() で多項式補間の係数を求めることができる。np.poly1d() で求めた係数から関数を作ることができる。

In [89]:
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

point = [
    [10, 1.],
    [50, 10.],
    [100, 30.],
    [200, 100.],
]

data = np.array(point)
data_x = data[:, 0]
data_y = data[:, 1]

n = 2
a = np.polyfit(data_x, data_y, n)
print("a = {}".format(a))
f = np.poly1d(a)

x = 100.
y = 0.
for i in range(n+1):
    y = y*x + a[i]

print("x = 100 : ", f(x), y)

x = np.linspace(0, 200, 100)
plt.plot(data_x, data_y, "o")
plt.plot(x, f(x), label="interpolate")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

a = [ 0.00198543  0.10395966 -0.21479268]
x = 100 :  30.03548748599179 30.03548748599179

AIC (赤池情報基準) を用いた補間式の決定

In [90]:
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

import sys

point = [
    [10, 1.],
    [50, 10.],
    [100, 30.],
    [200, 100.],
    [250, 190.],
    [400, 240.],
    [500, 420.],
]

data = np.array(point)
data_x = data[:, 0]
data_y = data[:, 1]

def calc_aic(x, y, yp, k):
    s2 = np.average((y - yp)**2)
    p = np.exp(-(y - yp)**2 / (2.*s2)) / np.sqrt(2.*np.pi*s2)
    return -2.*sum(np.log(p)) + 2.*k

min_aic = sys.float_info.max
min_n = 0
min_a = []
min_yp = []

for n in range(1, 4):
    a = np.polyfit(data_x, data_y, n)
    f = np.poly1d(a)
    yp = f(data_x)

    aic = calc_aic(data_x, data_y, yp, n + 2)

    print("{}: AIC = {:5.0f}".format(n, aic))

    if aic < min_aic:
        min_aic = aic
        min_n = n
        min_a = a
        min_yp = yp

print("-----")
print("n = {}".format(min_n))
print("a =")
print(min_a)

plt.plot(data_x, data_y, "o")

xp = np.linspace(0, 500, 100)
yp = np.poly1d(min_a)(xp)
plt.plot(xp, yp, label="interpolate {} degree".format(min_n))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

1: AIC =    74
2: AIC =    73
3: AIC =    75
-----
n = 2
a =
[ 8.12569022e-04  4.09239541e-01 -7.66167142e+00]

AIC が小さいほどデータに対して当てはまりがよい傾向にあると考えるが、小さければよいとも限らない。グラフの形や未知のデータとの当てはまり具合などを見て判断する。上記の例だとデータが少ないので、あまり次数を多くとってはいけない。データ点を通りさえすればよいなら、次数をデータ点数-1 とすればよいが、ふつうはぐにゃぐにゃのグラフになって使い物にならない。

データのスムージング

ノイズを含むデータをローパスフィルターでスムージングする。scipy の fftpack を用いる。

In [91]:
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

from scipy import fftpack

np.random.seed(0)
x = np.linspace(0, 2.*np.pi, 100)
y = np.sin(x) + np.random.normal(0, 0.1, x.size)

fc = 0.01 # cut off

datas = np.copy(data)

fs = fftpack.fftfreq(x.size, d=1)
F = fftpack.fft(y)
F[np.abs(fs) > fc] = 0.
f = np.real(fftpack.ifft(F))

plt.plot(x, y, label="data")
plt.plot(x, f, label="smoothed")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()