[解決済み] scipy.interpolateが入力範囲を超えて外挿された結果を与えるようにするには？

質問

数学者の同僚が開発した手動の補間器を使ったプログラムを、scipyが提供する補間器を使うように移植しようとしています。scipy のインターポレータを使用するか、ラップして、古いインターポレータにできるだけ近い動作をさせたいと考えています。

この2つの関数の重要な違いは、オリジナルのインターポレーターでは - 入力値が入力範囲より上または下にある場合、オリジナルのインターポレーターは結果を外挿するということです。これをscipyのインターポレータで試すと ValueError . このプログラムを例として考えてみましょう。

import numpy as np
from scipy import interpolate

x = np.arange(0,10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y)

print f(9)
print f(11) # Causes ValueError, because it's greater than max(x)

クラッシュする代わりに、最後の線が単に線形外挿を行い、最初と最後の2点で定義された勾配を無限大まで継続するようにする賢明な方法はないでしょうか。

実際のソフトウェアでは、exp関数は使っていないことに注意してください。

どのように解決するのですか？

1. 定数外挿

を使用することができます。 interp 関数は、左右の値を範囲外の定数として外挿します。

>>> from scipy import interp, arange, exp
>>> x = arange(0,10)
>>> y = exp(-x/3.0)
>>> interp([9,10], x, y)
array([ 0.04978707,  0.04978707])

2. 線形（または他のカスタム）外挿

線形補外を行う補間関数のラッパーを書くことができます。例えば

from scipy.interpolate import interp1d
from scipy import arange, array, exp

def extrap1d(interpolator):
    xs = interpolator.x
    ys = interpolator.y

    def pointwise(x):
        if x < xs[0]:
            return ys[0]+(x-xs[0])*(ys[1]-ys[0])/(xs[1]-xs[0])
        elif x > xs[-1]:
            return ys[-1]+(x-xs[-1])*(ys[-1]-ys[-2])/(xs[-1]-xs[-2])
        else:
            return interpolator(x)

    def ufunclike(xs):
        return array(list(map(pointwise, array(xs))))

    return ufunclike

extrap1d は、補間関数を受け取り、外挿もできる関数を返します。そして、このように使うことができる。

x = arange(0,10)
y = exp(-x/3.0)
f_i = interp1d(x, y)
f_x = extrap1d(f_i)

print f_x([9,10])

出力します。

[ 0.04978707  0.03009069]