Итеративно подобранная полиномиальная кривая

Question

Я хочу, чтобы итеративно подогнать кривую к данным в python следующим способом:

1. Подогнать полиномиальную кривую (или любой нелинейный подход)
2. Сбросить значения> 2стандартное отклонение от среднего значения кривой
3. повторять шаги 1 и 2 до тех пор, пока все значения не окажутся в пределах доверительного интервала кривой

Я могу подогнать полиномиальную кривую следующим образом:

vals = array([0.00441025, 0.0049001 , 0.01041189, 0.47368389, 0.34841961,
       0.3487533 , 0.35067096, 0.31142986, 0.3268407 , 0.38099566,
       0.3933048 , 0.3479948 , 0.02359819, 0.36329588, 0.42535543,
       0.01308297, 0.53873956, 0.6511364 , 0.61865282, 0.64750302,
       0.6630047 , 0.66744816, 0.71759617, 0.05965622, 0.71335208,
       0.71992683, 0.61635697, 0.12985441, 0.73410642, 0.77318621,
       0.75675988, 0.03003641, 0.77527201, 0.78673995, 0.05049178,
       0.55139476, 0.02665514, 0.61664748, 0.81121749, 0.05521697,
       0.63404375, 0.32649395, 0.36828268, 0.68981099, 0.02874863,
       0.61574739])
x_values = np.linspace(0, 1, len(vals))
poly_degree = 3

coeffs = np.polyfit(x_values, vals, poly_degree)
poly_eqn = np.poly1d(coeffs)
y_hat = poly_eqn(x_values)

Как мне выполнить шаги 2 и 3?

Answer 1

Поскольку точки удаления слишком далеки от ожидаемого решения, вы, вероятно, ищете RANSAC (RANdom SAmple Consensus), который соответствует кривой (или любой другой функции) для данные в определенных границах, как ваш случай с 2 ​​* STD.

Вы можете использовать оценку scikit-learn RANSAC , которая хорошо согласована с включенными регрессорами, такими как LinearRegression . Для вашего полиномиального случая вам нужно определить свой собственный класс регрессии:

from sklearn.metrics import mean_squared_error
class PolynomialRegression(object):
    def __init__(self, degree=3, coeffs=None):
        self.degree = degree
        self.coeffs = coeffs

    def fit(self, X, y):
        self.coeffs = np.polyfit(X.ravel(), y, self.degree)

    def get_params(self, deep=False):
        return {'coeffs': self.coeffs}

    def set_params(self, coeffs=None, random_state=None):
        self.coeffs = coeffs

    def predict(self, X):
        poly_eqn = np.poly1d(self.coeffs)
        y_hat = poly_eqn(X.ravel())
        return y_hat

    def score(self, X, y):
        return mean_squared_error(y, self.predict(X))

и тогда вы можете использовать RANSAC

from sklearn.linear_model import RANSACRegressor
ransac = RANSACRegressor(PolynomialRegression(degree=poly_degree),
                         residual_threshold=2 * np.std(y_vals),
                         random_state=0)
ransac.fit(np.expand_dims(x_vals, axis=1), y_vals)
inlier_mask = ransac.inlier_mask_

Обратите внимание, что переменная X преобразуется в массив 2d, как того требует реализация sklearn RANSAC, а в нашем пользовательском классе сглаживается из-за того, что функция nify polyfit работает с массивом 1d.

y_hat = ransac.predict(np.expand_dims(x_vals, axis=1))
plt.plot(x_vals, y_vals, 'bx', label='input samples')
plt.plot(x_vals[inlier_mask], y_vals[inlier_mask], 'go', label='inliers (2*STD)')
plt.plot(x_vals, y_hat, 'r-', label='estimated curve')

Более того, играя с полиномиальным порядком и остаточным расстоянием, я получил следующие результаты со степенью = 4 и диапазоном 1 * STD

* * 1030

Другой вариант - использовать регрессор более высокого порядка, например Гауссовский процесс

from sklearn.gaussian_process import GaussianProcessRegressor
ransac = RANSACRegressor(GaussianProcessRegressor(),
                         residual_threshold=np.std(y_vals))

Говоря об обобщении в DataFrame, вам просто нужно установить, что все столбцы, кроме одного, являются элементами, а оставшийся - это вывод, как здесь:

import pandas as pd
df = pd.DataFrame(np.array([x_vals, y_vals]).T)
ransac.fit(df[df.columns[:-1]], df[df.columns[-1]])
y_hat = ransac.predict(df[df.columns[:-1]])

