Python: Dokładne wskazanie liniowej części stoku

Question

Mam kilka działek, które wyglądają jak następuje:

Zastanawiam się, jakiego rodzaju metod nie może być za znalezienie nachylenie pomiędzy około 5,5 i 8 dla osi x. Tam, gdzie jest kilka takich stron, zastanawiam się, czy istnieje sposób automatycznego znalezienia wartości nachylenia.

Wszelkie sugestie?

Myślę o ployfit() lub regresji liniowej. Problem polega na tym, że nie jestem pewien, jak znaleźć wartości automatycznie.

Answer 1

Jeśli "model" danych składa się z danych, które w większości pasują do linii prostej, z kilkoma odstającymi lub bujnymi bitami na końcu, można wypróbować algorytm RANSAC.

(bardzo rozwlekły, przepraszam) pseudokod tutaj byłoby:

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 

Answer 2

Można użyć Ramer Douglas Peucker algorithm uprościć swoje dane do mniejszego zestawu odcinków. Algorytm pozwala określić wartość epsilon tak, aby każdy punkt danych znajdował się dalej niż epsilon z pewnego segmentu linii. Nachylenie segmentów linii dałoby przybliżoną ocenę nachylenia krzywej.

Jest Python implementation of the RDP algorithm here.

Answer 3

To jest tylko jednym z możliwych rozwiązań, znajdzie segment prostą punktów, który spełnia minimalne wartości chi^2, który jest dłuższy niż ustawionego minimum;

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() 

widzę to coraz większy problem dla zestawów danych chociaż ...

Answer 4

Ogólny sposób na znalezienie elementów liniowych w zbiorach danych jest obliczyć drugą pochodną funkcji, i zobacz, gdzie jest (blisko) zero. Jest kilka rzeczy do rozważenia na drodze do rozwiązania:

• Jak obliczyć drugą pochodną hałaśliwych danych? Jedną szybką i prostą metodą, którą można łatwo dostosować do różnych poziomów hałasu, rozmiarów zbioru danych i spodziewanych długości łaty liniowej, jest połączenie danych z jądrem splotu równym drugiej pochodnej Gaussa. Regulowana część to szerokość jądra.

• Co oznacza "bliski zeru" w twoim kontekście? Aby odpowiedzieć na to pytanie, musisz poeksperymentować z danymi.

• Wyniki tej metody mogą być wykorzystane jako dane wejściowe do opisanej powyżej metody chi^2, w celu zidentyfikowania regionów kandydujących w zbiorze danych.

Oto niektóre kod źródłowy, który będzie Ci zacząć:

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() 

Jest to wynikiem:

smooth_width jest szerokość jądra splotu. Aby dostosować wielkość szumu, zmień wartość 0.27 w random.normal na różne wartości. I proszę zauważyć, że ta metoda nie działa dobrze w pobliżu granicy przestrzeni danych.

Jak widać, wymóg "zbliżenia do zera" dla drugiej pochodnej (niebieska linia) utrzymuje się całkiem dobrze dla żółtych części, gdzie dane są liniowe.