线性回归

/*************************************************
Author:    Sayheyheyhey

Date:2019-7-4

Description:根据伪代码实现通用的RANSAC模板
        自定义线性模型，实现两种方式的直线拟合
**************************************************/

#include <random>
#include <iostream>
#include <time.h>
#include <set>
#include <cassert>
#include <limits.h>

using namespace std;
//数据点类型
struct st_point{
    st_point(){};
    st_point(double X, double Y) :x(X), y(Y){};
    double x;
    double y;
};
/**
  * @brief 线性模型
  *
  * Ax+By+C = 0;
*/
class linearModel{
public:
    //待估计参数
    double A, B, C;
public:
    linearModel(){};
    ~linearModel(){};
    
    //使用两个点对直线进行初始估计
    void Update(vector<st_point> &data, set<int> &maybe_inliers){
        assert(maybe_inliers.size() == 2);        //初始化的点不为2个，报错
        //根据索引读取数据
        vector<int> points(maybe_inliers.begin(), maybe_inliers.end());
        st_point pts1 = data[points[0]];
        st_point pts2 = data[points[1]];
        //根据两个点计算直线参数（得到其中一组解，可以任意比例缩放）
        double delta_x = pts2.x - pts1.x;
        double delta_y = pts2.y - pts1.y;
        A = delta_y;
        B = -delta_x;
        C = -delta_y*pts2.x + delta_x*pts2.y;
    }

    //返回点到直线的距离
    double computeError(st_point point){
        double numerator = abs(A*point.x + B*point.y + C);
        double denominator = sqrt(A*A + B*B);
        return numerator / denominator;
    }

    //根据一致点的集合对直线进行重新估计
    double Estimate(vector<st_point> &data, set<int> &consensus_set){
        assert(consensus_set.size() >= 2);
        //求均值 means
        double mX, mY;
        mX = mY = 0;
        for (auto &index : consensus_set){
            mX += data[index].x;
            mY += data[index].y;
        }
        mX /= consensus_set.size();
        mY /= consensus_set.size();
        
        //求二次项的和 sum
        double sXX, sYY, sXY;
        sXX = sYY = sXY = 0;
        for (auto &index : consensus_set){
            st_point point;
            point = data[index];
            sXX += (point.x - mX)*(point.x - mX);
            sYY += (point.y - mY)*(point.y - mY);
            sXY += (point.x - mX)*(point.y - mY);
        }
        
        //解法1：求y=kx+b的最小二乘估计，然后再转换成一般形式
        //参考 https://blog.csdn.net/hookie1990/article/details/91406309
        bool isVertical = sXY == 0 && sXX < sYY;
        bool isHorizontal = sXY == 0 && sXX > sYY;
        bool isIndeterminate = sXY == 0 && sXX == sYY;
        double k = NAN;
        double b = NAN;

        if (isVertical)
        {
            A = 1;
            B = 0;
            C = mX;
        }
        else if (isHorizontal)
        {
            A = 0;
            B = 1;
            C = mY;
        }
        else if (isIndeterminate)
        {
            A = NAN;
            B = NAN;
            C = NAN;
        }
        else
        {
            k = (sYY - sXX + sqrt((sYY - sXX) * (sYY - sXX) + 4.0 * sXY * sXY)) / (2.0 * sXY);    //斜率
            b = mY - k * mX;                                                            //截距
            //正则化项，使得A^2+B^2 = 1;
            double normFactor = 1 / sqrt(1 + k*k);
            A = normFactor * k;
            B = -normFactor;
            C = normFactor*b;
        }
        //返回残差
        if (isIndeterminate){
            return NAN;
        }
        double error = A*A*sXX + 2 * A*B*sXY + B*B*sYY;
                error /= consensus_set.size();
                return error;
        /*
        //解法2：
                if(sXX == 0){
                    A = 1;    
                    B = 0;
                    C = -mX;
                }
                else{
                    A = sXY/sXX;
                    B = -1;
                    C = mY - A*mX;
                    //归一化令A^2+B^2 = 1;
                    double normFactor = sqrt(A*A+B*B);
                    A /= normFactor;
                    B /= normFactor;
                    C /= normFactor;
                }
        double error = A*A*sXX + 2 * A*B*sXY + B*B*sYY;
                error /= consensus_set.size();    //求平均误差
        return error;
        */
    }
};


