python 에서 RANSAC 으로 polynomial fitting 방법

RANSAC은 scikit-learn 에 구현되어있고, line fitting 하는 example code 도 Robust linear model estimation using RANSAC에 친절하게 나와있다.

이걸 multiple order polynomial regression 으로 확장하기 위해서는 feature를 다항식에 맞게 확장만 해주면 된다.

예를 들어 $ax + b = y$ 에서 feature set 인 [X]를 $ax^2 + bx + c = y$ 에서는 feature set을 [X^2, X] 으로 확장시켜 주고 feature 가 2개인 linear regression 으로 각각의 coefficient 를 찾아주면 된다.

Robust linear model estimation using RANSAC 에서 fitting curve 와 dataset을 2nd order polynomial 로 바꾼 코드는 다음과 같다. 이 코드를 실행해보면 위 그림과 같이 2nd order polynomial fitting 을 RANSAC 으로 수행한 결과를 테스트할 수 있다.

outlier sample 이 많음에도 polynomial fitting 을 정확하게 수행하고 있다. !!

import numpy as np
from matplotlib import pyplot as plt

from sklearn import linear_model
np.random.seed(0)


n_samples = 1000
n_outliers = 50


def get_polynomial_samples(n_samples=1000):
    X = np.array(range(1000)) / 100.0
    np.random.shuffle(X)

    coeff = np.random.rand(2,) * 3

    y = coeff[0]*X**2 + coeff[1]*X + 10
    X = X.reshape(-1, 1)
    return coeff, X, y


def add_square_feature(X):
    X = np.concatenate([(X**2).reshape(-1,1), X], axis=1)
    return X


coef, X, y = get_polynomial_samples(n_samples)

# Add outlier data
X[:n_outliers] = 10 + 0.5 * np.random.normal(size=(n_outliers, 1))
y[:n_outliers] = -10 + 10 * np.random.normal(size=n_outliers)

# Fit line using all data
lr = linear_model.LinearRegression()
lr.fit(add_square_feature(X), y)

# Robustly fit linear model with RANSAC algorithm
ransac = linear_model.RANSACRegressor()
ransac.fit(add_square_feature(X), y)
inlier_mask = ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)

# Predict data of estimated models
line_X = np.arange(X.min(), X.max())[:, np.newaxis]
line_y = lr.predict(add_square_feature(line_X))
line_y_ransac = ransac.predict(add_square_feature(line_X))

# Compare estimated coefficients
print("Estimated coefficients (true, linear regression, RANSAC):")
print(coef, lr.coef_, ransac.estimator_.coef_)

lw = 2
plt.scatter(X[inlier_mask], y[inlier_mask], color='yellowgreen', marker='.',
            label='Inliers')
plt.scatter(X[outlier_mask], y[outlier_mask], color='gold', marker='.',
            label='Outliers')
plt.plot(line_X, line_y, color='navy', linewidth=lw, label='Linear regressor')
plt.plot(line_X, line_y_ransac, color='cornflowerblue', linewidth=lw,
         label='RANSAC regressor')
plt.legend(loc='lower right')
plt.xlabel("Input")
plt.ylabel("Response")
plt.show()