直线拟合、二次曲线拟合、折线拟合和KNN近邻（附代码）

冷不防 2022-08-27 06:49 323阅读 0赞

一个工程中的应用，需要对一组数据做上面四种形式的拟合回归，并且根据模型对输入做evaluation，就是做一个函数曲线拟合。

下面的RevPre定义了方法和结构，Util是使用案例，其中的Opt的type指示模型要拟合的是哪一种。

1）直线拟合： y = kx+b

2) 二次曲线拟合：y = AX^2 + BX + C

以上两种很典型，不多解释；

3） KNN

KNN的Opt阶段什么也没做，只是保存了所有数据对，然后检测阶段选取最近的K（MIN_KNN_K < K < MAX_KNN_K）个点来计算距离加权后的结果。

4） 折线拟合：

前面先介绍KNN，是因为折线拟合的结果与之很类似，同样保存了一系列的数据对，检测的时候判断处于哪一段。不同之处是做了简化，而不是全部存储起来。

首先对数据对按照x快排，然后以直线拟合判断总体是升序还是降序，接下来就是不断淘汰不符合顺序的点，发生“突起”或者“凹陷”时清除与拟合直线距离远的点，直到所有点都按序排列，形成一条单调折线。

代码：

RevPre.h

  1. #include <iostream>
  2. #include <math.h>
  4. /************************Model type*********************************/
  5. #define MAX_POL_DEPTH    3
  6. #define MAX_KNN_K            10
  7. #define MIN_KNN_K                1
  9. #define MAX(x,y)        (x) > (y) ? (x) : (y)
  10. #define MIN(x,y)            (x) < (y) ? (x) : (y)
  11. enum modelType{
  12.     StraightLine = 0, // default
  13.     CurveAt2,
  14.     BrokenLine,
  15.     KNNModel
  16. };
  17. typedef struct Model{
  18.     enum modelType type;
  19.     // Line parameters
  20.     double            lineParam[MAX_POL_DEPTH];
  21.     // Point model
  22.     double            *px, *py;
  23.     int                len;
  24. };
  25. Model* CreateModel();
  26. void ReleaseModel( Model** _ptr );
  27. bool SetOptData( Model* ptr, double *x, double *y, int len );
  28. bool Opt( Model *ptr, modelType type );
  29. double Predict( Model *ptr, double x );
  30. /**************************Polynomial*******************************/
  31. /* Internal */
  32. void CalculatePower(double *powers, int ptNum, int maxDepth, double *x );          //将初始x[i]的值的各幂次方存储在一个二维数组里面 
  33. void CalculateParams(double *powers, int ptNum, int maxDepth, 
  34.                      double *params, double *y);                                                                    //计算正规方程组的系数矩阵 
  35. void DirectLU( double *params, int ptNum, int maxDepth, double *x );                        //列主元LU分解
  36. inline void swap(double &,double &);                                                                        //交换两个变量的值
  37. /* External */
  38. bool PolynomialOpt( Model *ptr );
  39. /************************StraightLine********************************/
  40. bool StraightLineOpt( Model *ptr );
  41. /************************BrokenLine********************************/
  42. /*Internal*/
  43. int SingleSort( double *index, double *context, int start, int end );
  44. void QuickSort( double *index, double *context, int start, int end );
  45. int CheckSequence( double *context, int start, int end, bool upTrend );
  46. /*External*/
  47. bool BrokenLineOpt( Model *ptr );
  48. /********************KNN(Lazy-learning)****************************/
  49. bool KNNOpt( Model *ptr );

