RANSAC初识

RANSAC算法：随机抽样一致算法（random sample consensus,RANSAC）

一个简单的例子是从一组观测数据中找出合适的二维直线。假设观测数据中包含局内点和局外点，其中局内点近似的被直线所通过，而局外点远离于直线。简单的最小二乘法不能找到适应于局内点的直线，原因是最小二乘法尽量去适应包括局外点在内的所有点。相反，RANSAC能得出一个仅仅用局内点计算出模型，并且概率还足够高。但是，RANSAC并不能保证结果一定正确，为了保证算法有足够高的合理概率，我们必须小心的选择算法的参数。

RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点，并用下述方法进行验证：
    1.有一个模型适应于假设的局内点，即所有的未知参数都能从假设的局内点计算得出。
    2.用1中得到的模型去测试所有的其它数据，如果某个点适用于估计的模型，认为它也是局内点。
    3.如果有足够多的点被归类为假设的局内点，那么估计的模型就足够合理。
    4.然后，用所有假设的局内点去重新估计模型，因为它仅仅被初始的假设局内点估计过。
    5.最后，通过估计局内点与模型的错误率来评估模型。
    这个过程被重复执行固定的次数，每次产生的模型要么因为局内点太少而被舍弃，要么因为比现有的模型更好而被选用。

下面用一张图来模拟一下数据拟合的过程。

与最小二乘法的区别：

RANSA通过多次迭代，计算一个模型包含的局内点，模型好坏的判别依据是：模型包含越多的局内点，则这个模型越好。

而最小二乘法则是计算出一条直线，让所有的样本距离这条直线的距离最短。所以判别模型的标准就是样本距离直线的距离和。

优缺点：

RANSAC 算法的优点是能鲁棒的估计模型参数。例如，他能从包含大量局外点的数据集中估计出高精度的参数。

缺点是它计算参数的迭代次数没有上限，如果设置迭代次数的上限，得到的结果可能不是最优的结果，甚至可能得到错误的结果。

RANSAC只有一定的概率得到的可信的模型，概率与迭代次数成正比。另一个缺点是它要求设置跟问题相关的阈值，

RANSAC职能从特定的数据集中估计出一个模型，如果存在两个（或多个）模型，RANSAC不能找到别的模型。

输入：
data —— 一组观测数据
model —— 适应于数据的模型
n —— 适用于模型的最少数据个数
k —— 算法的迭代次数
t —— 用于决定数据是否适应于模型的阀值
d —— 判定模型是否适用于数据集的数据数目
输出：
best_model —— 跟数据最匹配的模型参数（如果没有找到好的模型，返回null）
best_consensus_set —— 估计出模型的数据点
best_error —— 跟数据相关的估计出的模型错误

附上一个可以实现的代码：