/**
* @brief 运行RANSAC算法
*
* @param[in]    data    一组观测数据
* @param[in]    n        适用于模型的最少数据个数
* @param[in]    k        算法的迭代次数
* @param[in]    t        用于决定数据是否适应于模型的阀值
* @param[in]    d        判定模型是否适用于数据集的数据数目 
* @param[in&out]    model    自定义的待估计模型，为该函数提供Update、computeError和Estimate三个成员函数
*                            运行结束后，模型参数被设置为最佳的估计值
* @param[out]    best_consensus_set    输出一致点的索引值
* @param[out]    best_error    输出最小损失函数
*/
template<typename T, typename U>
int ransac(vector<T> &data, int n, int k, double t, int d,
            U &best_model,set<int> &best_consensus_set, double &best_error){
    //1.初始化
    int  iterations = 0;    //迭代次数
    U maybe_model;            //使用随机选点初始化求得的模型
    U better_model;            //根据符合条件的一致点拟合出的模型
    
    int isFound = 0;                    //算法成功的标志
    set<int> maybe_inliers;                //初始随机选取的点（的索引值）

    //best_error = DBL_MAX;    //初始化为最大值
    best_error = 1.7976931348623158e+308;
    default_random_engine rng(time(NULL));                    //随机数生成器
    uniform_int_distribution<int> dist(0, data.size()-1);    //采用均匀分布
    
    //2.主循环
    while (iterations < k){
        //3.随机选点
        maybe_inliers.clear();    
        while (1){
            int index = dist(rng);
            maybe_inliers.insert(index);
            if (maybe_inliers.size() == n){
                break;
            }
        }
        //4.计算初始值
        maybe_model.Update(data, maybe_inliers);                                //自定义函数，更新模型
        set<int> consensus_set(maybe_inliers.begin(),maybe_inliers.end());        //选取模型后，根据误差阈值t选取的内点(的索引值)

        //5.根据初始模型和阈值t选择内点    
        for (int i = 0; i < data.size(); i++){
            double error_per_item = maybe_model.computeError(data[i]);
            if (error_per_item < t){
                consensus_set.insert(i);
            }
        }
        //6.根据全部的内点重新计算模型
        if (consensus_set.size() > d){
            double this_error = better_model.Estimate(data, consensus_set);        //自定义函数，（最小二乘）更新模型，返回计算出的误差
            //7.若当前模型更好，则更新输出量
            if (this_error < best_error){
                best_model = better_model;
                best_consensus_set = consensus_set;
                best_error = this_error;
            }
            isFound = 1;
        }
        ++iterations;
    }
    return isFound;
}


int main(){
    //1.读入数据
    int data_size;        //输入第一行表示数据大小
    cin >> data_size;    
    vector<st_point> Points(data_size);
    for (int i = 0; i < data_size; i++){
        cin >> Points[i].x >> Points[i].y;
    }
    //测试用
    //vector<st_point> Points{ st_point(3, 4), st_point(6, 8), st_point(9, 12), st_point(15, 20), st_point(10,-10)};
    //int data_size = Points.size();
    //2.设置输入量
    int k = 50;                //最大迭代次数
    int n = 2;                //适用于模型的最少数据个数
    double t = 0.01;        //用于决定数据是否适应于模型的阀值
    int d = data_size*0.5;    //判定模型是否适用于数据集的数据数目 
    //3.初始化输出量
    linearModel best_model;            //最佳线性模型
    set<int> best_consensus_set;    //记录一致点索引的set
    double best_error;                //最小残差
    //4.运行RANSAC            
    int status = ransac(Points, n, k, t, d, best_model, best_consensus_set, best_error);
    //5.输出
    cout << best_model.A << " " << best_model.B << " " << best_model.C << endl;
    return 0;
}
解法1可以全部测试通过;
解法2只能通过77%;
希望哪位大佬能告诉我这是为什么...想不出这两种公式有什么区别

一般直线拟合是2种，ransac和最小二乘法，该题用ransac