Как преобразовать подгонку сплайна в кусочную функцию?

Question

Допустим, у меня есть

import numpy as np
from scipy.interpolate import UnivariateSpline

# "true" data; I don't know this function
x = np.linspace(0, 100, 1000)
d = np.sin(x * 0.5) + 2 + np.cos(x * 0.1)

# sample data; that's what I actually measured
x_sample = x[::20]
d_sample = d[::20]

# fit spline
s = UnivariateSpline(x_sample, d_sample, k=3, s=0.005)

plt.plot(x, d)
plt.plot(x_sample, d_sample, 'o')
plt.plot(x, s(x))
plt.show()

Я получаю

То, что я хотел бы сейчас иметь, это функции между всеми оранжевые точки, что-то вроде

knots = s.get_knots()
f0 = <some expression> for knots[0] <= x < knots[1]
f1 = <some expression> for knots[1] <= x < knots[2]
...

. Таким образом, fi следует выбирать таким образом, чтобы он воспроизводил форму посадки сплайна.

Я нашел сообщение здесь , но созданный там сплайн кажется неверным для приведенного выше примера, и это также не совсем то, что мне нужно, поскольку он не возвращает выражения.

Как я могу превратить сплайн в кусочную функцию? Есть ли (простой) способ express каждого интервала, например, в виде полинома?

Answer 1

Короткий ответ: если вас интересуют коэффициенты полинома в стандартной степени мощности, вам лучше использовать CubicSpline (и см. это обсуждение ):

cu = scipy.interpolate.CubicSpline(x_sample, d_sample)

plt.plot(x_sample, y_sample, 'ko')
for i in range(len(cu.x)-1):
    xs = np.linspace(cu.x[i], cu.x[i+1], 100)
    plt.plot(xs, np.polyval(cu.c[:,i], xs - cu.x[i]))

И чтобы ответить на ваш вопрос, вы можете вместо этого создать кусочную функцию, используя numpy.piecewise, точки останова в cu.x и коэффициенты в cu.c, либо либо непосредственно кодируйте полиномиальные выражения самостоятельно, либо используйте numpy.polyval. Например,

cu.c[:,0]  # coeffs for 0th segment
# array([-0.01316353, -0.02680068,  0.51629024,  3.        ])

# equivalent ways to code polynomial for this segment
f0 = lambda x: cu.c[0,0]*(x-x[0])**3 + cu.c[1,0]*(x-x[0])**2 + cu.c[2,0]*(x-x[0]) + cu.c[3,0]
f0 = lambda x: [cu.c[i,0]*(x-x[0])**(3-i) for i in range(4)]

# ... or getting values directly from x's
y0 = np.polyval(cu.c[:,0], xs - cu.x[0])

БОЛЬШОЙ ОТВЕТ:

Здесь есть несколько точек возможной путаницы:

• UnivariateSpline подходит B-сплайн , поэтому коэффициенты не совпадают со стандартным полиномиальным базисом
• . Для преобразования из B-сплайна мы можем использовать PPoly.from_spline, но, к сожалению UnivariateSpline возвращает усеченный список узлов и коэффициентов, которые не будут играть с этой функцией. Мы можем решить эту проблему, обратившись к внутренним данным объекта сплайна, что является небольшим табу.
• Кроме того, матрица коэффициентов c (от UnivariateSpline или CubicSpline) имеет обратную степень заказать и предполагает, что вы «центрируете» себя, например, коэффициент в c[k,i] принадлежит c[k,i]*(x-x[i])^(3-k).

Учитывая ваши настройки, обратите внимание, что если вместо использования оболочки UnivariateSpline, мы напрямую подходим с splrep и без сглаживания (s=0), мы можем захватить кортеж tck (knots-coefficients coele) и отправить его в функцию PPoly.from_spline и получить нужные нам коэффициенты:

tck = scipy.interpolate.splrep(x_sample, d_sample, s=0)
tck
# (array([0.        , 0.        , 0.        , 0.        , 2.68456376,
#        4.02684564, 5.36912752, 6.7114094 , 9.39597315, 9.39597315,
#        9.39597315, 9.39597315]),
# array([3.        , 3.46200469, 4.05843704, 3.89649312, 3.33792889,
#        2.29435138, 1.65015175, 1.59021688, 0.        , 0.        ,
#        0.        , 0.        ]),
# 3)

p = scipy.interpolate.PPoly.from_spline(tck)
p.x.shape  # breakpoints in unexpected shape
# (12,)

p.c.shape  # polynomial coeffs in unexpected shape
# (4, 11)

Обратите внимание на странные повторяющиеся точки останова в tck и снова в p.x: это вещь FITPACK (алгоритм, выполняющий все это).

Если мы пытаемся отправить кортеж tck из UnivariateSpline с (s.get_knots(), s.get_coeffs(), 3), мы пропускаем эти повторы, поэтому from_spline не работает. Извлеките источник , хотя кажется, что полный вектор хранится в self._data, поэтому мы можем сделать

s = scipy.interpolate.UnivariateSpline(x_sample, d_sample, s=0)
tck = (s._data[8], s._data[9], 3)
p = scipy.interpolate.PPoly.from_spline(tck)

и получить то же, что и раньше. Чтобы проверить эти коэффициенты, работайте:

plt.plot(x_sample, d_sample, 'o')

for i in range(len(p.x)-1):
    xs = np.linspace(p.x[i], p.x[i+1], 10)
    plt.plot(xs, np.polyval(p.c[:,i], xs - p.x[i]))

Примечание numpy.polyval требует обратного порядка для коэффициентов, чтобы мы могли передать p.c как есть.

Answer 2

Должно быть в состоянии использовать кусочную функцию, как вы описали, что-то вроде:

import numpy as np
from scipy.interpolate import UnivariateSpline

# "true" data; I don't know this function
x = np.linspace(0, 100, 1000)
d = np.sin(x * 0.5) + 2 + np.cos(x * 0.1)

# sample data; that's what I actually measured
x_sample = x[::20]
d_sample = d[::20]

# fit spline
s = UnivariateSpline(x_sample, d_sample, k=3, s=0.005)

plt.plot(x, d)
plt.plot(x_sample, d_sample, 'o')
plt.plot(x, s(x))

knots = s.get_knots()

conditions = [x < knots[0], (x >= knots[0]) * (x < knots[1]), (x >= knots[1]) * (x < knots[10]), x >= knots[10]]
# need one function for each condition
fns = [0, lambda x :-x, lambda x: 0.01*x**2, lambda x: 0.2*x**0.5]
y = np.piecewise(x, conditions, fns)

plt.plot(x, y)
plt.show()

Я выбрал несколько случайных условий и функций - я уверен, что вы могли бы найти что-то более подходящее!