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