Сплайн-представление с scipy.interpolate: плохая интерполяция для низкоамплитудных, быстро колеблющихся функций

Question

Мне нужно (численно) вычислить первую и вторую производную функции, для которой я пытался использовать и splrep, и UnivariateSpline для создания сплайнов с целью интерполяции функции для получения производных.

Однако, похоже, что в самом представлении сплайна есть внутренняя проблема для функций, величина которых порядка 10 ^ -1 или ниже и (быстро) колеблются.

В качестве примера рассмотрим следующий код для создания сплайн-представления функции синуса на интервале (0,6 * pi) (поэтому функция колеблется только три раза):

import scipy
from scipy import interpolate
import numpy
from numpy import linspace
import math
from math import sin

k = linspace(0, 6.*pi, num=10000) #interval (0,6*pi) in 10'000 steps
y=[]
A = 1.e0 # Amplitude of sine function

for i in range(len(k)):

  y.append(A*sin(k[i]))

tck =interpolate.UnivariateSpline(x, y, w=None, bbox=[None, None], k=5, s=2)
M=tck(k)

Ниже приведены результаты для M для A = 1.e0 и A = 1.e-2

http://i.imgur.com/uEIxq.png Амплитуда = 1

http://i.imgur.com/zFfK0.png Амплитуда = 1/100

Очевидно, что интерполированная функция, созданная сплайнами, совершенно неверна! 2-й график даже не колеблется с правильной частотой.

У кого-нибудь есть понимание этой проблемы? Или знаете другой способ создания сплайнов внутри numpy / scipy?

Ура, Рори

Answer 1

Я предполагаю, что ваша проблема связана с псевдонимами.

Что такое x в вашем примере?

Если значения x, которые вы интерполируете, расположены менее близко, чем ваши исходные точки, вы по своей природе потеряете информацию о частоте. Это полностью независимо от любого типа интерполяции. Это присуще понижающей выборке.

Не берите в голову вышеупомянутый бит о псевдонимах. Это не применимо в этом случае (хотя я все еще не знаю, что такое x в вашем примере ...

Я только что понял, что вы оцениваете свои точки на исходных входных точках, когда вы используете ненулевой коэффициент сглаживания (s).

По определению, сглаживание не будет точно соответствовать данным. Попробуйте вместо этого ввести s=0.

В качестве быстрого примера:

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

x = np.linspace(0, 6.*np.pi, num=100) #interval (0,6*pi) in 10'000 steps
A = 1.e-4 # Amplitude of sine function
y = A*np.sin(x)

fig, axes = plt.subplots(nrows=2)
for ax, s, title in zip(axes, [2, 0], ['With', 'Without']):
    yinterp = interpolate.UnivariateSpline(x, y, s=s)(x)
    ax.plot(x, yinterp, label='Interpolated')
    ax.plot(x, y, 'bo',label='Original')
    ax.legend()
    ax.set_title(title + ' Smoothing')

plt.show()

Причина того, что вы четко видите эффекты сглаживания с низкой амплитудой, заключается в способе определения коэффициента сглаживания. См. Документацию для scipy.interpolate.UnivariateSpline для более подробной информации.

Даже с более высокой амплитудой интерполированные данные не будут соответствовать исходным данным, если вы используете сглаживание.

Например, если мы просто изменим амплитуду (A) на 1.0 в примере кода выше, мы все равно увидим эффекты сглаживания ...

Answer 2

Проблема заключается в выборе подходящих значений для параметра s.Его значения зависят от масштабирования данных.

Внимательно читая документацию, можно сделать вывод, что параметр следует выбирать в пределах s = len(y) * np.var(y), т. Е. Число точек данных * дисперсия.Например, s = 0.05 * len(y) * np.var(y) дает сглаживающий сплайн, который не зависит от масштабирования данных или количества точек данных.

EDIT : разумные значения для s зависят, конечнотакже об уровне шума в данных.Документы, похоже, рекомендуют выбрать s в диапазоне (m - sqrt(2*m)) * std**2 <= s <= (m + sqrt(2*m)) * std**2, где std - стандартное отклонение, связанное с "шумом", который вы хотите сгладить.