# -*- coding: utf-8 -*-
import numpy
import scipy  # use numpy if scipy unavailable
import scipy.linalg  # use numpy if scipy unavailable
import pylab
def ransac(data, model, n, k, t, d, debug=False, return_all=False):'''fit model parameters to data using the RANSAC algorithm
This implementation written from pseudocode found at
Given:data - a set of observed data points # 可观测数据点集model - a model that can be fitted to data points #n - the minimum number of data values required to fit the model # 拟合模型所需的最小数据点数目k - the maximum number of iterations allowed in the algorithm # 最大允许迭代次数t - a threshold value for determining when a data point fits a model #确认某一数据点是否符合模型的阈值d - the number of close data values required to assert that a model fits well to data
Return:bestfit - model parameters which best fit the data (or nil if no good model is found)'''iterations = 0bestfit = Nonebesterr = numpy.infbest_inlier_idxs = Nonewhile iterations < k:maybe_idxs, test_idxs = random_partition(n, data.shape[0])maybeinliers = data[maybe_idxs, :]test_points = data[test_idxs]maybemodel = model.fit(maybeinliers)test_err = model.get_error(test_points, maybemodel)also_idxs = test_idxs[test_err < t]  # select indices of rows with accepted pointsalsoinliers = data[also_idxs, :]if debug:print'test_err.min()', test_err.min()print'test_err.max()', test_err.max()print'numpy.mean(test_err)', numpy.mean(test_err)print'iteration %d:len(alsoinliers) = %d' % (iterations, len(alsoinliers))if len(alsoinliers) > d:betterdata = numpy.concatenate((maybeinliers, alsoinliers))bettermodel = model.fit(betterdata)better_errs = model.get_error(betterdata, bettermodel)thiserr = numpy.mean(better_errs)if thiserr < besterr:bestfit = bettermodelbesterr = thiserrbest_inlier_idxs = numpy.concatenate((maybe_idxs, also_idxs))iterations += 1if bestfit is None:raise ValueError("did not meet fit acceptance criteria")if return_all:return bestfit, {'inliers': best_inlier_idxs}else:return bestfitdef random_partition(n, n_data):#返回n行数据all_idxs = numpy.arange(n_data)numpy.random.shuffle(all_idxs)idxs1 = all_idxs[:n]idxs2 = all_idxs[n:]return idxs1, idxs2#线性系统可以使用最小二乘法求解，这个类是定义最小二乘法求解的一个示例
class LinearLeastSquaresModel:def __init__(self, input_columns, output_columns, debug=False):self.input_columns = input_columnsself.output_columns = output_columnsself.debug = debugdef fit(self, data):A = numpy.vstack([data[:, i] for i in self.input_columns]).TB = numpy.vstack([data[:, i] for i in self.output_columns]).Tx, resids, rank, s = scipy.linalg.lstsq(A, B)return xdef get_error(self, data, model):A = numpy.vstack([data[:, i] for i in self.input_columns]).TB = numpy.vstack([data[:, i] for i in self.output_columns]).TB_fit = scipy.dot(A, model)err_per_point = numpy.sum((B - B_fit) ** 2, axis=1)  # sum squared error per rowreturn err_per_pointdef test():#创建一些数据（局内点）使其在一条直线附近n_samples = 500n_inputs = 1n_outputs = 1A_exact = 20 * numpy.random.random((n_samples, n_inputs))  # x坐标perfect_fit = 60 * numpy.random.normal(size=(n_inputs, n_outputs))  # the model(斜率)B_exact = scipy.dot(A_exact, perfect_fit)  # y坐标#验证y坐标数组的大小,在数组计算中尺寸大小很重要，如果尺寸大小不匹配，则导致计算错误，或者无法计算assert B_exact.shape == (n_samples, n_outputs)#pylab.plot( A_exact, B_exact, 'b.', label='data' ) #这两行是图片显示创建出来的数据，不含噪声的数据#pylab.show()#添加一点点高斯噪声（仅用最小二乘的方法即可很好的去解决小噪声的线性回归问题）A_noisy = A_exact + numpy.random.normal(size=A_exact.shape)  # x坐标添加高斯噪声B_noisy = B_exact + numpy.random.normal(size=B_exact.shape)  # y坐标....#pylab.plot( A_noisy, B_noisy, 'b.', label='data')  #这两句是图片显示带有噪声的数据。#pylab.show()#添加局外点，对模型进行干扰n_outliers = 100  # 500个数据点有100个是putliersall_idxs = numpy.arange(A_noisy.shape[0])numpy.random.shuffle(all_idxs)  # 索引随机排列outlier_idxs = all_idxs[:n_outliers]  # 选取all_idxs前100个做outlier_idxsnon_outlier_idxs = all_idxs[n_outliers:]  # 后面的不是outlier_idxsA_noisy[outlier_idxs] = 20 * numpy.random.random((n_outliers, n_inputs))  # 外点的横坐标B_noisy[outlier_idxs] = 50 * numpy.random.normal(size=(n_outliers, n_outputs))  # 外点的纵坐标pylab.plot( A_noisy, B_noisy, 'b.', label='data' )pylab.show()#建立模型all_data = numpy.hstack((A_noisy, B_noisy))  # 组成坐标对input_columns = range(n_inputs)  # the first columns of the arrayoutput_columns = [n_inputs + i for i in range(n_outputs)]  # the last columns of the arraydebug = Falsemodel = LinearLeastSquaresModel(input_columns, output_columns, debug=debug)linear_fit, resids, rank, s = scipy.linalg.lstsq(all_data[:, input_columns],all_data[:, output_columns])#使用RANSAC算法ransac_fit, ransac_data = ransac(all_data, model,50, 1000, 7000, 300,  # misc. parametersdebug=debug, return_all=True)sort_idxs = numpy.argsort(A_exact[:, 0])  # 对A_exact排序， sort_idxs为排序索引A_col0_sorted = A_exact[sort_idxs]  # maintain as rank-2 arraypylab.plot(A_noisy[:, 0], B_noisy[:, 0], 'k.', label='data')pylab.plot(A_noisy[ransac_data['inliers'], 0], B_noisy[ransac_data['inliers'], 0], 'bx',label='RANSAC data')pylab.plot(A_col0_sorted[:, 0],numpy.dot(A_col0_sorted, ransac_fit)[:, 0] ,label='RANSAC fit')pylab.plot(A_col0_sorted[:, 0],numpy.dot(A_col0_sorted, perfect_fit)[:, 0],label='exact system')pylab.plot(A_col0_sorted[:, 0],numpy.dot(A_col0_sorted, linear_fit)[:, 0] ,label='linear fit')pylab.legend()pylab.show()
#运行主函数
if __name__ == '__main__':test()

这是我自己创建的数据，数据点多的那一大块是我创建的局内点，下面那一小排是我创建的局外点。