なまもの備忘録

気になったことをつらつらと書いていきます

scipyによる3次スプライン補完

未知の関数に出会った時3次スプライン補完をしたいと思うのは人類の普遍的欲求です。 普遍的欲求はだいたいpythonのライブラリが満たしてくれるもので、今回のケースではscipyがその役割を担ってくれます。

scipyの3次スプライン補完についてはいくつか日本語解説ページがあるのですが、

qiita.com

ここで使われているkind='cubic'interp1d関数だと、スプライン補完の係数取得に少々非自明な処理をする必要が出てきます。 人間が3次スプライン補完をするときって、欲しいものは大抵関数そのものではなく係数なので、これは困りますね。

さて、scipyではinterp1d以外にもCubicSpline関数で3次スプライン補完を実行することができ、こちらを使えばこの問題は解決します。 使い方は以下

import numpy
from scipy import interpolate
import matplotlib.pyplot as plt

if __name__ == '__main__':

    x_coordinates = []
    y_coordinates = []
    for line in open('test.txt', 'r'):
        point = [float(coord.split()[0]) for coord in line.split(',')]
        x_coordinates.append(point[0])
        y_coordinates.append(point[1])
        
    spline_func = interpolate.CubicSpline(x_coordinates, y_coordinates)
    print(spline_func.c) # ここで係数を取得しています
    
    x_res = numpy.linspace(1, 5, 200)
    y_res = spline_func(x_res)
   
    plt.plot(x_coordinates, y_coordinates, 'o', x_res, y_res)
    plt.show()

CubicSplineの返り値であるPPolyオブジェクトのcを呼び出すことで、係数にアクセスしています。

テストデータとしては下記を使用しました。

1.5, 1.65
1.677, 5.64
2.3, 5.53
3.0, 4.354
3.8, 3.555

出力グラフ

f:id:purple-apple:20180226130144p:plain

標準出力(コメント部分は後から追記)

[[ 25.17215028  25.17215028  -3.16687764  -3.16687764] # 3次の係数
 [-52.99186318 -39.62545138   7.42129748   0.77085443] # 2次の係数
 [ 31.13331437  14.74004969  -5.32313819   0.41136815]    # 1次の係数 
 [  1.65         5.64         5.53         4.354     ]]    # 定数項

バージョン情報

$ python -V
Python 3.5.2
>>> import scipy as sp
>>> sp.version.full_version
'0.18.1'

参考

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicSpline.html

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.PPoly.html#scipy.interpolate.PPoly