以上のアルゴリズムをPythonで実装してみた。
逆行列を求めるのには、逆行列を計算するのに便利なnumpyを使った。
import numpy as np
import sys
import matplotlib.pyplot as plt
class cubicSpline:
def __init__(self, x,y):
points = sorted(zip(x, y))
self.x, self.y = np.array(list(zip(*points)), dtype=float)
self.__initialize(self.x, self.y)
def __initialize(self,x,y):
xlen = len(x)
N = xlen - 1
if xlen != len(set(x)): raise ValueError("x must be different values")
matrix = np.zeros([4*N, 4*N])
Y = np.zeros([4*N])
equation = 0
for i in range(N):
for j in range(4):
matrix[equation, 4*i+j] = pow(x[i], j)
Y[equation] = y[i]
equation += 1
for i in range(N):
for j in range(4):
matrix[equation, 4*i+j] = pow(x[i+1], j)
Y[equation] = y[i+1]
equation += 1
for i in range(N-1):
for j in range(4):
matrix[equation, 4*i+j] = j*pow(x[i+1], j-1)
matrix[equation, 4*(i+1)+j] = -j*pow(x[i+1], j-1)
equation += 1
for i in range(N-1):
matrix[equation, 4*i+3] = 3*x[i+1]
matrix[equation, 4*i+2] = 1
matrix[equation, 4*(i+1)+3] = -3*x[i+1]
matrix[equation, 4*(i+1)+2] = -1
equation += 1
matrix[equation,3] = 3*x[0]
matrix[equation,2] = 1
equation += 1
matrix[4*N-1,4*N-1] = 3*x[N]
matrix[4*N-1,4*N-2] = 1
self.variables = np.dot(np.linalg.inv(matrix),Y)
def fit(self, x):
"""
引数xが該当する区間を見つけてきて補間後の値を返す
"""
xlen = len(self.x)
for index,j in enumerate(self.x):
if x < j:
index -= 1
break
if index == -1:
index += 1
elif index == xlen-1:
index -= 1
a3 = self.variables[4*index + 3]
a2 = self.variables[4*index + 2]
a1 = self.variables[4*index + 1]
a0 = self.variables[4*index + 0]
result = a3*pow(x,3) + a2*pow(x,2) + a1*x + a0
return result
if __name__=="__main__":
X = np.array([0,2,5,7,10])
y = np.array([2,10,3,2,4])
cubic = cubicSpline(X, y)
plt.scatter(X, y)
x = np.arange(-0.5,10.5,0.01)
plt.plot(x, list(map(cubic.fit,x)))
plt.show()
補足:ライブラリで用意している関数を使う
実際にPythonで3次スプライン補間するには、scipyで関数が用意されているので便利。
from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt
X = np.array([0,2,5,7,10])
y = np.array([2,10,3,2,4])
plt.scatter(X,y)
x = np.arange(0,10,0.01)
cubic = interp1d(X, y, kind="cubic")
plt.plot(x, cubic(x))
plt.show()
補足:回帰(最小二乗法)との比較
3次スプライン補間のような補間と最小二乗法による回帰との大きな違いは もともとの点群を通るかどうかだ。
したがってもともとの点自体が重要になるような場合は、補間を使う方が適していると言える。
しかし点の数が増えたり誤差や異常値が含まれている場合は回帰を使う方が良い。
from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt
X = np.array([0,2,5,6,7,8,10])
y = np.array([0,4,32,33,45,70,98])
plt.scatter(X,y)
x = np.arange(0,10,0.01)
cubic = interp1d(X, y, kind="cubic")
plt.plot(x, cubic(x) , c="r", label="cubic spline")
fit_params = np.polyfit(X, y, 2)
plt.plot(x, np.poly1d(fit_params)(x), color="g", label="polyfit")
plt.legend()
plt.show()
2次関数にフィッティングする線形回帰と3次スプライン補間の比較。
誤差を認める場合は線形回帰を使う方が良い。