Python: определение линейной части наклона

Ответ 1

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

• Как рассчитать вторую производную от шумных данных? Одним быстрым и простым методом, который может быть легко адаптирован к различным уровням шума, размерам данных и ожидаемым длинам линейного патча, является свертка данных с ядром свертки, равным второй производной гауссова. Регулируемая часть - это ширина ядра.

• Что означает "близко к нулю" в вашем контексте? Чтобы ответить на этот вопрос, вам нужно поэкспериментировать с вашими данными.

• Результаты этого метода могут быть использованы в качестве входных данных для описанного выше chi ^ 2-метода для идентификации областей-кандидатов в наборе данных.

Вот некоторые исходные тексты, которые помогут вам начать:

from matplotlib import pyplot as plt

import numpy as np

# create theoretical data
x_a = np.linspace(-8,0, 60)
y_a = np.sin(x_a)
x_b = np.linspace(0,4,30)[1:]
y_b = x_b[:]
x_c = np.linspace(4,6,15)[1:]
y_c = np.sin((x_c - 4)/4*np.pi)/np.pi*4. + 4
x_d = np.linspace(6,14,120)[1:]
y_d = np.zeros(len(x_d)) + 4 + (4/np.pi)

x = np.concatenate((x_a, x_b, x_c, x_d))
y = np.concatenate((y_a, y_b, y_c, y_d))


# make noisy data from theoretical data
y_n = y + np.random.normal(0, 0.27, len(x))

# create convolution kernel for calculating
# the smoothed second order derivative
smooth_width = 59
x1 = np.linspace(-3,3,smooth_width)
norm = np.sum(np.exp(-x1**2)) * (x1[1]-x1[0]) # ad hoc normalization
y1 = (4*x1**2 - 2) * np.exp(-x1**2) / smooth_width *8#norm*(x1[1]-x1[0])



# calculate second order deriv.
y_conv = np.convolve(y_n, y1, mode="same")

# plot data
plt.plot(x,y_conv, label = "second deriv")
plt.plot(x, y_n,"o", label = "noisy data")
plt.plot(x, y, label="theory")
plt.plot(x, x, "0.3", label = "linear data")
plt.hlines([0],-10, 20)
plt.axvspan(0,4, color="y", alpha=0.2)
plt.axvspan(6,14, color="y", alpha=0.2)
plt.axhspan(-1,1, color="b", alpha=0.2)
plt.vlines([0, 4, 6],-10, 10)
plt.xlim(-2.5,12)
plt.ylim(-2.5,6)
plt.legend(loc=0)
plt.show()

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

smooth_width - ширина ядра свертки. Чтобы отрегулировать уровень шума, измените значение 0.27 в случайном порядке. Нормально для разных значений. Обратите внимание, что этот метод не очень хорошо работает с границей пространства данных.

Как вы можете видеть, требование "близко к нулю" для второй производной (синяя линия) довольно хорошо выполняется для желтых частей, где данные являются линейными.

Ответ 2

Вы можете использовать алгоритм Ramer Douglas Peucker, чтобы упростить ваши данные до меньшего набора сегментов линии. Алгоритм позволяет указать epsilon, так что каждая точка данных будет не дальше, чем epsilon от некоторого сегмента линии. Наклон отрезков линии даст приблизительную оценку наклона кривой.

Существует Python реализация алгоритма RDP.

Ответ 3

Это всего лишь возможное решение, оно найдет отрезок прямой точек, у которого минимальное значение chi ^ 2 больше заданного минимума;

from matplotlib.pyplot import figure, show
from numpy import pi, sin, linspace, exp, polyfit
from matplotlib.mlab import stineman_interp

x = linspace(0,2*pi,20);
y = x + sin(x) + exp(-0.5*(x-2)**2);

num_points = len(x)

min_fit_length = 5

chi = 0

chi_min = 10000

i_best = 0
j_best = 0

for i in range(len(x) - min_fit_length):
    for j in range(i+min_fit_length, len(x)):

        coefs = polyfit(x[i:j],y[i:j],1)
        y_linear = x * coefs[0] + coefs[1]
        chi = 0
        for k in range(i,j):
            chi += ( y_linear[k] - y[k])**2

        if chi < chi_min:
            i_best = i
            j_best = j
            chi_min = chi
            print chi_min

coefs = polyfit(x[i_best:j_best],y[i_best:j_best],1)
y_linear = x[i_best:j_best] * coefs[0] + coefs[1]


fig = figure()
ax = fig.add_subplot(111)
ax.plot(x,y,'ro')
ax.plot(x[i_best:j_best],y_linear,'b-')


show()

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

Ответ 4

Если ваша "модель" данных состоит из данных, которые в основном подходят к прямой, с несколькими выбросами или фрагментами в конце, вы можете попробовать RANSAC.

Псевдокод (очень многословный, извините) был бы следующим:

choose a small threshold distance D

for N iterations:
    pick two random points from your data, a and b
    fit a straight line, L, to a and b
    count the inliers: data points within a distance D of the line L
    save the parameters of the line with the most inliers so far

estimate the final line using ALL the inliers of the best line