RevPre.cpp

  1. #include "Revise.h"
  3. Model* CreateModel(){
  4.     Model *ptr = new Model;
  5.     ptr->type = StraightLine;
  7.     ptr->px = ptr->py = NULL;
  8.     ptr->len = 0;
  9.     memset( ptr->lineParam, 0, sizeof(double)*MAX_POL_DEPTH );
  11.     return ptr;
  12. }
  14. void ReleaseModel( Model** _ptr ){
  15.     Model *ptr = *_ptr;
  16.     if ( ptr->px ) delete[] ptr->px;
  17.     if ( ptr->py ) delete[] ptr->py;
  18.     delete ptr;
  20.     *_ptr = NULL;
  21.     return ;
  22. }
  24. bool SetOptData( Model *ptr, double *x, double *y, int len ){
  25.     if ( !ptr || !x || !y ) return false;
  27.     if ( !ptr->px ) ptr->px = new double[len];
  28.     if ( !ptr->py ) ptr->py = new double[len];
  30.     ptr->len = len;
  31.     memcpy( ptr->px, x, sizeof(double)*len );
  32.     memcpy( ptr->py, y, sizeof(double)*len );
  33.     return true;
  34. }
  36. bool Opt( Model *ptr, modelType type ){
  37.     if ( !ptr ) return false;
  38.     switch( type )
  39.     {
  40.     case StraightLine:
  41.         return    StraightLineOpt( ptr );
  42.     case CurveAt2:
  43.         return    PolynomialOpt( ptr );
  44.     case BrokenLine:
  45.         return    BrokenLineOpt( ptr );
  46.     case KNNModel:
  47.         return    KNNOpt( ptr );
  48.     default:
  49.         return    false;
  50.     }
  51. }
  53. double Predict( Model *ptr, double x ){
  54.     if ( !ptr ) exit (-1);
  55.     switch( ptr->type )
  56.     {
  57.     case StraightLine:
  58.         return ptr->lineParam[0] + ptr->lineParam[1]*x;
  59.     case CurveAt2:
  60.         return ptr->lineParam[0] + ptr->lineParam[1]*x + ptr->lineParam[2]*x*x;
  61.     case BrokenLine:
  62.         {
  63.             if ( ptr->len < 3 ) exit(-2);
  64.             int first = 0;
  65.             if ( x <= ptr->px[0] )
  66.             {
  67.                 double x0 = ptr->px[0], x1 = ptr->px[1];
  68.                 double y0 = ptr->py[0], y1 = ptr->py[1];
  69.                 return y0 - (x0-x)*(y1-y0)/(x1-x0);
  70.             }
  71.             else if ( x >= ptr->px[ptr->len-1] )
  72.             {
  73.                 double x0 = ptr->px[ptr->len-2], y0 = ptr->py[ptr->len - 2];
  74.                 double x1 = ptr->px[ptr->len-1], y1 = ptr->py[ptr->len - 1];
  75.                 return y1 -(x-x1)*(y0-y1)/(x1-x0);
  76.             }
  77.             else
  78.             {
  79.                 while ( ptr->px[first] < x ) { first ++ ;}
  80.                 first --;
  81.                 double deltay = ptr->py[first+1] - ptr->py[first];
  82.                 double deltax = ptr->px[first+1] - ptr->px[first];
  83.                 return ptr->py[first] + deltay*(x-ptr->px[first])/deltax;
  84.             }
  85.         }
  86.     case KNNModel:
  87.         {
  88.             int K = MAX( MIN_KNN_K, MIN( int(ptr->len*0.1), MAX_KNN_K ) );
  89.             // Prepare the initial K neighbours
  90.             double *dist_team = new double[K];
  91.             int        *idx_team = new int[K];
  92.             int        farestIdt = -1;
  93.             double    farestDist = 0;
  95.             int id = 0;
  96.             for ( ; id < K; id ++ )
  97.             {
  98.                 idx_team[id] = id;
  99.                 dist_team[id] = abs( ptr->px[id] - x );
  100.                 if ( farestDist <= dist_team[id] )
  101.                 {
  102.                     farestIdt = id;
  103.                     farestDist = dist_team[id];
  104.                 }
  105.             }
  106.             // Looking for the K nearest neighbours
  107.             while ( id < ptr->len )
  108.             {
  109.                 if ( abs( ptr->px[id] -x ) < farestDist )
  110.                 {
  111.                     // Update the team
  112.                     idx_team[farestIdt] = id;
  113.                     dist_team[farestIdt] = abs( ptr->px[id] - x );
  114.                     // Update the farest record
  115.                     farestIdt = 0;
  116.                     farestDist = dist_team[0];
  117.                     for ( int searchIdt = 1; searchIdt < K; searchIdt ++ )
  118.                     {
  119.                         if ( dist_team[searchIdt] > farestDist )
  120.                         {
  121.                             farestDist = dist_team[searchIdt];
  122.                             farestIdt = searchIdt;
  123.                         }
  124.                     }
  125.                 }
  126.                 id ++;
  127.             }
  129.             // Calculate their contribution
  130.             double res = 0.0;
  131.             double weightSum = 0.0;
  132.             for ( int seachIdt = 0; seachIdt < K; seachIdt ++ )
  133.             {
  134.                 weightSum += 1.0/dist_team[seachIdt];
  135.                 res += 1.0/dist_team[seachIdt]*ptr->py[idx_team[seachIdt]];
  136.             }
  137.             delete[] dist_team;
  138.             delete[] idx_team;
  139.             return res/weightSum;
  140.         }
  141.     default:
  142.         exit(-2);
  143.     }
  144. }
  145. /**************************Polynomial*******************************/
  146. bool StraightLineOpt( Model *ptr )
  147. {
  148.     if ( !ptr ) return false; 
  149.     if ( !ptr->px || !ptr->py ) return false;
  151.     int outLen = 2;
  152.     int ptNum = ptr->len, maxDepth = outLen;
  154.     double *powers = new double[maxDepth*ptNum];
  155.     double *params = new double[maxDepth*(maxDepth+1)];
  157.     CalculatePower( powers, ptNum, maxDepth, ptr->px );
  159.     CalculateParams( powers, ptNum, maxDepth, params, ptr->py ); //计算正规方程组的系数矩阵
  161.     DirectLU( params, ptNum, maxDepth, ptr->lineParam ); //列主元LU分解
  162.     ptr->type = StraightLine;
  164.     std::cout<<"-------------------------"<<std::endl;
  165.     std::cout<<"拟合函数的系数分别为:\n";
  166.     for( int i=0;i<maxDepth;i++)
  167.         std::cout<<"a["<<i<<"]="<<ptr->lineParam[i]<<std::endl;
  168.     std::cout<<"-------------------------"<<std::endl;
  170.     delete[] powers;
  171.     delete[] params;
  173.     return true;
  174. }
  175. bool PolynomialOpt( Model *ptr )
  176. {
  177.     if ( !ptr ) return false; 
  178.     if ( !ptr->px || !ptr->py ) return false;
  180.     int outLen = MAX_POL_DEPTH;
  181.     int ptNum = ptr->len, maxDepth = outLen;
  183.     double *powers = new double[maxDepth*ptNum];
  184.     double *params = new double[maxDepth*(maxDepth+1)];
  186.     CalculatePower( powers, ptNum, maxDepth, ptr->px );
  188.     CalculateParams( powers, ptNum, maxDepth, params, ptr->py ); //计算正规方程组的系数矩阵
  190.     DirectLU( params, ptNum, maxDepth, ptr->lineParam ); //列主元LU分解
  191.     ptr->type = CurveAt2;
  193.     /*std::cout<<"-------------------------"<<std::endl;
  194.     std::cout<<"拟合函数的系数分别为:\n";
  195.     for( int i=0;i<maxDepth;i++)
  196.         std::cout<<"a["<<i<<"]="<<ptr->lineParam[i]<<std::endl;
  197.     std::cout<<"-------------------------"<<std::endl;*/
  199.     delete[] powers;
  200.     delete[] params;
  202.     return true;
  203. }
  205. void CalculatePower(double *powers, int ptNum, int maxDepth, double *x )
  206. {
  207.     if ( !powers || !x ) return ;
  209.     int            i, j, k;
  210.     double        temp;
  212.     for( i = 0; i < maxDepth; i ++ )
  213.         for( j = 0; j < ptNum; j ++ )
  214.         {
  215.             temp = 1;
  216.             for( k = 0; k < i; k ++ )
  217.                 temp *= x[j];
  218.             powers[i*ptNum+j] = temp;
  219.         }
  220.     return ;
  221. }
  223. void CalculateParams(double *powers, int ptNum, int maxDepth, 
  224.                      double *params, double *y)
  225. {
  226.     if ( !powers || !params || !y ) return ;
  228.     int            i, j, k;
  229.     double        temp;
  230.     int            step = maxDepth + 1;
  232.     for( i = 0; i < maxDepth; i ++ )
  233.     {
  234.         for(j = 0; j < maxDepth; j ++ )
  235.         {
  236.             temp = 0;
  237.             for( k = 0; k < ptNum; k ++ )
  238.                 temp += powers[i*ptNum+k]*powers[j*ptNum+k];
  239.             params[i*step+j] = temp;
  240.         }
  242.         temp = 0;
  243.         for( k = 0; k < ptNum; k ++ )
  244.         {
  245.             temp += y[k]*powers[i*ptNum+k];
  246.             params[i*step+maxDepth] = temp;
  247.         }
  248.     }
  250.     return ;
  251. }
  253. inline void swap(double &a,double &b)
  254. {
  255.     a=a+b;
  256.     b=a-b;
  257.     a=a-b;
  258. }
  260. void DirectLU( double *params, int ptNum, int maxDepth, double *x )
  261. {
  262.     int                i, r, k, j;
  263.     double            max;
  264.     int                step = maxDepth + 1;
  266.     double *s = new double[maxDepth];
  267.     double *t = new double[maxDepth];
  268.     // choose the main element
  269.     for( r = 0; r < maxDepth; r ++ )
  270.     {
  271.         max = 0;
  272.         j = r;
  273.         for( i = r; i < maxDepth; i ++ ) 
  274.         {
  275.             s[i] = params[i*step+r];
  276.             for( k = 0; k < r; k ++ )
  277.                 s[i] -= params[i*step+k] * params[k*step+r];
  278.             s[i] = abs(s[i]);
  280.             if( s[i] > max ){
  281.                 j = i;
  282.                 max = s[i];
  283.             }
  284.         }
  285.         // if the "main"element is not @ row r, swap the corresponding element 
  286.         if( j != r ) 
  287.         {
  288.             for( i = 0; i < maxDepth + 1; i ++ )
  289.                 swap( params[r*step+i], params[j*step+i] );
  290.         }
  291.         for( i = r; i < step; i ++ ) 
  292.             for( k = 0; k < r; k ++ ){
  293.                 params[r*step+i] -= params[r*step+k] * params[k*step+i];
  294.             }
  295.             for(i = r+1; i < maxDepth; i ++ ) 
  296.             {
  297.                 for ( k = 0; k < r; k ++ )
  298.                     params[i*step+r] -= params[i*step+k] * params[k*step+r];
  299.                 params[i*step+r] /= params[r*step+r];
  300.             }
  301.     }
  302.     for( i = 0; i < maxDepth; i ++ )
  303.         t[i] = params[ i*step + maxDepth ];
  304.     for ( i = maxDepth - 1; i >= 0; i -- ) //利用回代法求最终解
  305.     {
  306.         for ( r = maxDepth - 1; r > i; r -- )
  307.             t[i] -= params[ i*step + r ] * x[r];
  308.         x[i] = t[i]/params[i*step+i];
  309.     }
  311.     delete[] s;
  312.     delete[] t;
  314.     return ;
  315. }
  317. /**********************Broken Line***************************/
  318. // Quick Sort
  319. int SingleSort( double *index, double *context, int start, int end )
  320. {
  321.     if ( end - start < 1 ) return start;
  323.     int i = start, j = end;
  324.     double key = index[i];
  325.     double key_ = context[i];
  326.     while ( i < j )
  327.     {
  328.         while ( index[j] > key && j > i ) j --;
  329.         if ( index[j] < key )
  330.         {
  331.             index[i] = index[j];
  332.             context[i] = context[j];
  333.         }
  334.         while ( index[i] < key && j > i ) i ++;
  335.         if ( index[i] > key )
  336.         {
  337.             index[j] = index[i];
  338.             context[j] = context[i];
  339.         }
  340.     }
  342.     index[i] = key;
  343.     context[j] = key_;
  345.     return i;
  346. }
  348. void QuickSort( double *index, double *context, int start, int end )
  349. {
  350.     if ( end - start < 1 ) return ; // important
  351.     int mid = SingleSort( index, context, start, end );
  352.     QuickSort( index, context, start, mid - 1 );
  353.     QuickSort( index, context, mid+ 1, end );
  354. }
  355. int CheckSequence( double *context, int start, int end, bool upTrend )
  356. {
  357.     int i = start;
  358.     for ( ; i < end; i ++ )
  359.     {
  360.         if ( upTrend && context[i+1] < context[i] )
  361.         {
  362.             return i;
  363.         }
  364.         if ( !upTrend && context[i] < context[i+1] )
  365.         {
  366.             return i;
  367.         }
  368.     }
  369.     return -1;
  370. }
  371. // Form the broken line
  372. bool BrokenLineOpt( Model *ptr )
  373. {
  374.     if ( !ptr ) return false;
  375.     if ( !ptr->len || !ptr->px || !ptr->py ) return false;
  377.     // analyse the trend of points and get its approximate line
  378.     StraightLineOpt( ptr );
  379.     double k = ptr->lineParam[1], b = ptr->lineParam[0];
  380.     bool upTrend = ( k > 0 );
  381.     // sort the sequence by py
  382.     QuickSort( ptr->px, ptr->py, 0, ptr->len - 1 );
  384.     int oddPoint = 0;
  385.     while ( (oddPoint = CheckSequence( ptr->py, oddPoint, ptr->len -1, upTrend ) ) != -1 )
  386.     {
  387.         double formerErr = abs( k*ptr->px[oddPoint] + b - ptr->py[oddPoint] );
  388.         double laterErr = abs( k*ptr->px[oddPoint+1] + b - ptr->py[oddPoint+1] );
  389.         oddPoint = formerErr > laterErr ? oddPoint : oddPoint + 1;
  390.         // remove the odd point
  391.         memcpy( ptr->py + oddPoint, ptr->py + oddPoint + 1, sizeof(double) );
  392.         memcpy( ptr->px + oddPoint, ptr->px + oddPoint + 1, sizeof(double) );
  393.         ptr->len --;
  394.         oddPoint --;
  395.     }
  396.     memset( ptr->lineParam, 0, sizeof(double)*MAX_POL_DEPTH );
  397.     ptr->type = BrokenLine;
  398.     return true;
  399. }
  401. /**********************Lazy Learning***************************/
  402. bool KNNOpt( Model *ptr )
  403. {
  404.     // We do nothing as we say it's a lazy-learning method
  405.     // Only when predict() is called, the learning process is invoked
  406.     memset( ptr->lineParam, 0, sizeof(double)*MAX_POL_DEPTH );
  407.     ptr->type = KNNModel;
  408.     return true;
  409. }

