python机器学习手写算法系列

2023-10-14 20:42| 来源: 网络整理| 查看: 265

本文旨在通过在最简单的一元一次方程中运用RANSAC，编写代码，来学习RANSAC算法。

RANSAC算法

RANSAC[1]算法，全程是Random Sample Consensus（随机抽样一致）。它采用迭代的方式从一组包含outliers的数据中估算出回归模型。RANSAC是一个非确定性算法，因为它是“随机抽样”的。

它原来是计算机视觉的算法，后来被用来做线性回归。

RANSAC回归用最少的数据点训练模型，基本算法总结如下：

选取的最少的数据点，作为样本。用样本和基础回归器训练模型。根据 ϵ \epsilon ϵ，选出inliers和outliers如果inliers的数量达到了阈值 τ \tau τ，用inliers和基础回归器训练模型判断终止条件。达到，退出。否则，返回1。（最多 N N N次）

迭代次数 N N N，表示在 N N N次迭代以后，至少有一次迭代，任意选取的 m m m个点，都是inliers的概率是 p p p（通常为0.99）。如果inliers的比例是 u u u，我们就有以下公式。 1 − p = ( 1 − u m ) N 1-p=(1-u^m)^N 1−p=(1−um)N 变换后得到： N = l o g ( 1 − p ) l o g ( 1 − u m ) N = \frac{log(1-p)}{log(1-u^m)} N=log(1−um)log(1−p)

其实，这里我有个疑问。就算，我选择的两个点都是inliers，也无法保证他们他们组成的直线能够穿过其他的inliers。这个就说到RANSAC的原始用途了

在这里插入图片描述

在原始的RANSAC算法里，所谓的样本，是指两张照片上面的相同点。如果找到了这样的样本，当然就可以确定他们之间的向量了。但是回归不是这么回事。如果我的inliers比例是95%，那么 N = l o g ( 1 − 0.99 ) l o g ( 1 − 0.9 5 2 ) = l o g ( 0.01 ) l o g ( 1 − 0.9025 ) = − 2 l o g ( 0.0975 ) = − 2 − 1.01 = 2 N = \frac{log(1-0.99)}{log(1-0.95^2)} = \frac{log(0.01)}{log(1-0.9025)} = \frac{-2}{log(0.0975)}=\frac{-2}{-1.01}=2 N=log(1−0.952)log(1−0.99)=log(1−0.9025)log(0.01)=log(0.0975)−2=−1.01−2=2 3次就能找到模型了？这个答案，我们要去看Scikit-Learn的源代码。

Scikit-Learn中的RANSAC算法

Talk is cheap, show me the code. 为了弄明白RANSAC的具体实现，我找到了sklearn的RANSAC[2]。通过其源代码，可以了解其具体实现。

第一步，选取数据点。sklearn的默认值是X.shape[1] + 1。我们的手写代码是一元一次方程，所以，对于我们来说，这个值就是2。

第二步，用样本训练模型。这里，我们的样本数量是2，可以直接用两点决定一条直线的方法，确定这个模型。

第三步，根据 ϵ \epsilon ϵ，选出inliers和outliers。这里scikit-learn把 ϵ \epsilon ϵ叫做residual_threshold，其默认值是MAD (median absolute deviation，偏差的绝对值的中位数)

第四步，如果inliers的数量达到了阈值 τ \tau τ，用inliers和基础回归器训练模型。这里，Scikit-Learn是默认是没有阈值的。（你可以传一个函数进去，这个函数的默认值是None）。我们这里就设置成50%。也就是说，至少50%的点都是inliers，才用他们训练模型。

第五步，终止条件。这里，Scikit-Learn用的是动态值，也就是说，inliers的比例 u u u，是随着找到的inliers的数量变化的。这时，即使我们找到的两个点，都是inliers，但是他们组成的直线却和真实值大相径庭，也不用担心。因为这时的 u u u很小，所以 N N N，比较大。随着我们找到越来越多的inliers， u u u变大，相应的 N N N变小。

手写RANSAC算法

通过几篇文章和scikit-learn的源代码，大致对RANSAC有所了解了。接下来就是通过自己写算法来验证自己的理解。

Scikit Learn示例

首先，我们复用scikit-learn的示例。

https://scikit-learn.org/stable/auto_examples/linear_model/plot_ransac.html

import numpy as np from matplotlib import pyplot as plt from sklearn import linear_model, datasets n_samples = 1000 n_outliers = 50 X, y, coef = datasets.make_regression(n_samples=n_samples, n_features=1, n_informative=1, noise=10, coef=True, random_state=0) # Add outlier data np.random.seed(0) X[:n_outliers] = 3 + 0.5 * np.random.normal(size=(n_outliers, 1)) y[:n_outliers] = -3 + 10 * np.random.normal(size=n_outliers) # Fit line using all data lr = linear_model.LinearRegression() lr.fit(X, y) # Robustly fit linear model with RANSAC algorithm ransac = linear_model.RANSACRegressor() ransac.fit(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(line_X) line_y_ransac = ransac.predict(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()

在这里插入图片描述

接着，我们就可以手写算法了。

探索

我们先计算阈值 ϵ \epsilon ϵ，也就是scikit-learn里面的residual_threshold

from scipy import stats import numba residual_threshold = stats.median_abs_deviation(y) print(f"Residual threshold is {residual_threshold}.")

接着，我们取两个outliers，两个inliers，以及分别取一个。

探索代码：

n_samples = len(y) @numba.njit() def one_run(i, j): # build a model with Point I and Point J slope = (y[i] - y[j]) / (X[i][0] - X[j][0]) # y = ax + b # b = y - ax intercept = y[i] - slope * X[i][0] # Count number of inliers inliers = [] for k in range(n_samples): y_hat = slope * X[k][0] + intercept if y[k] > y_hat - residual_threshold and y[k]

【本文地址】

公司简介

联系我们