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]:
2 次元の配列を生成する場合には次のようにする。 In [3]:
np.array([[1., 2.], [3., 4.]])
Out[3]:
2 次元の配列の特殊な場合として、matrix クラスが用意されている。 In [4]:
np.matrix([[1., 2.], [3., 4.]])
Out[4]:
結合リストや配列どうしを結合して配列を作ることができる。 In [5]:
np.hstack(([1, 2],[3, 4]))
Out[5]:
In [6]:
a = np.array([1., 2.])
b = np.array([3., 4.])
np.hstack((a, b))
Out[6]:
縦に結合することもできる。 In [7]:
a = np.array([[1.], [2.]])
b = np.array([[3.], [4.]])
np.vstack((a, b))
Out[7]:
ベクトル定義ベクトルは 1 次元の配列で表す。 In [8]:
a = np.array([1., 2., 3.])
a
Out[8]:
行列で表してもよい。 In [9]:
a = np.matrix([1., 2., 3.]).T
a
Out[9]:
行列の要素は文字列で指定することもできる。 In [10]:
a = np.matrix('1., 2., 3.').T
a
Out[10]:
np.linspace() を使えば、次のように 0 から 1 まで 11 個の値をもつ配列といったものを作ることができる。 In [11]:
a = np.linspace(0, 1, 11)
a
Out[11]:
ある範囲の整数の要素を持つ配列を作る場合は、np.arange() が使える。 In [12]:
a = np.arange(5, 11)
a
Out[12]:
成分In [13]:
a = np.array([1., 2., 3.])
a[1]
Out[13]:
添え字の範囲を指定して配列の部分を取り出すことができる。 In [14]:
a[1:3]
Out[14]:
省略もできる。 In [15]:
a[:2]
Out[15]:
In [16]:
a[:]
Out[16]:
代入In [17]:
a = np.array([1., 2., 3.])
b = np.array(a)
b
Out[17]:
In [18]:
b = a.copy()
b
Out[18]:
スカラーとの和・差In [19]:
a = np.array([1., 2., 3.])
a + 1
Out[19]:
In [20]:
a - 1
Out[20]:
スカラーとの積・除算In [21]:
a = np.array([1., 2., 3.])
a*2
Out[21]:
In [22]:
a/2.
Out[22]:
ベクトルどうしの和・差In [23]:
a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])
a + b
Out[23]:
In [24]:
b - a
Out[24]:
ベクトルの成分どうしの積・除算In [25]:
a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])
a*b
Out[25]:
In [26]:
b/a
Out[26]:
ノルムIn [27]:
a = np.array([1., 2., 3.])
np.linalg.norm(a)
Out[27]:
内積In [28]:
a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])
np.dot(a, b)
Out[28]:
次のようにしてもよい。 In [29]:
a.dot(b)
Out[29]:
内積を使うと、ノルムを次のように計算できる。 In [30]:
np.sqrt(np.dot(a, a))
Out[30]:
外積In [31]:
a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])
np.cross(a, b)
Out[31]:
総和・平均In [32]:
a = np.array([1., 2., 3.])
np.sum(a)
Out[32]:
In [33]:
a.sum()
Out[33]:
In [34]:
np.average(a)
Out[34]:
最大・最小In [35]:
np.max(a)
Out[35]:
In [36]:
np.min(a)
Out[36]:
In [37]:
a.max()
Out[37]:
In [38]:
a.min()
Out[38]:
論理配列に比較演算子を適用することができる。 In [39]:
a = np.linspace(1, 10, 10)
a % 2 == 0
Out[39]:
これを用いて、条件を満たす要素だけを取り出すことができる。 In [40]:
a[a % 2 == 0]
Out[40]:
and などは直接使えないので、条件を段階的に適用する必要がある。たとえば、5 以上の偶数を取り出すには次のようにする。 In [41]:
a2 = a[a > 5]
a2[a2 % 2 == 0]
Out[41]:
np.all() ですべて真かどうか、np.any() で真があるかどうかを調べることができる。 In [42]:
np.all(a % 2 == 0)
Out[42]:
In [43]:
np.any(a % 2 == 0)
Out[43]:
行列定義In [44]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
A
Out[44]:
行列の成分を文字列で指定することもできる。 In [45]:
A = np.matrix('''
1., 2., 3.;
4., 5., 6.;
7., 8., 9.
''')
A
Out[45]:
matrix の代りに短く mat でもよい。 In [46]:
A = np.mat([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
A
Out[46]:
単位行列In [47]:
np.eye(3)
Out[47]:
In [48]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
np.eye(3)*A
Out[48]:
ゼロ行列In [49]:
np.zeros([3, 3])
Out[49]:
In [50]:
np.zeros([2, 3])
Out[50]:
行列のサイズは shape() で得られる。 In [51]:
A.shape
Out[51]:
成分In [52]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
A[:, 1]
Out[52]:
添え字の範囲を指定して行列の部分を取り出すことができる。 In [53]:
A[0:2, 1:3]
Out[53]:
添え字のリストを指定することもできる。 In [54]:
A[:, [0, 2]]
Out[54]:
代入In [55]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
B = np.matrix(A)
B
Out[55]:
In [56]:
B = A.copy()
B
Out[56]:
スカラーとの和・差In [57]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
A + 1
Out[57]:
In [58]:
A - 1
Out[58]:
スカラーとの積・除算In [59]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
A*2.
Out[59]:
In [60]:
A/2.
Out[60]:
行列どうしの和・差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]:
In [62]:
B - A
Out[62]:
行列の成分どうしの積・除算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]:
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 の成分どうしで積を取りたい場合は、array に変換してやればよい。 In [65]:
np.matrix(np.array(A)*np.array(B))
Out[65]:
除算は普通にできる。 In [66]:
B/A
Out[66]:
行列の積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]:
dot() 関数を使ってもよい。 In [68]:
np.dot(A, B)
Out[68]:
次のような書き方もできる。 In [69]:
A.dot(B)
Out[69]:
array の場合、"*" の代わりに "@" を使えば行列の積になる。 In [70]:
a = np.array(A)
b = np.array(B)
a*b
Out[70]:
In [71]:
a@b
Out[71]:
dot() 関数はそのまま使える。 In [72]:
np.dot(a, b)
Out[72]:
行列の累乗In [73]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
A**2
Out[73]:
転置In [74]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
A.T
Out[74]:
トレースIn [75]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
np.trace(A)
Out[75]:
対角成分In [76]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]
])
np.diag(A)
Out[76]:
対角の和と取るとトレースになる。 In [77]:
np.sum(np.diag(A))
Out[77]:
行列式In [78]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]
])
np.linalg.det(A)
Out[78]:
逆行列In [79]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]
])
np.linalg.inv(A)
Out[79]:
In [80]:
A*np.linalg.inv(A)
Out[80]:
余因子行列を用いた逆行列の計算 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]:
連立方程式の解In [82]:
b = np.array([1., 2., 3.])
x = np.linalg.solve(A, b)
x
Out[82]:
In [83]:
A.dot(x)
Out[83]:
固有値・固有ベクトルIn [84]:
A = np.matrix([
[1., 2., 3.],
[4., 5., 6.],
[7., 8., 10.]
])
lam, P = np.linalg.eig(A)
P, lam
Out[84]:
対角化 In [85]:
np.linalg.inv(P)*A*P
Out[85]:
最適化非線形方程式を解くscipy.optimize の newton() で Newton 法により非線形方程式を解くことができる。 In [86]:
from scipy.optimize import newton
newton(lambda x: x**2 - 2., 1.)
Out[86]:
最小化scipy.optimize の fmin() 関数を最小化する変数を求めることができる。 In [87]:
from scipy.optimize import fmin
fmin(lambda x: abs(x**2 - 2.), 1.)
Out[87]:
カーブフィッティング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()
多項式補間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()
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()
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()
| |
PENGUINITIS |