NumPy & SciPy メモ
2022年11月13日
はじめに
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]:
-2.9999999999999996
逆行列
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, 0.00000000e+00, 1.11022302e-16],
[0.00000000e+00, 1.00000000e+00, 2.22044605e-16],
[8.88178420e-16, 0.00000000e+00, 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([-0.33333333, 0.66666667, -0. ])
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, -1.64862768e-15, 1.42386260e-16],
[-6.37415816e-16, -9.05740180e-01, -5.09618857e-16],
[ 1.78090322e-15, -8.31326629e-16, 1.98246863e-01]])
最適化
非線形方程式を解く
scipy.optimize の newton() で Newton 法により非線形方程式を解くことができる。
In [86]:
from scipy.optimize import newton
newton(lambda x: x**2 - 2., 1.)
Out[86]:
1.414213562373095
最小化
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
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.grid()
plt.legend()
plt.show()
a = 0.011393733856615188, b = 1.7135856775825407
線形補間
2つのデータ点から線形補間するには、scipy.interpolate.interp1d() を用いる。
In [89]:
%matplotlib inline
from matplotlib import pyplot as plt
from scipy import interpolate
point = [
[50, 10.],
[100, 30.],
]
data = np.array(point)
data_x = data[:, 0]
data_y = data[:, 1]
f = interpolate.interp1d(data_x, data_y, fill_value="extrapolate")
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.grid()
plt.legend()
plt.show()
データの間を線形補間したい場合は、np.interp() を用いる。
In [90]:
%matplotlib inline
from matplotlib import pyplot as plt
point = [
[10, 1.],
[50, 10.],
[100, 30.],
[200, 100.],
]
data = np.array(point)
data_x = data[:, 0]
data_y = data[:, 1]
f = lambda x: np.interp(x, data_x, data_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.grid()
plt.legend()
plt.show()
多項式補間
np.polyfit() で多項式補間の係数を求めることができる。np.poly1d() で求めた係数から関数を作ることができる。
In [91]:
%matplotlib inline
from matplotlib import pyplot as plt
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.grid()
plt.legend()
plt.show()
a = [ 0.00198543 0.10395966 -0.21479268] x = 100 : 30.035487485991776 30.035487485991776
AIC (赤池情報基準) を用いた補間式の決定
In [92]:
%matplotlib inline
from matplotlib import pyplot as plt
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.grid()
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 とすればよいが、ふつうはぐにゃぐにゃのグラフになって使い物にならない。
積分
In [93]:
from scipy import integrate
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)
def F(x):
y = 0.
for i in range(n+1):
y += 1/(n+1 - i)*a[i]*x**(n+1 - i)
return y
x = np.linspace(0, 200, 100)
print(integrate.trapz(f(x), x=x))
print(F(200.) - F(0.))
a = [ 0.00198543 0.10395966 -0.21479268] 7330.988565111638 7330.718465944461
データのスムージング
ノイズを含むデータをローパスフィルターでスムージングする。scipy の fftpack を用いる。
In [94]:
%matplotlib inline
from matplotlib import pyplot as plt
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
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.grid()
plt.legend()
plt.show()
移動平均
In [95]:
%matplotlib inline
from matplotlib import pyplot as plt
np.random.seed(0)
x = np.linspace(0, 2.*np.pi, 100)
y = np.sin(x) + np.random.normal(0, 0.1, x.size)
def moveing_average(y, n):
return np.convolve(y, np.ones(n), "same")/n
f = moveing_average(y, 10)
plt.plot(x, y, label="data")
plt.plot(x, f, label="smoothed")
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.legend()
plt.show()
包絡線
In [96]:
%matplotlib inline
from matplotlib import pyplot as plt
from scipy import signal
x = np.linspace(0, 2.*np.pi, 1000)
y = np.sin(2.*np.pi*3.*x)*np.sin(2.*x)
f = np.abs(signal.hilbert(y))
plt.plot(x, y, label="data")
plt.plot(x, f, label="envelope")
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.legend()
plt.show()