Answer 2

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

с учетом сказанного, вот как ответить на ваш вопрос:

начните с поиска в библиотеках и получения некоторых данных:

import matplotlib.pyplot as plt
import numpy as np

Y = np.array([
    0.00441025, 0.0049001 , 0.01041189, 0.47368389, 0.34841961,
    0.3487533 , 0.35067096, 0.31142986, 0.3268407 , 0.38099566,
    0.3933048 , 0.3479948 , 0.02359819, 0.36329588, 0.42535543,
    0.01308297, 0.53873956, 0.6511364 , 0.61865282, 0.64750302,
    0.6630047 , 0.66744816, 0.71759617, 0.05965622, 0.71335208,
    0.71992683, 0.61635697, 0.12985441, 0.73410642, 0.77318621,
    0.75675988, 0.03003641, 0.77527201, 0.78673995, 0.05049178,
    0.55139476, 0.02665514, 0.61664748, 0.81121749, 0.05521697,
    0.63404375, 0.32649395, 0.36828268, 0.68981099, 0.02874863,
    0.61574739])
X = np.linspace(0, 1, len(Y))

затем сделайте начальный график данных:

plt.plot(X, Y, '.')

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

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

poly_degree = 5
sd_cutoff = 1 # 2 keeps everything

coeffs = np.polyfit(X, Y, poly_degree)
poly_eqn = np.poly1d(coeffs)

Y_hat = poly_eqn(X)
delta = Y - Y_hat
sd_p = np.std(delta)

ok = abs(delta) < sd_p * sd_cutoff

, надеюсь, это имеет смысл, я использую полином более высокой степени и только отсечение в 1SD, потому что иначе ничего не будет выброшено.массив ok содержит True значения для тех точек, которые находятся в пределах sd_cutoff стандартных отклонений

, чтобы проверить это, я бы тогда сделал другой график.что-то вроде:

plt.scatter(X, Y, color=np.where(ok, 'k', 'r'))
plt.fill_between(
    X,
    Y_hat - sd_p * sd_cutoff, 
    Y_hat + sd_p * sd_cutoff,
    color='#00000020')
plt.plot(X, Y_hat)

, что дает мне:

, поэтому черные точки являются точками для сохранения (то есть X[ok] возвращает мне это, а np.where(ok) дает вам указание).

вы можете поиграться с параметрами, но вы, вероятно, захотите распределение с более толстыми хвостами (например, T-распределение Стьюдента), но, как ясказал выше, использование Google для обнаружения выбросов будет моим предложением

Answer 3

Для решения этой проблемы необходимы три функции.Во-первых, необходима функция подгонки линии для подгонки линии к набору точек:

def fit_line(x_values, vals, poly_degree):
    coeffs = np.polyfit(x_values, vals, poly_degree)
    poly_eqn = np.poly1d(coeffs)
    y_hat = poly_eqn(x_values)
    return poly_eqn, y_hat

Нам необходимо знать стандартное отклонение от точек к линии.Эта функция вычисляет это стандартное отклонение:

def compute_sd(x_values, vals, y_hat):
    distances = []
    for x,y, y1 in zip(x_values, vals, y_hat): distances.append(abs(y - y1))
    return np.std(distances)

Наконец, нам нужно сравнить расстояние от точки до линии.Точка должна быть выброшена, если расстояние от точки до линии превышает стандартное отклонение более чем в два раза.

def compare_distances(x_values, vals):    
    new_vals, new_x_vals = [],[]
    for x,y in zip(x_values, vals):    
        y1 = np.polyval(poly_eqn, x)
        distance = abs(y - y1)
        if distance < 2*sd:
            plt.plot((x,x),(y,y1), c='g')
            new_vals.append(y)
            new_x_vals.append(x)
        else:
            plt.plot((x,x),(y,y1), c='r')
            plt.scatter(x,y, c='r')
    return new_vals, new_x_vals

Как видно на следующих графиках, этот метод плохо работает дляПодгонка линии к данным, которые имеют много выбросов.Все точки в конечном итоге удаляются из-за того, что они находятся слишком далеко от установленной линии.

while len(vals)>0:
    poly_eqn, y_hat = fit_line(x_values, vals, poly_degree)
    plt.scatter(x_values, vals)
    plt.plot(x_values, y_hat)
    sd = compute_sd(x_values, vals, y_hat)
    new_vals, new_x_vals = compare_distances(x_values, vals)
    plt.show()
    vals, x_values = np.array(new_vals), np.array(new_x_vals)