scipy.interpolate.UnivariateSpline не сглаживает независимо от параметров

Question

У меня проблемы с получением scipy.interpolate.UnivariateSpline для использования любого сглаживания при интерполяции. Основываясь на странице функции , а также на некоторых предыдущих постах , я считаю, что она должна обеспечивать сглаживание с параметром s.

Вот мой код:

# Imports
import scipy
import pylab

# Set up and plot actual data
x = [0, 5024.2059124920379, 7933.1645067836089, 7990.4664106277542, 9879.9717114947653, 13738.60563208926, 15113.277958924193]
y = [0.0, 3072.5653360000988, 5477.2689107965398, 5851.6866463790966, 6056.3852496014106, 7895.2332350173638, 9154.2956175610598]
pylab.plot(x, y, "o", label="Actual")

# Plot estimates using splines with a range of degrees
for k in range(1, 4):
    mySpline = scipy.interpolate.UnivariateSpline(x=x, y=y, k=k, s=2)
    xi = range(0, 15100, 20)
    yi = mySpline(xi)
    pylab.plot(xi, yi, label="Predicted k=%d" % k)

# Show the plot
pylab.grid(True)
pylab.xticks(rotation=45)
pylab.legend( loc="lower right" )
pylab.show()

Вот результат:

Я пробовал это с диапазоном значений s (0,01, 0,1, 1, 2, 5, 50), а также с явными весами, установленными на одно и то же (1,0) или случайным образом. Я все еще не могу получить сглаживание, и количество узлов всегда совпадает с количеством точек данных. В частности, я ищу сглаживающие выбросы, подобные этой 4-й точке (7990.4664106277542, 5851.6866463790966).

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

Answer 1

Краткий ответ: вам нужно более тщательно выбрать значение для s.

В документации для UnivariateSpline указано, что:

Positive smoothing factor used to choose the number of knots. Number of 
knots will be increased until the     smoothing condition is satisfied:
sum((w[i]*(y[i]-s(x[i])))**2,axis=0) <= s

Из этогоМожно сделать вывод, что «разумные» значения для сглаживания, если вы не передаете явные веса, составляют около s = m * v, где m - это число точек данных, а v - дисперсия данных.В этом случае s_good ~ 5e7.

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

Answer 2

@ Ответ Жени о ручной установке узлов между точками данных был слишком грубым, чтобы давать хорошие результаты при наличии зашумленных данных без избирательного подхода к применению этого метода.Однако, вдохновленный его / ее предложением, я добился успеха с Mean-Shift кластеризацией из пакета scikit-learn.Он выполняет автоматическое определение количества кластеров и, по-видимому, выполняет довольно хорошую работу по сглаживанию (на самом деле очень сглаживает).

Answer 3

Хотя я не знаю ни одной библиотеки, которая сделает это за вас, я бы попробовал немного больше самодельного подхода: я бы начал с создания сплайна с узлами между точками исходных данныхx и y.В вашем конкретном примере наличие единого узла между 4-м и 5-м точками должно помочь, так как это удалит огромную производную в районе x=8000.

Answer 4

У меня были проблемы с запуском ответа BigChef, вот вариант, который работает на python 3.6:

# Imports
import pylab
import scipy
import sklearn.cluster

# Set up original data - note that it's monotonically increasing by X value!
data = {}
data['original'] = {}
data['original']['x'] = [0, 5024.2059124920379, 7933.1645067836089, 7990.4664106277542, 9879.9717114947653, 13738.60563208926, 15113.277958924193]
data['original']['y'] = [0.0, 3072.5653360000988, 5477.2689107965398, 5851.6866463790966, 6056.3852496014106, 7895.2332350173638, 9154.2956175610598]

# Cluster data, sort it and and save
import numpy
inputNumpy = numpy.array([[data['original']['x'][i], data['original']['y'][i]] for i in range(0, len(data['original']['x']))])
meanShift = sklearn.cluster.MeanShift()
meanShift.fit(inputNumpy)
clusteredData = [[pair[0], pair[1]] for pair in meanShift.cluster_centers_]

clusteredData.sort(key=lambda li: li[0])
data['clustered'] = {}
data['clustered']['x'] = [pair[0] for pair in clusteredData]
data['clustered']['y'] = [pair[1] for pair in clusteredData]

# Build a spline using the clustered data and predict
mySpline = scipy.interpolate.UnivariateSpline(x=data['clustered']['x'], y=data['clustered']['y'], k=1)
xi = range(0, int(round(max(data['original']['x']), -3)) + 3000, 20)
yi = mySpline(xi)

# Plot the datapoints
pylab.plot(data['clustered']['x'], data['clustered']['y'], "D", label="Datapoints (%s)" % 'clustered')
pylab.plot(xi, yi, label="Predicted (%s)" %  'clustered')
pylab.plot(data['original']['x'], data['original']['y'], "o", label="Datapoints (%s)" % 'original')

# Show the plot
pylab.grid(True)
pylab.xticks(rotation=45)
pylab.show()