Posted at

MスプラインとIスプラインをNumpyで

monotone(単調増加 or 単調減少)なスプライン、それがIスプライン。

Iスプラインを作るために使う、負にならないスプライン、それがMスプライン。

Iスプラインは、monotoneなはずのデータの回帰分析やスムージングに使えそうに見える。しかしどちらも、どうも世の中でほとんど使われていないようなので、数式をコードに起こすのに苦労した。

import numpy as np

import matplotlib.pyplot as plt

def main():
K = 3 # 次数
L = 4 # ノットの数
x_min = 0.1
x_max = 1.0
ts = np.concatenate([
np.full(K, x_min),
np.linspace(x_min, x_max, L + 2)[1:-1], # ノットをuniformに配置。non-uniformでもいい
np.full(K + 1, x_max)
])

xs = np.linspace(x_min, x_max - 1e-5, 100)

def function_m(k, i, where=np.ones_like(xs, np.bool)):
lxs = xs[where]
ys = np.zeros_like(lxs)
if k == 1:
hit = (ts[i] <= lxs) & (lxs < ts[i + 1])
ys[hit] = 1. / (ts[i + 1] - ts[i])
elif ts[i + k] - ts[i] != 0.:
_a = (k - 1) * (ts[i + k] - ts[i])
_b = lxs - ts[i]
_c = function_m(k - 1, i, where)
_d = ts[i + k] - lxs
_e = function_m(k - 1, i + 1, where)
ys = k * (_b * _c + _d * _e) / _a
return ys

for i in np.arange(K + L):
plt.plot(xs, function_m(K, i), label=str(i))
plt.legend(bbox_to_anchor=(1,1), loc='upper right')
plt.title('M-spline')
plt.show()

def function_i(i):
hit = (ts[:-1,np.newaxis] <= xs) & (xs < ts[1:,np.newaxis])
js = np.argmax(hit, axis=0)
def f_a(m):
mask = js >= m
_a = K + 1
_b = ts[m + K + 1] - ts[m]
_c = function_m(K + 1, m, where=mask)
ys = np.zeros_like(xs)
ys[mask] = _b * _c / _a
return ys
ufunc_a = np.vectorize(f_a, signature='()->(n)')
yss = ufunc_a(np.arange(i, K + L))
return yss.sum(axis=0)

for i in np.arange(K + L):
plt.plot(xs, function_i(i), label=str(i))
plt.legend(bbox_to_anchor=(1,1), loc='upper right')
plt.title('I-spline')
plt.show()

if __name__=='__main__':
main()

参考:

Wikipedia I-spline

spline

INFERENCE USING SHAPE-RESTRICTED REGRESSION SPLINES