Util.cpp

  1. #include "Revise.h"
  3. int _tmain(int argc, _TCHAR* argv[])
  4. {
  5.     double x[10], y[10];
  6.     for ( int i = 0; i < 10; i ++ )
  7.     {
  8.         x[i] = 10 - i;
  9.         y[i] = (i-2)*(i-2);
  10.     }
  12.     Model *model = CreateModel();
  13.     SetOptData( model, x, y, 10 );
  14.     Opt( model, BrokenLine );
  15.     double result = Predict( model, 1.5 );
  16.     ReleaseModel( &model );
  18.     return 0;
  19. }

发表评论

表情：
评论列表 （有 0 条评论，323人围观）

还没有评论，来说两句吧...

相关阅读

    相关 拟合和抵抗过拟合

    回归拟合有三种情况： ![在这里插入图片描述][20210221113528609.png] (1)欠拟合就是模型不能正确预测出数据的分布情况。 (2)正确拟合就是

    朱雀朱雀/ 0 赞/ 266 阅读

    相关 过拟合和欠拟合

    开始我是很难弄懂什么是过拟合，什么是欠拟合以及造成两者的各自原因以及相应的解决办法，学习了一段时间机器学习和深度学习后，分享下自己的观点，方便初学者能很好很形象地理解上面的问题

    朴灿烈づ我的快乐病毒、朴灿烈づ我的快乐病毒、/ 0 赞/ 340 阅读

    相关 MATLAB 拟合曲线

    最近做一个实验 ，代码中需要找出数据之间的函数关系，所以需要把数据进行拟合，找到关系式。听别人说MATLAB做 拟合的时候特别方便，所以就测试了一下，果真方便，现在将大

    野性酷女野性酷女/ 0 赞/ 318 阅读