In [53]:
# 科学計算モジュール Scipy
"""
補間 スプライン補間 interplate linalg
特異値分解 LU分解 コレスキー分解
数値積分 微分方程式 integrate
最適化 二分法 ブレント法 ニュートン法 optimize
高速フーリエ変換 信号処理 画像処理
"""
Out[53]:
In [54]:
import numpy as np
import numpy.random as random
import scipy as sp
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline
%precision 3
Out[54]:
In [13]:
# 補間
# 0から10 を11個のノードで等分する
x = np.linspace(0, 10, num=11, endpoint=True)
x
"""
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
"""
y = np.cos(x)
y
"""
array([ 1. , 0.54 , -0.416, -0.99 , -0.654, 0.284, 0.96 , 0.754,
-0.146, -0.911, -0.839])
"""
plt.plot(x, y, marker='o', linestyle='', color='blue')
plt.grid(True)
plt.show()
In [17]:
# 線形補間 sp.interplolate.interp1d
# ※ピーいちデー なので注意
import scipy as sp
import scipy.interpolate # 明示的にinterplateをインポートしないと動かない
# sp.interplolate.interp1d インスタンスを生成する
f = sp.interpolate.interp1d(x, y, 'linear')
plt.plot(x, f(x), '-')
plt.grid(True)
# 一次式で繋がれているので fに任意のxを入れれば値が帰ってくる。
print(f(2.2))
print(f(7.5))
print(f(9.3))
In [18]:
# dataプロット
plt.plot(x, y, 'o')
#リニア１次補間
f_linear = sp.interpolate.interp1d(x, y, 'linear')
# スプライン２次補間
f_quadratic = sp.interpolate.interp1d(x, y, 'quadratic')
# スプライン３次補間
f_cubic = sp.interpolate.interp1d(x, y, 'cubic')
# スプライン描写用
x_new = np.linspace(0,10,101)
plt.plot(x , f_linear(x),linestyle='-')
plt.plot(x_new, f_quadratic(x_new),linestyle='--')
plt.plot(x_new, f_cubic(x_new),linestyle="dashdot")
plt.grid(True)
plt.legend(['data', 'linear', 'quadratic', 'cubic'])
plt.show()
In [28]:
# 線形大数
# 特異値分解 Singular value decomposition SVD
"""
Ax = λA で
Aが正方行列のとき xは固有ベクトル λは固有値 になる。
Aが正方行列でないときも特異値分解をすればでxとλが算出できる。
A = UΣV*
A:行列(m, n)
U:AA'の固有ベクトルを列ベクトルとして並べた行列
（A'はAの共役転置行列）
Σ:特異値を対角に並べた行列
V:A'A の固有ベクトルを列ベクトルとして並べた行列
A'A の固有ベクトル:min(m,n)
正の固有値をσi^2としたとき σiを特異値と呼ぶ
"""
A = np.array([[1,2,3,4,5],[6,7,8,9,10]])
#特異値分解の関数 linalg.svd
U, s, Vs = sp.linalg.svd(A)
print(U)
print(s)
print(Vs)
"""
[[-0.37 -0.929]
[-0.929 0.37 ]]
[19.538 1.81 ]
[[-0.304 -0.371 -0.437 -0.504 -0.57 ]
[ 0.712 0.403 0.094 -0.215 -0.524]
[-0.374 -0.008 0.862 -0.206 -0.274]
[-0.365 0.371 -0.159 0.665 -0.512]
[-0.357 0.75 -0.179 -0.464 0.251]]
"""
m, n = A.shape
print(m, n) #2 5
S = sp.linalg.diagsvd(s, m, n)
print(S)
"""
[[19.538 0. 0. 0. 0. ]
[ 0. 1.81 0. 0. 0. ]]
"""
print('U.S.V* = \n', U @ S @Vs)
"""
U.S.V* =
[[ 1. 2. 3. 4. 5.]
[ 6. 7. 8. 9. 10.]]
"""
In [42]:
# LU分解
"""
A：正方行列
Ax = b を解く代わりに PLUx =b を解く。
置換行列P,
対角成分が全て１の下三角行列L
上三角行列UをA = PLUとなるようにおく
"""
A = np.identity(5) # 単位行列を生成する
A[0,:] = 1
A[:,0] = 1
A[0,0] = 5
"""
[[5. 1. 1. 1. 1.]
[1. 1. 0. 0. 0.]
[1. 0. 1. 0. 0.]
[1. 0. 0. 1. 0.]
[1. 0. 0. 0. 1.]]
"""
b = np.ones(5)
print(b) # [1. 1. 1. 1. 1.]
#正方行列をLU分解する
(LU, piv) = sp.linalg.lu_factor(A)
L = np.identity(5) + np.tril(LU, -1)
U = np.triu(LU)
P = np.identity(5)[piv]
# 解を求める
x = sp.linalg.lu_solve((LU,piv),b)
x
"""
array([-3., 4., 4., 4., 4.])
"""
Out[42]:
In [50]:
# コレスキー分解
A = np.array([[7,-1,0,1],
[-1,9,-2,2],
[0,-2,8,-3],
[1,2,-3,10]])
b = np.array([5,20,0,20])
L = sp.linalg.cholesky(A)
print(L)
"""
[[ 2.646 -0.378 0. 0.378]
[ 0. 2.976 -0.672 0.72 ]
[ 0. 0. 2.747 -0.916]
[ 0. 0. 0. 2.915]]
"""
t = sp.linalg.solve(L.T.conj(),b)
print(t)
"""
[1.89 6.96 1.702 5.431]
"""
x = sp.linalg.solve(L, t)
print(x)
"""
[0.758 2.168 1.241 1.863]
"""
# 確認
np.dot(A, x)
"""
array([5.000e+00, 2.000e+01, 1.776e-15, 2.000e+01])
bと同値になっているのでOK
"""
In [52]:
"""
QR分解
非負値行列因子分解：Non-negative Matrix Factorization NMF
リコメンデーションに使用する
"""
Out[52]:
In [64]:
# NMF
from sklearn.decomposition import NMF
X = np.array([[1,1,1],
[2,2,2],
[3,3,3],
[4,4,4]])
model = NMF(n_components=2, init='random', random_state=0)
W = model.fit_transform(X)
H = model.components_
print(W)
print(H)
"""
[[0.425 0.222]
[0.698 0.537]
[0.039 1.434]
[2.377 0.463]]
[[1.281 1.281 1.282]
[2.058 2.058 2.058]]
"""
# 確認
np.dot(W, H)
"""
array([[1., 1., 1.],
[2., 2., 2.],
[3., 3., 3.],
[4., 4., 4.]])
"""
Out[64]:
In [73]:
# 積分計算
from scipy import integrate
import math
def calc(x):
return 4/(1+x**2)
result = integrate.quad(calc, 0, 1) # 積分範囲0 ,1
print(result[0]) # 結果：3.1415926535897936
print(result[1]) # 推定誤差：3.4878684980086326e-14
# 無名関数も使用可能
result = integrate.quad(lambda x: 4/(1+x**2), 0, 1)
In [75]:
# sin関数を求める
result = integrate.quad(np.sin, 0, 1)
print(result)
In [79]:
# ２重積分
def I(n):
return integrate.dblquad(lambda t, x:np.exp(-x*t)/t**n, 0, np.inf, lambda x: 1, lambda x: np.inf)
print(I(1))
print(I(2))
print(I(3))
print(I(4))
print(I(5))
In [95]:
# 微分方程式
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# ローレンツ方程式
def lorenz_eq(v,t,p,r,b):
return [-p*v[0] + p*v[1], -v[0]*v[2] + r*v[0]-v[1], v[0]*v[1] - b*v[2]]
# パラメータの設定
p = 10
r = 28
b = 8/3
v0 = [0.1, 0.1, 0.1]
t = np.arange(0, 100, 0.01)
# 関数の呼び出し
v = odeint(lorenz_eq, v0, t, args=(p,r, b))
# 可視化
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(v[:,0], v[:,1], v[:,2])
# ラベル
plt.title('Lorenz equation')
plt.grid(True)
In [107]:
# 最適化
from scipy.optimize import fsolve
# ２次関数の最適化
def f(x):
y = 2 * x**2 +2*x -10
return y
#可視化
x =np.linspace(-4,4,101)
plt.plot(x, f(x))
plt.grid(True)
"""
解は-3と2付近にある
"""
# 解く
x = fsolve(f, -3)# 関数名, xの初期値
print(x) #[-2.791]
x = fsolve(f, 2)# 関数名, xの初期値
print(x) # [1.791]
In [114]:
# 最適化問題を解く
from scipy.optimize import minimize
# 目的となる関数
def objective_func(x):
x1 = x[0]
x2 = x[1]
x3 = x[2]
x4 = x[3]
return x1 * x4 * (x1+x2+x3) + x3
# 制約条件１
def constraint_1(x):
return x[0]*x[1]*x[2]*x[3] - 25.0
# 成約条件2
def constraint_2(x):
sum_sq = 40
for i in range(4):
suj_sq = sum_sq - x[i]**2
return sum_sq
# 初期値
x0 = [1,5,5,1]
print(objective_func(x0)) # 16
# パラメータ
b = (1.0, 5.0)
bnds = (b,b,b,b)
con1 = {'type':'ineq', 'fun':constraint_1}
con2 = {'type':'ineq', 'fun':constraint_2}
cons = [con1, con2]
solve = minimize(objective_func,x0, method='SLSQP', bounds=bnds, constraints=cons)
print(solve)
In [116]:
print('Y:', sol.fun)
print('X:', sol.x)
"""
Y: 16.0
X: [1. 5. 5. 1.]
"""